|
0
|
1 import re
|
|
|
2 import argparse
|
|
2
|
3 import time
|
|
|
4 starttime= int(time.time() * 1000)
|
|
0
|
5
|
|
|
6 parser = argparse.ArgumentParser()
|
|
|
7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file")
|
|
|
8 parser.add_argument("--outdir", help="Output directory, 7 output files will be written here")
|
|
|
9
|
|
|
10 args = parser.parse_args()
|
|
|
11
|
|
|
12 infile = args.input
|
|
|
13 #infile = "test_VH-Ca_Cg_25nt/1_Summary_test_VH-Ca_Cg_25nt_241013.txt"
|
|
|
14 outdir = args.outdir
|
|
|
15 #outfile = "identified.txt"
|
|
|
16
|
|
|
17 dic = dict()
|
|
|
18 total = 0
|
|
|
19
|
|
|
20 first = True
|
|
|
21 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
|
|
|
22 for line in f:
|
|
|
23 total += 1
|
|
|
24 if first:
|
|
|
25 first = False
|
|
|
26 continue
|
|
|
27 linesplt = line.split("\t")
|
|
|
28 if linesplt[2] == "No results":
|
|
|
29 continue
|
|
|
30 ID = linesplt[1]
|
|
|
31 seq = linesplt[28]
|
|
|
32 dic[ID] = seq
|
|
|
33
|
|
|
34 #lambda/kappa reference sequence
|
|
|
35 searchstrings = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccagatgggaacgtggtcatcgcctgcctggtccagggcttcttcccccaggagccactcagtgtgacctggagcgaaag",
|
|
|
36 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccagcggcgtgcacaccttcc",
|
|
|
37 "cm": "gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc"}
|
|
|
38
|
|
|
39 compiledregex = {"ca": [],
|
|
|
40 "cg": [],
|
|
|
41 "cm": []}
|
|
|
42
|
|
|
43 #lambda/kappa reference sequence variable nucleotides
|
|
|
44 ca1 = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'}
|
|
|
45 ca2 = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'}
|
|
|
46 cg1 = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
|
|
|
47 cg2 = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'}
|
|
|
48 cg3 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
|
|
|
49 cg4 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'}
|
|
|
50
|
|
|
51 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap
|
|
|
52 chunklength = 8
|
|
|
53
|
|
|
54 #create the chunks of the reference sequence with regular expressions for the variable nucleotides
|
|
|
55 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength / 2):
|
|
|
56 pos = i
|
|
|
57 chunk = searchstrings["ca"][i:i+chunklength]
|
|
|
58 result = ""
|
|
|
59 varsInResult = 0
|
|
|
60 for c in chunk:
|
|
|
61 if pos in ca1.keys():
|
|
|
62 varsInResult += 1
|
|
|
63 result += "[" + ca1[pos] + ca2[pos] + "]"
|
|
|
64 else:
|
|
|
65 result += c
|
|
|
66 pos += 1
|
|
|
67 compiledregex["ca"].append((re.compile(result), varsInResult))
|
|
|
68
|
|
|
69 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength / 2):
|
|
|
70 pos = i
|
|
|
71 chunk = searchstrings["cg"][i:i+chunklength]
|
|
|
72 result = ""
|
|
|
73 varsInResult = 0
|
|
|
74 for c in chunk:
|
|
|
75 if pos in cg1.keys():
|
|
|
76 varsInResult += 1
|
|
|
77 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]"
|
|
|
78 else:
|
|
|
79 result += c
|
|
|
80 pos += 1
|
|
|
81 compiledregex["cg"].append((re.compile(result), varsInResult))
|
|
|
82
|
|
|
83 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength / 2):
|
|
|
84 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False))
|
|
|
85
|
|
|
86
|
|
|
87
|
|
|
88 def removeAndReturnMaxIndex(x): #simplifies a list comprehension
|
|
|
89 m = max(x)
|
|
|
90 index = x.index(m)
|
|
|
91 x[index] = 0
|
|
|
92 return index
|
|
|
93
|
|
|
94
|
|
|
95 start_location = dict()
|
|
|
96 hits = dict()
|
|
|
97 alltotal = 0
|
|
|
98 for key in compiledregex.keys(): #for ca/cg/cm
|
|
|
99 regularexpressions = compiledregex[key] #get the compiled regular expressions
|
|
|
100 for ID in dic.keys()[0:]: #for every ID
|
|
|
101 if ID not in hits.keys(): #ensure that the dictionairy that keeps track of the hits for every gene exists
|
|
|
102 hits[ID] = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0}
|
|
|
103 currentIDHits = hits[ID]
|
|
|
104 seq = dic[ID]
|
|
|
105 lastindex = 0
|
|
|
106 start = [0] * len(seq)
|
|
|
107 for i, regexp in enumerate(regularexpressions): #for every regular expression
|
|
|
108 regex, hasVar = regexp
|
|
|
109 matches = regex.finditer(seq[lastindex:])
|
|
|
110 for match in matches: #for every match with the current regex, only uses the first hit
|
|
|
111 lastindex += match.start()
|
|
|
112 start[lastindex - chunklength / 2 * i] += 1
|
|
|
113 if hasVar: #if the regex has a variable nt in it
|
|
|
114 chunkstart = chunklength / 2 * i #where in the reference does this chunk start
|
|
|
115 chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end
|
|
|
116 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
|
|
|
117 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]])
|
|
|
118 currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]])
|
|
|
119 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper
|
|
|
120 currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]])
|
|
|
121 currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]])
|
|
|
122 currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]])
|
|
|
123 currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]])
|
|
|
124 else: #key == "cm" #no variable regions in 'cm'
|
|
|
125 pass
|
|
|
126 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped
|
|
|
127 else: #only runs if there were no hits
|
|
|
128 continue
|
|
|
129 #print "found ", regex.pattern , "at", lastindex, "adding one to", (lastindex - chunklength / 2 * i), "to the start array of", ID, "gene", key, "it's now:", start[lastindex - chunklength / 2 * i]
|
|
|
130 currentIDHits[key + "_hits"] += 1
|
|
|
131 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1) for x in range(5) if max(start) > 1])
|
|
|
132 #start_location[ID + "_" + key] = str(start.index(max(start)))
|
|
|
133
|
|
|
134
|
|
|
135 chunksInCA = len(compiledregex["ca"])
|
|
|
136 chunksInCG = len(compiledregex["cg"])
|
|
|
137 chunksInCM = len(compiledregex["cm"])
|
|
|
138 requiredChunkPercentage = 0.7
|
|
|
139 varsInCA = float(len(ca1.keys()) * 2)
|
|
|
140 varsInCG = float(len(cg1.keys()) * 2) + 1
|
|
|
141 varsInCM = 0
|
|
|
142 requiredVarPercentage = 0.7
|
|
|
143
|
|
|
144 ca = 0
|
|
|
145 ca1 = 0
|
|
|
146 ca2 = 0
|
|
|
147 cg = 0
|
|
|
148 cg1 = 0
|
|
|
149 cg2 = 0
|
|
|
150 cg3 = 0
|
|
|
151 cg4 = 0
|
|
|
152 cm = 0
|
|
|
153 try:
|
|
|
154 cafile = open(outdir + "/ca.txt", 'w')
|
|
|
155 ca1file = open(outdir + "/ca1.txt", 'w')
|
|
|
156 ca2file = open(outdir + "/ca2.txt", 'w')
|
|
|
157 cgfile = open(outdir + "/cg.txt", 'w')
|
|
|
158 cg1file = open(outdir + "/cg1.txt", 'w')
|
|
|
159 cg2file = open(outdir + "/cg2.txt", 'w')
|
|
|
160 cg3file = open(outdir + "/cg3.txt", 'w')
|
|
|
161 cg4file = open(outdir + "/cg4.txt", 'w')
|
|
|
162 cmfile = open(outdir + "/cm.txt", 'w')
|
|
|
163 unmatchedfile = open(outdir + "/unmatched.txt", 'w')
|
|
|
164 cafile.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
165 ca1file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
166 ca2file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
167 cgfile.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
168 cg1file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
169 cg2file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
170 cg3file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
171 cg4file.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
172 cmfile.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
|
|
|
173 unmatchedfile.write("ID\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\tbest_match\n")
|
|
|
174 for ID in hits.keys():
|
|
|
175 currentIDHits = hits[ID]
|
|
|
176 possibleca = float(len(compiledregex["ca"]))
|
|
|
177 possiblecg = float(len(compiledregex["cg"]))
|
|
|
178 possiblecm = float(len(compiledregex["cm"]))
|
|
|
179 cahits = currentIDHits["ca_hits"]
|
|
|
180 cghits = currentIDHits["cg_hits"]
|
|
|
181 cmhits = currentIDHits["cm_hits"]
|
|
|
182 if cahits > cghits and cahits > cmhits: #its a ca gene
|
|
|
183 if cahits <= int(chunksInCA * requiredChunkPercentage):
|
|
|
184 unmatchedfile.write(ID + "\tNA\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\tca\n")
|
|
|
185 continue
|
|
|
186 ca += 1
|
|
|
187 ca1hits = currentIDHits["ca1"]
|
|
|
188 ca2hits = currentIDHits["ca2"]
|
|
|
189 cafile.write(ID + "\tNA\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
|
|
|
190 if ca1hits > ca2hits:
|
|
|
191 #print ID, "is ca1 with", (ca1hits / 2), "hits for ca1 and", (ca2hits / 2), "hits for ca2", (int((ca1hits / varsInCA) * 100)), "percent hit"
|
|
|
192 if ca1hits <= int(varsInCA * requiredVarPercentage):
|
|
|
193 unmatchedfile.write(ID + "\t" + str(int(ca1hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\tca1\n")
|
|
|
194 continue
|
|
|
195 ca1 += 1
|
|
|
196 ca1file.write(ID + "\t" + str(int(ca1hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
|
|
|
197 else:
|
|
|
198 #print ID, "is ca2 with", (ca1hits / 2), "hits for ca1 and", (ca2hits / 2), "hits for ca2", (int((ca2hits / varsInCA) * 100)), "percent hit"
|
|
|
199 if ca2hits <= int(varsInCA * requiredVarPercentage):
|
|
|
200 unmatchedfile.write(ID + "\t" + str(int(ca2hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\tca1\n")
|
|
|
201 continue
|
|
|
202 ca2 += 1
|
|
|
203 ca2file.write(ID + "\t" + str(int(ca2hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
|
|
|
204 elif cghits > cahits and cghits > cmhits: #its a cg gene
|
|
|
205 if cghits <= int(chunksInCG * requiredChunkPercentage):
|
|
|
206 unmatchedfile.write(ID + "\tNA\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_ca"] + "\tcg\n")
|
|
|
207 continue
|
|
|
208 cg += 1
|
|
|
209 cg1hits = currentIDHits["cg1"]
|
|
|
210 cg2hits = currentIDHits["cg2"]
|
|
|
211 cg3hits = currentIDHits["cg3"]
|
|
|
212 cg4hits = currentIDHits["cg4"]
|
|
|
213 cgfile.write(ID + "\tNA\t" + str(int(cghits / possibleca * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
|
|
|
214 if cg1hits > cg2hits and cg1hits > cg3hits and cg1hits > cg4hits: #cg1 gene
|
|
|
215 if cg1hits <= int(varsInCG * requiredVarPercentage):
|
|
|
216 unmatchedfile.write(ID + "\t" + str(int(cg1hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\tcg1\n")
|
|
|
217 continue
|
|
|
218 cg1 += 1
|
|
|
219 cg1file.write(ID + "\t" + str(int(cg1hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
|
|
|
220 elif cg2hits > cg1hits and cg2hits > cg3hits and cg2hits > cg4hits: #cg2 gene
|
|
|
221 if cg2hits <= int(varsInCG * requiredVarPercentage):
|
|
|
222 unmatchedfile.write(ID + "\t" + str(int(cg2hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\tcg2\n")
|
|
|
223 continue
|
|
|
224 cg2 += 1
|
|
|
225 cg2file.write(ID + "\t" + str(int(cg2hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
|
|
|
226 elif cg3hits > cg1hits and cg3hits > cg2hits and cg3hits > cg4hits: #cg3 gene
|
|
|
227 if cg3hits <= int(varsInCG * requiredVarPercentage):
|
|
|
228 unmatchedfile.write(ID + "\t" + str(int(cg3hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\tcg3\n")
|
|
|
229 continue
|
|
|
230 cg3 += 1
|
|
|
231 cg3file.write(ID + "\t" + str(int(cg3hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
|
|
|
232 else: #cg4 gene
|
|
|
233 if cg4hits <= int(varsInCG * requiredVarPercentage):
|
|
|
234 unmatchedfile.write(ID + "\t" + str(int(cg4hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\tcg4\n")
|
|
|
235 continue
|
|
|
236 cg4 += 1
|
|
|
237 cg4file.write(ID + "\t" + str(int(cg4hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
|
|
|
238 else: #its a cm gene
|
|
|
239 if cmhits <= int(chunksInCM * requiredChunkPercentage):
|
|
|
240 unmatchedfile.write(ID + "\tNA\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_ca"] + "\tcm\n")
|
|
|
241 continue
|
|
|
242 cm += 1
|
|
|
243 cmfile.write(ID + "\tNA\t" + str(int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n")
|
|
|
244 finally:
|
|
|
245 cafile.close()
|
|
|
246 ca1file.close()
|
|
|
247 ca2file.close()
|
|
|
248 cgfile.close()
|
|
|
249 cg1file.close()
|
|
|
250 cg2file.close()
|
|
|
251 cg3file.close()
|
|
|
252 cg4file.close()
|
|
|
253 cmfile.close()
|
|
|
254 unmatchedfile.close()
|
|
|
255
|
|
|
256
|
|
|
257 #print ca,cg,cm,(ca+cg+cm)
|
|
|
258
|
|
2
|
259 print "Time: %i" % (int(time.time() * 1000) - starttime)
|
|
0
|
260
|
|
|
261
|
|
|
262
|
|
|
263
|
|
|
264
|
|
|
265
|
|
|
266
|
