Mercurial > repos > davidvanzessen > shm_csr

comparison gene_identification.py @ 83:729738462297 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
author rhpvorderman
date Wed, 15 Sep 2021 12:24:06 +0000
parents b6f9a640e098
children cf8ad181628f
comparison
equal deleted inserted replaced
82:a103134ee6e0 83:729738462297
21 first = True 21 first = True
22 IDIndex = 0 22 IDIndex = 0
23 seqIndex = 0 23 seqIndex = 0
24 24
25 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence 25 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
26 for line in f: 26 for line in f:
27 total += 1 27 total += 1
28 linesplt = line.split("\t") 28 linesplt = line.split("\t")
29 if first: 29 if first:
30 print "linesplt", linesplt 30 print("linesplt", linesplt)
31 IDIndex = linesplt.index("Sequence ID") 31 IDIndex = linesplt.index("Sequence ID")
32 seqIndex = linesplt.index("Sequence") 32 seqIndex = linesplt.index("Sequence")
33 first = False 33 first = False
34 continue 34 continue
35 35
36 ID = linesplt[IDIndex] 36 ID = linesplt[IDIndex]
37 if len(linesplt) < 28: #weird rows without a sequence 37 if len(linesplt) < 28: #weird rows without a sequence
38 dic[ID] = "" 38 dic[ID] = ""
39 else: 39 else:
40 dic[ID] = linesplt[seqIndex] 40 dic[ID] = linesplt[seqIndex]
41 41
42 print "Number of input sequences:", len(dic) 42 print("Number of input sequences:", len(dic))
43 43
44 #old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc 44 #old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc
45 #old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag 45 #old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag
46 46
47 #lambda/kappa reference sequence 47 #lambda/kappa reference sequence
71 71
72 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap 72 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap
73 chunklength = 8 73 chunklength = 8
74 74
75 #create the chunks of the reference sequence with regular expressions for the variable nucleotides 75 #create the chunks of the reference sequence with regular expressions for the variable nucleotides
76 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength / 2): 76 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength // 2):
77 pos = i 77 pos = i
78 chunk = searchstrings["ca"][i:i+chunklength] 78 chunk = searchstrings["ca"][i:i+chunklength]
79 result = "" 79 result = ""
80 varsInResult = 0 80 varsInResult = 0
81 for c in chunk: 81 for c in chunk:
82 if pos in ca1.keys(): 82 if pos in list(ca1.keys()):
83 varsInResult += 1 83 varsInResult += 1
84 result += "[" + ca1[pos] + ca2[pos] + "]" 84 result += "[" + ca1[pos] + ca2[pos] + "]"
85 else: 85 else:
86 result += c 86 result += c
87 pos += 1 87 pos += 1
88 compiledregex["ca"].append((re.compile(result), varsInResult)) 88 compiledregex["ca"].append((re.compile(result), varsInResult))
89 89
90 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength / 2): 90 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength // 2):
91 pos = i 91 pos = i
92 chunk = searchstrings["cg"][i:i+chunklength] 92 chunk = searchstrings["cg"][i:i+chunklength]
93 result = "" 93 result = ""
94 varsInResult = 0 94 varsInResult = 0
95 for c in chunk: 95 for c in chunk:
96 if pos in cg1.keys(): 96 if pos in list(cg1.keys()):
97 varsInResult += 1 97 varsInResult += 1
98 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]" 98 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]"
99 else: 99 else:
100 result += c 100 result += c
101 pos += 1 101 pos += 1
102 compiledregex["cg"].append((re.compile(result), varsInResult)) 102 compiledregex["cg"].append((re.compile(result), varsInResult))
103 103
104 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength / 2): 104 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength // 2):
105 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False)) 105 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False))
106 106
107 for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength / 2): 107 for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength // 2):
108 compiledregex["ce"].append((re.compile(searchstrings["ce"][i:i+chunklength]), False)) 108 compiledregex["ce"].append((re.compile(searchstrings["ce"][i:i+chunklength]), False))
109 109
110 def removeAndReturnMaxIndex(x): #simplifies a list comprehension 110 def removeAndReturnMaxIndex(x): #simplifies a list comprehension
111 m = max(x) 111 m = max(x)
112 index = x.index(m) 112 index = x.index(m)
115 115
116 116
117 start_location = dict() 117 start_location = dict()
118 hits = dict() 118 hits = dict()
119 alltotal = 0 119 alltotal = 0
120 for key in compiledregex.keys(): #for ca/cg/cm/ce 120 for key in compiledregex: #for ca/cg/cm/ce
121 regularexpressions = compiledregex[key] #get the compiled regular expressions 121 regularexpressions = compiledregex[key] # get the compiled regular expressions
122 for ID in dic.keys()[0:]: #for every ID 122 for ID in list(dic.keys())[0:]: #for every ID
123 if ID not in hits.keys(): #ensure that the dictionairy that keeps track of the hits for every gene exists 123 if ID not in list(hits.keys()): #ensure that the dictionairy that keeps track of the hits for every gene exists
124 hits[ID] = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ce_hits": 0, "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0} 124 hits[ID] = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ce_hits": 0, "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0}
125 currentIDHits = hits[ID] 125 currentIDHits = hits[ID]
126 seq = dic[ID] 126 seq = dic[ID]
127 lastindex = 0 127 lastindex = 0
128 start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0) 128 start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0)
129 start = [0] * (len(seq) + start_zero) 129 start = [0] * (len(seq) + start_zero)
130 for i, regexp in enumerate(regularexpressions): #for every regular expression 130 for i, regexp in enumerate(regularexpressions): #for every regular expression
131 relativeStartLocation = lastindex - (chunklength / 2) * i 131 relativeStartLocation = lastindex - (chunklength // 2) * i
132 if relativeStartLocation >= len(seq): 132 if relativeStartLocation >= len(seq):
133 break 133 break
134 regex, hasVar = regexp 134 regex, hasVar = regexp
135 matches = regex.finditer(seq[lastindex:]) 135 matches = regex.finditer(seq[lastindex:])
136 for match in matches: #for every match with the current regex, only uses the first hit because of the break at the end of this loop 136 for match in matches: #for every match with the current regex, only uses the first hit because of the break at the end of this loop
137 lastindex += match.start() 137 lastindex += match.start()
138 start[relativeStartLocation + start_zero] += 1 138 start[relativeStartLocation + start_zero] += 1
139 if hasVar: #if the regex has a variable nt in it 139 if hasVar: #if the regex has a variable nt in it
140 chunkstart = chunklength / 2 * i #where in the reference does this chunk start 140 chunkstart = chunklength // 2 * i #where in the reference does this chunk start
141 chunkend = chunklength / 2 * i + chunklength #where in the reference does this chunk end 141 chunkend = chunklength // 2 * i + chunklength #where in the reference does this chunk end
142 if key == "ca": #just calculate the variable nt score for 'ca', cheaper 142 if key == "ca": #just calculate the variable nt score for 'ca', cheaper
143 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]]) 143 currentIDHits["ca1"] += len([1 for x in ca1 if chunkstart <= x < chunkend and ca1[x] == seq[lastindex + x - chunkstart]])
144 currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]]) 144 currentIDHits["ca2"] += len([1 for x in ca2 if chunkstart <= x < chunkend and ca2[x] == seq[lastindex + x - chunkstart]])
145 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper 145 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper
146 currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]]) 146 currentIDHits["cg1"] += len([1 for x in cg1 if chunkstart <= x < chunkend and cg1[x] == seq[lastindex + x - chunkstart]])
147 currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]]) 147 currentIDHits["cg2"] += len([1 for x in cg2 if chunkstart <= x < chunkend and cg2[x] == seq[lastindex + x - chunkstart]])
148 currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]]) 148 currentIDHits["cg3"] += len([1 for x in cg3 if chunkstart <= x < chunkend and cg3[x] == seq[lastindex + x - chunkstart]])
149 currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]]) 149 currentIDHits["cg4"] += len([1 for x in cg4 if chunkstart <= x < chunkend and cg4[x] == seq[lastindex + x - chunkstart]])
150 else: #key == "cm" #no variable regions in 'cm' or 'ce' 150 else: #key == "cm" #no variable regions in 'cm' or 'ce'
151 pass 151 pass
152 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped 152 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped
153 else: #only runs if there were no hits 153 else: #only runs if there were no hits
154 continue 154 continue
155 #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] 155 #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]
156 currentIDHits[key + "_hits"] += 1 156 currentIDHits[key + "_hits"] += 1
157 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1]) 157 start_location[ID + "_" + key] = str([(removeAndReturnMaxIndex(start) + 1 - start_zero) for x in range(5) if len(start) > 0 and max(start) > 1])
158 #start_location[ID + "_" + key] = str(start.index(max(start))) 158 #start_location[ID + "_" + key] = str(start.index(max(start)))
159 159
160 160
161 varsInCA = float(len(ca1.keys()) * 2) 161 varsInCA = float(len(list(ca1.keys())) * 2)
162 varsInCG = float(len(cg1.keys()) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice 162 varsInCG = float(len(list(cg1.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice
163 varsInCM = 0 163 varsInCM = 0
164 varsInCE = 0 164 varsInCE = 0
165 165
166 def round_int(val): 166 def round_int(val):
167 return int(round(val)) 167 return int(round(val))
168 168
169 first = True 169 first = True
170 seq_write_count=0 170 seq_write_count=0
171 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence 171 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence
172 with open(output, 'w') as o: 172 with open(output, 'w') as o:
173 for line in f: 173 for line in f:
174 total += 1 174 total += 1
175 if first: 175 if first:
176 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n") 176 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n")
177 first = False 177 first = False
178 continue 178 continue
179 linesplt = line.split("\t") 179 linesplt = line.split("\t")
180 if linesplt[2] == "No results": 180 if linesplt[2] == "No results":
181 pass 181 pass
182 ID = linesplt[1] 182 ID = linesplt[1]
183 currentIDHits = hits[ID] 183 currentIDHits = hits[ID]
184 possibleca = float(len(compiledregex["ca"])) 184 possibleca = float(len(compiledregex["ca"]))
185 possiblecg = float(len(compiledregex["cg"])) 185 possiblecg = float(len(compiledregex["cg"]))
186 possiblecm = float(len(compiledregex["cm"])) 186 possiblecm = float(len(compiledregex["cm"]))
187 possiblece = float(len(compiledregex["ce"])) 187 possiblece = float(len(compiledregex["ce"]))
188 cahits = currentIDHits["ca_hits"] 188 cahits = currentIDHits["ca_hits"]
189 cghits = currentIDHits["cg_hits"] 189 cghits = currentIDHits["cg_hits"]
190 cmhits = currentIDHits["cm_hits"] 190 cmhits = currentIDHits["cm_hits"]
191 cehits = currentIDHits["ce_hits"] 191 cehits = currentIDHits["ce_hits"]
192 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene 192 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene
193 ca1hits = currentIDHits["ca1"] 193 ca1hits = currentIDHits["ca1"]
194 ca2hits = currentIDHits["ca2"] 194 ca2hits = currentIDHits["ca2"]
195 if ca1hits >= ca2hits: 195 if ca1hits >= ca2hits:
196 o.write(ID + "\tIGA1\t" + str(round_int(ca1hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") 196 o.write(ID + "\tIGA1\t" + str(round_int(ca1hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
197 else: 197 else:
198 o.write(ID + "\tIGA2\t" + str(round_int(ca2hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n") 198 o.write(ID + "\tIGA2\t" + str(round_int(ca2hits / varsInCA * 100)) + "\t" + str(round_int(cahits / possibleca * 100)) + "\t" + start_location[ID + "_ca"] + "\n")
199 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene 199 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene
200 cg1hits = currentIDHits["cg1"] 200 cg1hits = currentIDHits["cg1"]
201 cg2hits = currentIDHits["cg2"] 201 cg2hits = currentIDHits["cg2"]
202 cg3hits = currentIDHits["cg3"] 202 cg3hits = currentIDHits["cg3"]
203 cg4hits = currentIDHits["cg4"] 203 cg4hits = currentIDHits["cg4"]
204 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene 204 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene
205 o.write(ID + "\tIGG1\t" + str(round_int(cg1hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 205 o.write(ID + "\tIGG1\t" + str(round_int(cg1hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
206 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene 206 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene
207 o.write(ID + "\tIGG2\t" + str(round_int(cg2hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 207 o.write(ID + "\tIGG2\t" + str(round_int(cg2hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
208 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene 208 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene
209 o.write(ID + "\tIGG3\t" + str(round_int(cg3hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 209 o.write(ID + "\tIGG3\t" + str(round_int(cg3hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
210 else: #cg4 gene 210 else: #cg4 gene
211 o.write(ID + "\tIGG4\t" + str(round_int(cg4hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n") 211 o.write(ID + "\tIGG4\t" + str(round_int(cg4hits / varsInCG * 100)) + "\t" + str(round_int(cghits / possiblecg * 100)) + "\t" + start_location[ID + "_cg"] + "\n")
212 else: #its a cm or ce gene 212 else: #its a cm or ce gene
213 if cmhits >= cehits: 213 if cmhits >= cehits:
214 o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n") 214 o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n")
215 else: 215 else:
216 o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n") 216 o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n")
217 seq_write_count += 1 217 seq_write_count += 1
218 218
219 print "Time: %i" % (int(time.time() * 1000) - starttime) 219 print("Time: %i" % (int(time.time() * 1000) - starttime))
220 220
221 print "Number of sequences written to file:", seq_write_count 221 print("Number of sequences written to file:", seq_write_count)
222 222
223 223
224 224
225 225
226 226