Mercurial > repos > davidvanzessen > mutation_analysis

comparison gene_identification.py @ 19:c518cf0d4adb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Wed, 01 Apr 2015 05:09:59 -0400
parents b0242cd1da34
children c9f9623f1f76
comparison
equal deleted inserted replaced
18:9dd557f17153 19:c518cf0d4adb
3 import time 3 import time
4 starttime= int(time.time() * 1000) 4 starttime= int(time.time() * 1000)
5 5
6 parser = argparse.ArgumentParser() 6 parser = argparse.ArgumentParser()
7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file") 7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file")
8 parser.add_argument("--output", help="The annotated summary output file") 8 parser.add_argument("--output", help="The annotated output file to be merged back with the summary file")
9 9
10 args = parser.parse_args() 10 args = parser.parse_args()
11 11
12 infile = args.input 12 infile = args.input
13 #infile = "test_VH-Ca_Cg_25nt/1_Summary_test_VH-Ca_Cg_25nt_241013.txt" 13 #infile = "test_VH-Ca_Cg_25nt/1_Summary_test_VH-Ca_Cg_25nt_241013.txt"
110 if ID not in hits.keys(): #ensure that the dictionairy that keeps track of the hits for every gene exists 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} 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] 112 currentIDHits = hits[ID]
113 seq = dic[ID] 113 seq = dic[ID]
114 lastindex = 0 114 lastindex = 0
115 start = [0] * len(seq) 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)
116 for i, regexp in enumerate(regularexpressions): #for every regular expression 117 for i, regexp in enumerate(regularexpressions): #for every regular expression
117 relativeStartLocation = lastindex - (chunklength / 2) * i 118 relativeStartLocation = lastindex - (chunklength / 2) * i
118 if relativeStartLocation < 0 or relativeStartLocation >= len(seq): 119 if relativeStartLocation >= len(seq):
119 break 120 break
120 regex, hasVar = regexp 121 regex, hasVar = regexp
121 matches = regex.finditer(seq[lastindex:]) 122 matches = regex.finditer(seq[lastindex:])
122 for match in matches: #for every match with the current regex, only uses the first hit 123 for match in matches: #for every match with the current regex, only uses the first hit
123 lastindex += match.start() 124 lastindex += match.start()
124 start[relativeStartLocation] += 1 125 start[relativeStartLocation + start_zero] += 1
125 if hasVar: #if the regex has a variable nt in it 126 if hasVar: #if the regex has a variable nt in it
126 chunkstart = chunklength / 2 * i #where in the reference does this chunk start 127 chunkstart = chunklength / 2 * i #where in the reference does this chunk start
127 chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end 128 chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end
128 if key == "ca": #just calculate the variable nt score for 'ca', cheaper 129 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
129 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]]) 130 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]])
138 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped 139 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped
139 else: #only runs if there were no hits 140 else: #only runs if there were no hits
140 continue 141 continue
141 #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] 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]
142 currentIDHits[key + "_hits"] += 1 143 currentIDHits[key + "_hits"] += 1
143 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1) for x in range(5) if len(start) > 0 and max(start) > 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])
144 #start_location[ID + "_" + key] = str(start.index(max(start))) 145 #start_location[ID + "_" + key] = str(start.index(max(start)))
145 146
146 147
147 chunksInCA = len(compiledregex["ca"]) 148 chunksInCA = len(compiledregex["ca"])
148 chunksInCG = len(compiledregex["cg"]) 149 chunksInCG = len(compiledregex["cg"])
159 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence 160 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
160 with open(output, 'w') as o: 161 with open(output, 'w') as o:
161 for line in f: 162 for line in f:
162 total += 1 163 total += 1
163 if first: 164 if first:
164 o.write(line.rstrip() + "\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n") 165 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
165 first = False 166 first = False
166 continue 167 continue
167 linesplt = line.split("\t") 168 linesplt = line.split("\t")
168 if linesplt[2] == "No results": 169 if linesplt[2] == "No results":
169 pass 170 pass
170 ID = linesplt[1] 171 ID = linesplt[1]
171 currentIDHits = hits[ID] 172 currentIDHits = hits[ID]
173 print currentIDHits
172 possibleca = float(len(compiledregex["ca"])) 174 possibleca = float(len(compiledregex["ca"]))
173 possiblecg = float(len(compiledregex["cg"])) 175 possiblecg = float(len(compiledregex["cg"]))
174 possiblecm = float(len(compiledregex["cm"])) 176 possiblecm = float(len(compiledregex["cm"]))
175 cahits = currentIDHits["ca_hits"] 177 cahits = currentIDHits["ca_hits"]
176 cghits = currentIDHits["cg_hits"] 178 cghits = currentIDHits["cg_hits"]
177 cmhits = currentIDHits["cm_hits"] 179 cmhits = currentIDHits["cm_hits"]
178 if cahits >= cghits and cahits >= cmhits: #its a ca gene 180 if cahits >= cghits and cahits >= cmhits: #its a ca gene
179 ca1hits = currentIDHits["ca1"] 181 ca1hits = currentIDHits["ca1"]
180 ca2hits = currentIDHits["ca2"] 182 ca2hits = currentIDHits["ca2"]
181 if ca1hits >= ca2hits: 183 if ca1hits >= ca2hits:
182 o.write(line.rstrip() + "\tca1\t" + str(int(ca1hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") 184 o.write(ID + "\tca1\t" + str(int(ca1hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
183 else: 185 else:
184 o.write(line.rstrip() + "\tca2\t" + str(int(ca2hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") 186 o.write(ID + "\tca2\t" + str(int(ca2hits / varsInCA * 100)) + "\t" + str(int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
185 elif cghits >= cahits and cghits >= cmhits: #its a cg gene 187 elif cghits >= cahits and cghits >= cmhits: #its a cg gene
186 cg1hits = currentIDHits["cg1"] 188 cg1hits = currentIDHits["cg1"]
187 cg2hits = currentIDHits["cg2"] 189 cg2hits = currentIDHits["cg2"]
188 cg3hits = currentIDHits["cg3"] 190 cg3hits = currentIDHits["cg3"]
189 cg4hits = currentIDHits["cg4"] 191 cg4hits = currentIDHits["cg4"]
190 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene 192 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene
191 o.write(line.rstrip() + "\tcg1\t" + str(int(cg1hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 193 o.write(ID + "\tcg1\t" + str(int(cg1hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
192 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene 194 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene
193 o.write(line.rstrip() + "\tcg2\t" + str(int(cg2hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 195 o.write(ID + "\tcg2\t" + str(int(cg2hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
194 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene 196 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene
195 o.write(line.rstrip() + "\tcg3\t" + str(int(cg3hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 197 o.write(ID + "\tcg3\t" + str(int(cg3hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
196 else: #cg4 gene 198 else: #cg4 gene
197 o.write(line.rstrip() + "\tcg4\t" + str(int(cg4hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 199 o.write(ID + "\tcg4\t" + str(int(cg4hits / varsInCG * 100)) + "\t" + str(int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
198 else: #its a cm gene 200 else: #its a cm gene
199 o.write(line.rstrip() + "\tcm\t0\t" + str(int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 201 o.write(ID + "\tcm\t0\t" + str(int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
200 seq_write_count += 1 202 seq_write_count += 1
201 203
202 print "Time: %i" % (int(time.time() * 1000) - starttime) 204 print "Time: %i" % (int(time.time() * 1000) - starttime)
203 205
204 print "Number of sequences written to file:", seq_write_count 206 print "Number of sequences written to file:", seq_write_count