Mercurial > repos > davidvanzessen > imgt_loader_igg

comparison imgtconvert.py @ 0:f8ba6ab5514b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 23 Jan 2014 08:09:30 -0500
parents
children 4405849aa053
comparison
equal deleted inserted replaced
-1:000000000000 0:f8ba6ab5514b
1 import pandas as pd
2 import re
3 import argparse
4 import os
5
6 def stop_err( msg, ret=1 ):
7 sys.stderr.write( msg )
8 sys.exit( ret )
9
10 #docs.python.org/dev/library/argparse.html
11 parser = argparse.ArgumentParser()
12 parser.add_argument("--input", help="Input folder with files")
13 parser.add_argument("--output", help="Output file")
14
15 args = parser.parse_args()
16
17 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']
18 old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT']
19 old_junction_columns = [u'JUNCTION']
20
21 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']
22 added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT']
23 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"]
24
25 inputFolder = args.input
26
27 dirContents = os.listdir(inputFolder)
28 if len(dirContents) == 1:
29 inputFolder = os.path.join(inputFolder, dirContents[0])
30 if os.path.isdir(inputFolder):
31 print "is dir"
32 dirContents = os.listdir(inputFolder)
33 files = sorted([os.path.join(inputFolder, f) for f in dirContents])
34
35 if len(files) % 3 is not 0:
36 stop_err("Files in zip not a multiple of 3, it should contain the all the 1_, 5_ and 6_ files for a sample")
37 import sys
38 sys.exit()
39
40 triplets = []
41 step = len(files) / 3
42 for i in range(0, step):
43 triplets.append((files[i], files[i + step], files[i + step + step]))
44
45 outFile = args.output
46
47 fSummary = pd.read_csv(triplets[0][0], sep="\t")
48 fSequence = pd.read_csv(triplets[0][1], sep="\t")
49 fJunction = pd.read_csv(triplets[0][2], sep="\t")
50 tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
51
52 tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
53 tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
54
55 tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
56 tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
57
58 tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
59 tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
60
61 tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
62 tmp["CDR3 Length DNA"] = '1'
63 tmp["Strand"] = fSummary["Orientation"]
64 tmp["CDR3 Found How"] = 'a'
65
66 for col in added_summary_columns:
67 tmp[col] = fSummary[col]
68
69 for col in added_sequence_columns:
70 tmp[col] = fSequence[col]
71
72 for col in added_junction_columns:
73 tmp[col] = fJunction[col]
74
75 outFrame = tmp
76
77 for triple in triplets[1:]:
78 fSummary = pd.read_csv(triple[0], sep="\t")
79 fSequence = pd.read_csv(triple[1], sep="\t")
80 fJunction = pd.read_csv(triple[2], sep="\t")
81
82 tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
83
84 tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
85 tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
86
87 tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
88 tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
89
90 tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
91 tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
92
93 tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
94 tmp["CDR3 Length DNA"] = '1'
95 tmp["Strand"] = fSummary["Orientation"]
96 tmp["CDR3 Found How"] = 'a'
97
98 for col in added_summary_columns:
99 tmp[col] = fSummary[col]
100
101 for col in added_sequence_columns:
102 tmp[col] = fSequence[col]
103
104 for col in added_junction_columns:
105 tmp[col] = fJunction[col]
106
107 outFrame = outFrame.append(tmp)
108
109 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']
110
111 vPattern = re.compile(r"IGHV[1-9]-[0-9ab]+-?[1-9]?")
112 dPattern = re.compile(r"IGHD[1-9]-[0-9ab]+")
113 jPattern = re.compile(r"IGHJ[1-9]")
114
115 def filterGenes(s, pattern):
116 if type(s) is not str:
117 return "NA"
118 res = pattern.search(s)
119 if res:
120 return res.group(0)
121 return "NA"
122
123
124 outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern))
125 outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern))
126 outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern))
127
128
129
130 tmp = outFrame["VDJ Frame"]
131 tmp = tmp.replace("in-frame", "In-frame")
132 tmp = tmp.replace("null", "Out-of-frame")
133 tmp = tmp.replace("out-of-frame", "Out-of-frame")
134 outFrame["VDJ Frame"] = tmp
135 outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len)
136 safeLength = lambda x: len(x) if type(x) == str else 0
137 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?
138 outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")