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