annotate imgt_loader/imgt_loader.py @ 13:d3ebaa2d2fe0 draft

Uploaded
author davidvanzessen
date Tue, 20 Dec 2016 06:02:44 -0500
parents afe85eb6572e
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
1 import pandas as pd
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
2 try:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
3 pd.options.mode.chained_assignment = None # default='warn'
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
4 except:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
5 pass
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
6 import re
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
7 import argparse
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
8 import os
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
9
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
10 def stop_err( msg, ret=1 ):
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
11 sys.stderr.write( msg )
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
12 sys.exit( ret )
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
13
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
14 #docs.python.org/dev/library/argparse.html
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
15 parser = argparse.ArgumentParser()
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
16 parser.add_argument("--summ", help="The 1_Summary file from the imgt output")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
17 parser.add_argument("--aa", help="The 5_AA-Sequence file from the imgt output")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
18 parser.add_argument("--junction", help="The 6_Junction file from the imgt output")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
19 parser.add_argument("--output", help="Output file")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
20
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
21 args = parser.parse_args()
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
22
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
23 old_summary_columns = [u'Sequence ID', u'JUNCTION frame', u'V-GENE and allele', u'D-GENE and allele', u'J-GENE and allele', u'CDR1-IMGT length', u'CDR2-IMGT length', u'CDR3-IMGT length', u'Orientation']
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
24 old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT']
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
25 old_junction_columns = [u'JUNCTION']
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
26
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
27 added_summary_columns = [u'Functionality', u'V-REGION identity %', u'V-REGION identity nt', u'D-REGION reading frame', u'AA JUNCTION', u'Functionality comment', u'Sequence']
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
28 added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT']
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
29 added_junction_columns = [u"P3'V-nt nb", u'N-REGION-nt nb', u'N1-REGION-nt nb', u"P5'D-nt nb", u"P3'D-nt nb", u'N2-REGION-nt nb', u"P5'J-nt nb", u"3'V-REGION trimmed-nt nb",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
30 u"5'D-REGION trimmed-nt nb", u"3'D-REGION trimmed-nt nb", u"5'J-REGION trimmed-nt nb", u"N-REGION", u"N1-REGION", u"N2-REGION"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
31
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
32 outFile = args.output
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
33
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
34 #fSummary = pd.read_csv(triplets[0][0], sep="\t", low_memory=False)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
35 fSummary = pd.read_csv(args.summ, sep="\t", dtype=object)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
36 #fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
37 fSequence = pd.read_csv(args.aa, sep="\t", dtype=object)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
38 #fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
39 fJunction = pd.read_csv(args.junction, sep="\t", dtype=object)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
40 tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
41
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
42 tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
43 tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
44
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
45 tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
46 tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
47
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
48 tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
49 tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
50
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
51 tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
52 tmp["CDR3 Length DNA"] = '1'
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
53 tmp["Strand"] = fSummary["Orientation"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
54 tmp["CDR3 Found How"] = 'a'
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
55
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
56 for col in added_summary_columns:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
57 tmp[col] = fSummary[col]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
58
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
59 for col in added_sequence_columns:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
60 tmp[col] = fSequence[col]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
61
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
62 for col in added_junction_columns:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
63 tmp[col] = fJunction[col]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
64
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
65 outFrame = tmp
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
66
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
67 outFrame.columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', u'CDR3 Seq', u'CDR3 Length',
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
68 u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'Functionality', 'V-REGION identity %', 'V-REGION identity nt', 'D-REGION reading frame',
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
69 'AA JUNCTION', 'Functionality comment', 'Sequence', 'FR1-IMGT', 'FR2-IMGT', 'FR3-IMGT', 'CDR3-IMGT', 'JUNCTION', 'J-REGION', 'FR4-IMGT', 'P3V-nt nb',
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
70 'N-REGION-nt nb', 'N1-REGION-nt nb', 'P5D-nt nb', 'P3D-nt nb', 'N2-REGION-nt nb', 'P5J-nt nb', '3V-REGION trimmed-nt nb', '5D-REGION trimmed-nt nb', '3D-REGION trimmed-nt nb',
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
71 '5J-REGION trimmed-nt nb', "N-REGION", "N1-REGION", "N2-REGION"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
72
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
73 """
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
74 IGHV[0-9]-[0-9ab]+-?[0-9]?D?
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
75 TRBV[0-9]{1,2}-?[0-9]?-?[123]?
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
76 IGKV[0-3]D?-[0-9]{1,2}
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
77 IGLV[0-9]-[0-9]{1,2}
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
78 TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
79 TRGV[234589]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
80 TRDV[1-3]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
81
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
82 IGHD[0-9]-[0-9ab]+
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
83 TRBD[12]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
84 TRDD[1-3]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
85
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
86 IGHJ[1-6]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
87 TRBJ[12]-[1-7]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
88 IGKJ[1-5]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
89 IGLJ[12367]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
90 TRAJ[0-9]{1,2}
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
91 TRGJP?[12]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
92 TRDJ[1-4]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
93 """
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
94
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
95 vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
96 r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
97 r"(IGKV[0-3]D?-[0-9]{1,2})",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
98 r"(IGLV[0-9]-[0-9]{1,2})",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
99 r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
100 r"(TRGV[234589])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
101 r"(TRDV[1-3])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
102 r"(IGHV[0-9]S[0-9]+)"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
103
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
104 dPattern = [r"(IGHD[0-9]-[0-9ab]+)",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
105 r"(TRBD[12])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
106 r"(TRDD[1-3])"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
107
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
108 jPattern = [r"(IGHJ[1-6])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
109 r"(TRBJ[12]-[1-7])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
110 r"(IGKJ[1-5])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
111 r"(IGLJ[12367])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
112 r"(TRAJ[0-9]{1,2})",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
113 r"(TRGJP?[12])",
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
114 r"(TRDJ[1-4])"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
115
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
116 vPattern = re.compile(r"|".join(vPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
117
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
118 dPattern = re.compile(r"|".join(dPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
119
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
120 jPattern = re.compile(r"|".join(jPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
121
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
122
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
123 def filterGenes(s, pattern):
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
124 if type(s) is not str:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
125 return "NA"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
126 res = pattern.search(s)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
127 if res:
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
128 return res.group(0)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
129 return "NA"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
130
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
131
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
132
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
133 outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
134 outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
135 outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern))
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
136
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
137
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
138 tmp = outFrame["VDJ Frame"]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
139 tmp = tmp.replace("in-frame", "In-frame")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
140 tmp = tmp.replace("null", "Out-of-frame")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
141 tmp = tmp.replace("out-of-frame", "Out-of-frame")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
142 outFrame["VDJ Frame"] = tmp
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
143 outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
144 safeLength = lambda x: len(x) if type(x) == str else 0
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
145 #outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows?
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
146 #outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows?
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
147 outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")