Mercurial > repos > davidvanzessen > shm_csr
comparison gene_identification.py @ 83:729738462297 draft
"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 |