Mercurial > repos > davidvanzessen > imgt_loader_igg

comparison imgtconvert.py @ 7:04e72fc8b2c4 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Fri, 05 Sep 2014 04:21:48 -0400
parents 5b030e48b308
children
comparison
equal deleted inserted replaced
6:5b030e48b308 7:04e72fc8b2c4
1 import pandas as pd
2 try:
3 pd.options.mode.chained_assignment = None # default='warn'
4 except:
5 pass
6 import re
7 import argparse
8 import os
9
10 def stop_err( msg, ret=1 ):
11 sys.stderr.write( msg )
12 sys.exit( ret )
13
14 #docs.python.org/dev/library/argparse.html
15 parser = argparse.ArgumentParser()
16 parser.add_argument("--input", help="Input folder with files")
17 parser.add_argument("--output", help="Output file")
18
19 args = parser.parse_args()
20
21 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']
22 old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT']
23 old_junction_columns = [u'JUNCTION']
24
25 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']
26 added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT']
27 added_junction_columns = [u"P3'V-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", u"5'D-REGION trimmed-nt nb", u"3'D-REGION trimmed-nt nb", u"5'J-REGION trimmed-nt nb"]
28
29 inputFolder = args.input
30
31 dirContents = os.listdir(inputFolder)
32 if len(dirContents) == 1:
33 inputFolder = os.path.join(inputFolder, dirContents[0])
34 if os.path.isdir(inputFolder):
35 dirContents = os.listdir(inputFolder)
36 files = sorted([os.path.join(inputFolder, f) for f in dirContents if os.path.isfile(os.path.join(inputFolder, f))])
37
38 if len(files) % 3 is not 0:
39 stop_err("Files in zip not a multiple of 3, it should contain the all the 1_, 5_ and 6_ files for a sample")
40 import sys
41 sys.exit()
42
43 triplets = []
44 step = len(files) / 3
45 for i in range(0, step):
46 triplets.append((files[i], files[i + step], files[i + step + step]))
47
48 outFile = args.output
49
50 #fSummary = pd.read_csv(triplets[0][0], sep="\t", low_memory=False)
51 fSummary = pd.read_csv(triplets[0][0], sep="\t")
52 #fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False)
53 fSequence = pd.read_csv(triplets[0][1], sep="\t")
54 #fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False)
55 fJunction = pd.read_csv(triplets[0][2], sep="\t")
56 tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
57
58 tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
59 tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
60
61 tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
62 tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
63
64 tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
65 tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
66
67 tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
68 tmp["CDR3 Length DNA"] = '1'
69 tmp["Strand"] = fSummary["Orientation"]
70 tmp["CDR3 Found How"] = 'a'
71
72 for col in added_summary_columns:
73 tmp[col] = fSummary[col]
74
75 for col in added_sequence_columns:
76 tmp[col] = fSequence[col]
77
78 for col in added_junction_columns:
79 tmp[col] = fJunction[col]
80
81 outFrame = tmp
82
83
84
85 for triple in triplets[1:]:
86 fSummary = pd.read_csv(triple[0], sep="\t")
87 fSequence = pd.read_csv(triple[1], sep="\t")
88 fJunction = pd.read_csv(triple[2], sep="\t")
89
90 tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
91
92 tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
93 tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
94
95 tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
96 tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
97
98 tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
99 tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
100
101 tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
102 tmp["CDR3 Length DNA"] = '1'
103 tmp["Strand"] = fSummary["Orientation"]
104 tmp["CDR3 Found How"] = 'a'
105
106 for col in added_summary_columns:
107 tmp[col] = fSummary[col]
108
109 for col in added_sequence_columns:
110 tmp[col] = fSequence[col]
111
112 for col in added_junction_columns:
113 tmp[col] = fJunction[col]
114
115 outFrame = outFrame.append(tmp)
116
117
118 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', 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', 'AA JUNCTION', 'Functionality comment', 'Sequence', 'FR1-IMGT', 'FR2-IMGT', 'FR3-IMGT', 'CDR3-IMGT', 'JUNCTION', 'J-REGION', 'FR4-IMGT', 'P3V-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', '5J-REGION trimmed-nt nb']
119
120 """
121 IGHV[0-9]-[0-9ab]+-?[0-9]?D?
122 TRBV[0-9]{1,2}-?[0-9]?-?[123]?
123 IGKV[0-3]D?-[0-9]{1,2}
124 IGLV[0-9]-[0-9]{1,2}
125 TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?
126 TRGV[234589]
127 TRDV[1-3]
128
129 IGHD[0-9]-[0-9ab]+
130 TRBD[12]
131 TRDD[1-3]
132
133 IGHJ[1-6]
134 TRBJ[12]-[1-7]
135 IGKJ[1-5]
136 IGLJ[12367]
137 TRAJ[0-9]{1,2}
138 TRGJP?[12]
139 TRDJ[1-4]
140 """
141
142 vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)",
143 r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)",
144 r"(IGKV[0-3]D?-[0-9]{1,2})",
145 r"(IGLV[0-9]-[0-9]{1,2})",
146 r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)",
147 r"(TRGV[234589])",
148 r"(TRDV[1-3])"]
149
150 dPattern = [r"(IGHD[0-9]-[0-9ab]+)",
151 r"(TRBD[12])",
152 r"(TRDD[1-3])"]
153
154 jPattern = [r"(IGHJ[1-6])",
155 r"(TRBJ[12]-[1-7])",
156 r"(IGKJ[1-5])",
157 r"(IGLJ[12367])",
158 r"(TRAJ[0-9]{1,2})",
159 r"(TRGJP?[12])",
160 r"(TRDJ[1-4])"]
161
162 vPattern = re.compile(r"|".join(vPattern))
163
164 dPattern = re.compile(r"|".join(dPattern))
165
166 jPattern = re.compile(r"|".join(jPattern))
167
168
169 def filterGenes(s, pattern):
170 if type(s) is not str:
171 return "NA"
172 res = pattern.search(s)
173 if res:
174 return res.group(0)
175 return "NA"
176
177
178
179 outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern))
180 outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern))
181 outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern))
182
183 print outFrame
184
185 tmp = outFrame["VDJ Frame"]
186 tmp = tmp.replace("in-frame", "In-frame")
187 tmp = tmp.replace("null", "Out-of-frame")
188 tmp = tmp.replace("out-of-frame", "Out-of-frame")
189 outFrame["VDJ Frame"] = tmp
190 outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len)
191 safeLength = lambda x: len(x) if type(x) == str else 0
192 outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows?
193 #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?
194 outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")