Mercurial > repos > davidvanzessen > argalaxy_tools

comparison gene_identification.py @ 11:0510cf1f7cbc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 04 Aug 2015 09:59:26 -0400
parents
children
comparison
equal deleted inserted replaced
10:edbf4fba5fc7 11:0510cf1f7cbc
1 import re
2 import argparse
3 import time
4 starttime= int(time.time() * 1000)
5
6 parser = argparse.ArgumentParser()
7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file")
8 parser.add_argument("--output", help="The annotated output file to be merged back with the summary file")
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 output = args.output
15 #outfile = "identified.txt"
16
17 dic = dict()
18 total = 0
19
20
21 first = True
22 IDIndex = 0
23 seqIndex = 0
24
25 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
26 for line in f:
27 total += 1
28 if first:
29 linesplt = line.split("\t")
30 IDIndex = linesplt.index("Sequence ID")
31 seqIndex = linesplt.index("Sequence")
32 first = False
33 continue
34 linesplt = line.split("\t")
35 ID = linesplt[IDIndex]
36 if len(linesplt) < 28: #weird rows without a sequence
37 dic[ID] = ""
38 else:
39 dic[ID] = linesplt[seqIndex]
40
41 print "Number of input sequences:", len(dic)
42
43 #lambda/kappa reference sequence
44 searchstrings = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccagatgggaacgtggtcatcgcctgcctgg",
45 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag",
46 "cm": "gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc"}
47
48 compiledregex = {"ca": [],
49 "cg": [],
50 "cm": []}
51
52 #lambda/kappa reference sequence variable nucleotides
53 ca1 = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'}
54 ca2 = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'}
55 cg1 = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
56 cg2 = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'}
57 cg3 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'}
58 cg4 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'}
59
60 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap
61 chunklength = 8
62
63 #create the chunks of the reference sequence with regular expressions for the variable nucleotides
64 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength / 2):
65 pos = i
66 chunk = searchstrings["ca"][i:i+chunklength]
67 result = ""
68 varsInResult = 0
69 for c in chunk:
70 if pos in ca1.keys():
71 varsInResult += 1
72 result += "[" + ca1[pos] + ca2[pos] + "]"
73 else:
74 result += c
75 pos += 1
76 compiledregex["ca"].append((re.compile(result), varsInResult))
77
78 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength / 2):
79 pos = i
80 chunk = searchstrings["cg"][i:i+chunklength]
81 result = ""
82 varsInResult = 0
83 for c in chunk:
84 if pos in cg1.keys():
85 varsInResult += 1
86 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]"
87 else:
88 result += c
89 pos += 1
90 compiledregex["cg"].append((re.compile(result), varsInResult))
91
92 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength / 2):
93 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False))
94
95
96
97 def removeAndReturnMaxIndex(x): #simplifies a list comprehension
98 m = max(x)
99 index = x.index(m)
100 x[index] = 0
101 return index
102
103
104 start_location = dict()
105 hits = dict()
106 alltotal = 0
107 for key in compiledregex.keys(): #for ca/cg/cm
108 regularexpressions = compiledregex[key] #get the compiled regular expressions
109 for ID in dic.keys()[0:]: #for every ID
110 if ID not in hits.keys(): #ensure that the dictionairy that keeps track of the hits for every gene exists
111 hits[ID] = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0}
112 currentIDHits = hits[ID]
113 seq = dic[ID]
114 lastindex = 0
115 start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0)
116 start = [0] * (len(seq) + start_zero)
117 for i, regexp in enumerate(regularexpressions): #for every regular expression
118 relativeStartLocation = lastindex - (chunklength / 2) * i
119 if relativeStartLocation >= len(seq):
120 break
121 regex, hasVar = regexp
122 matches = regex.finditer(seq[lastindex:])
123 for match in matches: #for every match with the current regex, only uses the first hit
124 lastindex += match.start()
125 start[relativeStartLocation + start_zero] += 1
126 if hasVar: #if the regex has a variable nt in it
127 chunkstart = chunklength / 2 * i #where in the reference does this chunk start
128 chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end
129 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
130 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]])
131 currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]])
132 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper
133 currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]])
134 currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]])
135 currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]])
136 currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]])
137 else: #key == "cm" #no variable regions in 'cm'
138 pass
139 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped
140 else: #only runs if there were no hits
141 continue
142 #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]
143 currentIDHits[key + "_hits"] += 1
144 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1])
145 #start_location[ID + "_" + key] = str(start.index(max(start)))
146
147
148 chunksInCA = len(compiledregex["ca"])
149 chunksInCG = len(compiledregex["cg"])
150 chunksInCM = len(compiledregex["cm"])
151 requiredChunkPercentage = 0.7
152 varsInCA = float(len(ca1.keys()) * 2)
153 varsInCG = float(len(cg1.keys()) * 2) + 1
154 varsInCM = 0
155 requiredVarPercentage = 0.7
156
157
158 first = True
159 seq_write_count=0
160 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
161 with open(output, 'w') as o:
162 for line in f:
163 total += 1
164 if first:
165 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
166 first = False
167 continue
168 linesplt = line.split("\t")
169 if linesplt[2] == "No results":
170 pass
171 ID = linesplt[1]
172 currentIDHits = hits[ID]
173 possibleca = float(len(compiledregex["ca"]))
174 possiblecg = float(len(compiledregex["cg"]))
175 possiblecm = float(len(compiledregex["cm"]))
176 cahits = currentIDHits["ca_hits"]
177 cghits = currentIDHits["cg_hits"]
178 cmhits = currentIDHits["cm_hits"]
179 if cahits >= cghits and cahits >= cmhits: #its a ca gene
180 ca1hits = currentIDHits["ca1"]
181 ca2hits = currentIDHits["ca2"]
182 if ca1hits >= ca2hits:
183 o.write(ID + "\tca1\t" + str(int(ca1hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
184 else:
185 o.write(ID + "\tca2\t" + str(int(ca2hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
186 elif cghits >= cahits and cghits >= cmhits: #its a cg gene
187 cg1hits = currentIDHits["cg1"]
188 cg2hits = currentIDHits["cg2"]
189 cg3hits = currentIDHits["cg3"]
190 cg4hits = currentIDHits["cg4"]
191 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene
192 o.write(ID + "\tcg1\t" + str(int(cg1hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
193 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene
194 o.write(ID + "\tcg2\t" + str(int(cg2hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
195 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene
196 o.write(ID + "\tcg3\t" + str(int(cg3hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
197 else: #cg4 gene
198 o.write(ID + "\tcg4\t" + str(int(cg4hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
199 else: #its a cm gene
200 o.write(ID + "\tcm\t0\t" + str(int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
201 seq_write_count += 1
202
203 print "Time: %i" % (int(time.time() * 1000) - starttime)
204
205 print "Number of sequences written to file:", seq_write_count
206
207
208
209
210