Mercurial > repos > davidvanzessen > shm_csr
comparison gene_identification.py @ 92:cf8ad181628f draft
planemo upload commit 36be3b053802693392f935e6619ba3f2b1704e3c
author | rhpvorderman |
---|---|
date | Mon, 12 Dec 2022 12:32:44 +0000 |
parents | 729738462297 |
children |
comparison
equal
deleted
inserted
replaced
91:f387cc1580c6 | 92:cf8ad181628f |
---|---|
1 #!/usr/bin/env python3 | |
2 | |
3 import argparse | |
1 import re | 4 import re |
2 import argparse | 5 from typing import Dict, Iterator, List, Tuple |
3 import time | 6 |
4 starttime= int(time.time() * 1000) | 7 |
5 | 8 def generate_sequence_and_id_from_summary(summary_file: str |
6 parser = argparse.ArgumentParser() | 9 ) -> Iterator[Tuple[str, str]]: |
7 parser.add_argument("--input", help="The 1_Summary file from an IMGT zip file") | 10 with open(summary_file, "rt") as summary: |
8 parser.add_argument("--output", help="The annotated output file to be merged back with the summary file") | 11 header = next(summary) |
9 | 12 column_names = header.strip("\n").split("\t") |
10 args = parser.parse_args() | 13 id_column = column_names.index("Sequence ID") |
11 | 14 sequence_column = column_names.index("Sequence") |
12 infile = args.input | 15 for line in summary: |
13 #infile = "test_VH-Ca_Cg_25nt/1_Summary_test_VH-Ca_Cg_25nt_241013.txt" | 16 values = line.strip("\n").split("\t") |
14 output = args.output | 17 id = values[id_column] |
15 #outfile = "identified.txt" | 18 try: |
16 | 19 sequence = values[sequence_column] |
17 dic = dict() | 20 except IndexError: # weird rows without a sequence |
18 total = 0 | 21 sequence = "" |
19 | 22 yield id, sequence |
20 | 23 |
21 first = True | 24 |
22 IDIndex = 0 | 25 # old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc |
23 seqIndex = 0 | 26 # old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctg |
24 | 27 # ggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcagg |
25 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence | 28 # cgccctgaccag |
26 for line in f: | 29 SEARCHSTRINGS = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccag" |
27 total += 1 | 30 "atgggaacgtggtcatcgcctgcctgg", |
28 linesplt = line.split("\t") | 31 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctc" |
29 if first: | 32 "tgggggcacagcggcc", |
30 print("linesplt", linesplt) | 33 "ce": "gcctccacacagagcccatccgtcttccccttgacccgctgctgcaaaaacatt" |
31 IDIndex = linesplt.index("Sequence ID") | 34 "ccctcc", |
32 seqIndex = linesplt.index("Sequence") | |
33 first = False | |
34 continue | |
35 | |
36 ID = linesplt[IDIndex] | |
37 if len(linesplt) < 28: #weird rows without a sequence | |
38 dic[ID] = "" | |
39 else: | |
40 dic[ID] = linesplt[seqIndex] | |
41 | |
42 print("Number of input sequences:", len(dic)) | |
43 | |
44 #old cm sequence: gggagtgcatccgccccaacccttttccccctcgtctcctgtgagaattccc | |
45 #old cg sequence: ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggccctgggctgcctggtcaaggactacttccccgaaccggtgacggtgtcgtggaactcaggcgccctgaccag | |
46 | |
47 #lambda/kappa reference sequence | |
48 searchstrings = {"ca": "catccccgaccagccccaaggtcttcccgctgagcctctgcagcacccagccagatgggaacgtggtcatcgcctgcctgg", | |
49 "cg": "ctccaccaagggcccatcggtcttccccctggcaccctcctccaagagcacctctgggggcacagcggcc", | |
50 "ce": "gcctccacacagagcccatccgtcttccccttgacccgctgctgcaaaaacattccctcc", | |
51 "cm": "gggagtgcatccgccccaacc"} #new (shorter) cm sequence | 35 "cm": "gggagtgcatccgccccaacc"} #new (shorter) cm sequence |
52 | 36 |
53 compiledregex = {"ca": [], | 37 #lambda/kappa referesearchstringsnce sequence variable nucleotides |
54 "cg": [], | 38 CA1_MUTATIONS = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'} |
55 "ce": [], | 39 CA2_MUTATIONS = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'} |
56 "cm": []} | 40 CG1_MUTATIONS = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'} |
57 | 41 CG2_MUTATIONS = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'} |
58 #lambda/kappa reference sequence variable nucleotides | 42 CG3_MUTATIONS = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'} |
59 ca1 = {38: 't', 39: 'g', 48: 'a', 49: 'g', 51: 'c', 68: 'a', 73: 'c'} | 43 CG4_MUTATIONS = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'} |
60 ca2 = {38: 'g', 39: 'a', 48: 'c', 49: 'c', 51: 'a', 68: 'g', 73: 'a'} | |
61 cg1 = {0: 'c', 33: 'a', 38: 'c', 44: 'a', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'} | |
62 cg2 = {0: 'c', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'g', 132: 't'} | |
63 cg3 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 't', 56: 'g', 58: 'g', 66: 'g', 132: 'c'} | |
64 cg4 = {0: 't', 33: 'g', 38: 'g', 44: 'g', 54: 'c', 56: 'a', 58: 'a', 66: 'c', 132: 'c'} | |
65 | 44 |
66 #remove last snp for shorter cg sequence --- note, also change varsInCG | 45 #remove last snp for shorter cg sequence --- note, also change varsInCG |
67 del cg1[132] | 46 del CG1_MUTATIONS[132] |
68 del cg2[132] | 47 del CG2_MUTATIONS[132] |
69 del cg3[132] | 48 del CG3_MUTATIONS[132] |
70 del cg4[132] | 49 del CG4_MUTATIONS[132] |
71 | 50 |
72 #reference sequences are cut into smaller parts of 'chunklength' length, and with 'chunklength' / 2 overlap | 51 # reference sequences are cut into smaller parts of 'chunklength' length, |
73 chunklength = 8 | 52 # and with 'chunklength' / 2 overlap |
74 | 53 CHUNK_LENGTH = 8 |
75 #create the chunks of the reference sequence with regular expressions for the variable nucleotides | 54 |
76 for i in range(0, len(searchstrings["ca"]) - chunklength, chunklength // 2): | 55 |
77 pos = i | 56 def create_compiled_regexes() -> Dict[str, List[Tuple[re.Pattern, int]]]: |
78 chunk = searchstrings["ca"][i:i+chunklength] | 57 |
79 result = "" | 58 compiledregex: Dict[str, List[Tuple[re.Pattern, int]]] = { |
80 varsInResult = 0 | 59 "ca": [], |
81 for c in chunk: | 60 "cg": [], |
82 if pos in list(ca1.keys()): | 61 "ce": [], |
83 varsInResult += 1 | 62 "cm": [] |
84 result += "[" + ca1[pos] + ca2[pos] + "]" | 63 } |
85 else: | 64 |
86 result += c | 65 for i in range(0, len(SEARCHSTRINGS["ca"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2): |
87 pos += 1 | 66 pos = i |
88 compiledregex["ca"].append((re.compile(result), varsInResult)) | 67 chunk = SEARCHSTRINGS["ca"][i:i + CHUNK_LENGTH] |
89 | 68 result = "" |
90 for i in range(0, len(searchstrings["cg"]) - chunklength, chunklength // 2): | 69 varsInResult = 0 |
91 pos = i | 70 for c in chunk: |
92 chunk = searchstrings["cg"][i:i+chunklength] | 71 if pos in list(CA1_MUTATIONS.keys()): |
93 result = "" | 72 varsInResult += 1 |
94 varsInResult = 0 | 73 result += "[" + CA1_MUTATIONS[pos] + CA2_MUTATIONS[pos] + "]" |
95 for c in chunk: | 74 else: |
96 if pos in list(cg1.keys()): | 75 result += c |
97 varsInResult += 1 | 76 pos += 1 |
98 result += "[" + "".join(set([cg1[pos], cg2[pos], cg3[pos], cg4[pos]])) + "]" | 77 compiledregex["ca"].append((re.compile(result), varsInResult)) |
99 else: | 78 |
100 result += c | 79 for i in range(0, len(SEARCHSTRINGS["cg"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2): |
101 pos += 1 | 80 pos = i |
102 compiledregex["cg"].append((re.compile(result), varsInResult)) | 81 chunk = SEARCHSTRINGS["cg"][i:i + CHUNK_LENGTH] |
103 | 82 result = "" |
104 for i in range(0, len(searchstrings["cm"]) - chunklength, chunklength // 2): | 83 varsInResult = 0 |
105 compiledregex["cm"].append((re.compile(searchstrings["cm"][i:i+chunklength]), False)) | 84 for c in chunk: |
106 | 85 if pos in list(CG1_MUTATIONS.keys()): |
107 for i in range(0, len(searchstrings["ce"]) - chunklength + 1, chunklength // 2): | 86 varsInResult += 1 |
108 compiledregex["ce"].append((re.compile(searchstrings["ce"][i:i+chunklength]), False)) | 87 result += "[" + "".join(set([CG1_MUTATIONS[pos], CG2_MUTATIONS[pos], CG3_MUTATIONS[pos], CG4_MUTATIONS[pos]])) + "]" |
88 else: | |
89 result += c | |
90 pos += 1 | |
91 compiledregex["cg"].append((re.compile(result), varsInResult)) | |
92 | |
93 for i in range(0, len(SEARCHSTRINGS["cm"]) - CHUNK_LENGTH, CHUNK_LENGTH // 2): | |
94 compiledregex["cm"].append((re.compile(SEARCHSTRINGS["cm"][i:i + CHUNK_LENGTH]), 0)) | |
95 | |
96 for i in range(0, len(SEARCHSTRINGS["ce"]) - CHUNK_LENGTH + 1, CHUNK_LENGTH // 2): | |
97 compiledregex["ce"].append((re.compile(SEARCHSTRINGS["ce"][i:i + CHUNK_LENGTH]), 0)) | |
98 | |
99 return compiledregex | |
100 | |
109 | 101 |
110 def removeAndReturnMaxIndex(x): #simplifies a list comprehension | 102 def removeAndReturnMaxIndex(x): #simplifies a list comprehension |
111 m = max(x) | 103 m = max(x) |
112 index = x.index(m) | 104 index = x.index(m) |
113 x[index] = 0 | 105 x[index] = 0 |
114 return index | 106 return index |
115 | 107 |
116 | 108 |
117 start_location = dict() | 109 def match_sequence(seq, compiledregex): |
118 hits = dict() | 110 currentIDHits = {"ca_hits": 0, "cg_hits": 0, "cm_hits": 0, "ce_hits": 0, |
119 alltotal = 0 | 111 "ca1": 0, "ca2": 0, "cg1": 0, "cg2": 0, "cg3": 0, "cg4": 0} |
120 for key in compiledregex: #for ca/cg/cm/ce | 112 alltotal = 0 |
121 regularexpressions = compiledregex[key] # get the compiled regular expressions | 113 start_location = dict() |
122 for ID in list(dic.keys())[0:]: #for every ID | 114 for key in compiledregex: # for ca/cg/cm/ce |
123 if ID not in list(hits.keys()): #ensure that the dictionairy that keeps track of the hits for every gene exists | 115 regularexpressions = compiledregex[key] |
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] | |
126 seq = dic[ID] | |
127 lastindex = 0 | 116 lastindex = 0 |
128 start_zero = len(searchstrings[key]) #allows the reference sequence to start before search sequence (start_locations of < 0) | 117 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) | 118 start = [0] * (len(seq) + start_zero) |
130 for i, regexp in enumerate(regularexpressions): #for every regular expression | 119 for i, regexp in enumerate(regularexpressions): #for every regular expression |
131 relativeStartLocation = lastindex - (chunklength // 2) * i | 120 relativeStartLocation = lastindex - (CHUNK_LENGTH // 2) * i |
132 if relativeStartLocation >= len(seq): | 121 if relativeStartLocation >= len(seq): |
133 break | 122 break |
134 regex, hasVar = regexp | 123 regex, hasVar = regexp |
135 matches = regex.finditer(seq[lastindex:]) | 124 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 | 125 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() | 126 lastindex += match.start() |
138 start[relativeStartLocation + start_zero] += 1 | 127 start[relativeStartLocation + start_zero] += 1 |
139 if hasVar: #if the regex has a variable nt in it | 128 if hasVar: #if the regex has a variable nt in it |
140 chunkstart = chunklength // 2 * i #where in the reference does this chunk start | 129 chunkstart = CHUNK_LENGTH // 2 * i #where in the reference does this chunk start |
141 chunkend = chunklength // 2 * i + chunklength #where in the reference does this chunk end | 130 chunkend = CHUNK_LENGTH // 2 * i + CHUNK_LENGTH #where in the reference does this chunk end |
142 if key == "ca": #just calculate the variable nt score for 'ca', cheaper | 131 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]]) | 132 currentIDHits["ca1"] += len([1 for x in CA1_MUTATIONS if chunkstart <= x < chunkend and CA1_MUTATIONS[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]]) | 133 currentIDHits["ca2"] += len([1 for x in CA2_MUTATIONS if chunkstart <= x < chunkend and CA2_MUTATIONS[x] == seq[lastindex + x - chunkstart]]) |
145 elif key == "cg": #just calculate the variable nt score for 'cg', cheaper | 134 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]]) | 135 currentIDHits["cg1"] += len([1 for x in CG1_MUTATIONS if chunkstart <= x < chunkend and CG1_MUTATIONS[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]]) | 136 currentIDHits["cg2"] += len([1 for x in CG2_MUTATIONS if chunkstart <= x < chunkend and CG2_MUTATIONS[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]]) | 137 currentIDHits["cg3"] += len([1 for x in CG3_MUTATIONS if chunkstart <= x < chunkend and CG3_MUTATIONS[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]]) | 138 currentIDHits["cg4"] += len([1 for x in CG4_MUTATIONS if chunkstart <= x < chunkend and CG4_MUTATIONS[x] == seq[lastindex + x - chunkstart]]) |
150 else: #key == "cm" #no variable regions in 'cm' or 'ce' | 139 else: #key == "cm" #no variable regions in 'cm' or 'ce' |
151 pass | 140 pass |
152 break #this only breaks when there was a match with the regex, breaking means the 'else:' clause is skipped | 141 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 | 142 else: #only runs if there were no hits |
154 continue | 143 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] | 144 #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 | 145 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]) | 146 start_location[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))) | 147 |
159 | 148 cahits = currentIDHits["ca_hits"] |
160 | 149 cghits = currentIDHits["cg_hits"] |
161 varsInCA = float(len(list(ca1.keys())) * 2) | 150 cmhits = currentIDHits["cm_hits"] |
162 varsInCG = float(len(list(cg1.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice | 151 cehits = currentIDHits["ce_hits"] |
163 varsInCM = 0 | 152 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: # its a ca gene |
164 varsInCE = 0 | 153 ca1hits = currentIDHits["ca1"] |
165 | 154 ca2hits = currentIDHits["ca2"] |
166 def round_int(val): | 155 if ca1hits >= ca2hits: |
167 return int(round(val)) | 156 # TODO: All variants with 0 matched are matched to IGA1 with 0 hits |
168 | 157 # TODO: these are later turned into unmatched by the merge_and_filter.R |
169 first = True | 158 # TODO: script |
170 seq_write_count=0 | 159 return "IGA1", ca1hits, cahits, start_location["ca"] |
171 with open(infile, 'r') as f: #read all sequences into a dictionary as key = ID, value = sequence | 160 else: |
172 with open(output, 'w') as o: | 161 return "IGA2", ca2hits, cahits, start_location["ca"] |
173 for line in f: | 162 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: # its a cg gene |
174 total += 1 | 163 cg1hits = currentIDHits["cg1"] |
175 if first: | 164 cg2hits = currentIDHits["cg2"] |
176 o.write("Sequence ID\tbest_match\tnt_hit_percentage\tchunk_hit_percentage\tstart_locations\n") | 165 cg3hits = currentIDHits["cg3"] |
177 first = False | 166 cg4hits = currentIDHits["cg4"] |
178 continue | 167 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: # cg1 gene |
179 linesplt = line.split("\t") | 168 return "IGG1", cg1hits, cghits, start_location["cg"] |
180 if linesplt[2] == "No results": | 169 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: # cg2 gene |
181 pass | 170 return "IGG2", cg2hits, cghits, start_location["cg"] |
182 ID = linesplt[1] | 171 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: # cg3 gene |
183 currentIDHits = hits[ID] | 172 return "IGG3", cg3hits, cghits, start_location["cg"] |
184 possibleca = float(len(compiledregex["ca"])) | 173 else: # cg4 gene |
185 possiblecg = float(len(compiledregex["cg"])) | 174 return "IGG4", cg4hits, cghits, start_location["cg"] |
186 possiblecm = float(len(compiledregex["cm"])) | 175 else: # its a cm or ce gene |
187 possiblece = float(len(compiledregex["ce"])) | 176 if cmhits >= cehits: |
188 cahits = currentIDHits["ca_hits"] | 177 return "IGM", 0, cmhits, start_location["cm"] |
189 cghits = currentIDHits["cg_hits"] | 178 else: |
190 cmhits = currentIDHits["cm_hits"] | 179 return "IGE", 0, cehits, start_location["ce"] |
191 cehits = currentIDHits["ce_hits"] | 180 |
192 if cahits >= cghits and cahits >= cmhits and cahits >= cehits: #its a ca gene | 181 |
193 ca1hits = currentIDHits["ca1"] | 182 def main(): |
194 ca2hits = currentIDHits["ca2"] | 183 parser = argparse.ArgumentParser() |
195 if ca1hits >= ca2hits: | 184 parser.add_argument("--input", |
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") | 185 help="The 1_Summary file from an IMGT zip file") |
197 else: | 186 parser.add_argument("--output", |
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") | 187 help="The annotated output file to be merged back " |
199 elif cghits >= cahits and cghits >= cmhits and cghits >= cehits: #its a cg gene | 188 "with the summary file") |
200 cg1hits = currentIDHits["cg1"] | 189 args = parser.parse_args() |
201 cg2hits = currentIDHits["cg2"] | 190 varsInCA = float(len(list(CA1_MUTATIONS.keys())) * 2) |
202 cg3hits = currentIDHits["cg3"] | 191 varsInCG = float(len(list( |
203 cg4hits = currentIDHits["cg4"] | 192 CG1_MUTATIONS.keys())) * 2) - 2 # -2 because the sliding window doesn't hit the first and last nt twice |
204 if cg1hits >= cg2hits and cg1hits >= cg3hits and cg1hits >= cg4hits: #cg1 gene | 193 subclass_vars = { |
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") | 194 "IGA1": varsInCA, "IGA2": varsInCA, |
206 elif cg2hits >= cg1hits and cg2hits >= cg3hits and cg2hits >= cg4hits: #cg2 gene | 195 "IGG1": varsInCG, "IGG2": varsInCG, "IGG3": varsInCG, "IGG4": varsInCG, |
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") | 196 "IGE": 0, |
208 elif cg3hits >= cg1hits and cg3hits >= cg2hits and cg3hits >= cg4hits: #cg3 gene | 197 "IGM": 0, |
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") | 198 } |
210 else: #cg4 gene | 199 compiledregex = create_compiled_regexes() |
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") | 200 possibleca = float(len(compiledregex["ca"])) |
212 else: #its a cm or ce gene | 201 possiblecg = float(len(compiledregex["cg"])) |
213 if cmhits >= cehits: | 202 possiblecm = float(len(compiledregex["cm"])) |
214 o.write(ID + "\tIGM\t100\t" + str(round_int(cmhits / possiblecm * 100)) + "\t" + start_location[ID + "_cm"] + "\n") | 203 possiblece = float(len(compiledregex["ce"])) |
215 else: | 204 class_chunks = { |
216 o.write(ID + "\tIGE\t100\t" + str(round_int(cehits / possiblece * 100)) + "\t" + start_location[ID + "_ce"] + "\n") | 205 "IGA1": possibleca, "IGA2": possibleca, |
217 seq_write_count += 1 | 206 "IGE": possiblece, |
218 | 207 "IGG1": possiblecg, "IGG2": possiblecg, "IGG3": possiblecg, |
219 print("Time: %i" % (int(time.time() * 1000) - starttime)) | 208 "IGG4": possiblecg, |
220 | 209 "IGM": possiblecm |
221 print("Number of sequences written to file:", seq_write_count) | 210 } |
222 | 211 with open(args.output, "wt") as output: |
223 | 212 output.write("Sequence ID\tbest_match\tnt_hit_percentage\t" |
224 | 213 "chunk_hit_percentage\tstart_locations\n") |
225 | 214 for id, sequence in generate_sequence_and_id_from_summary(args.input): |
226 | 215 best_match, subclass_hits, class_hits, start_locations = \ |
216 match_sequence(sequence, compiledregex) | |
217 variable_nucs = subclass_vars[best_match] | |
218 if variable_nucs: | |
219 subclass_percentage = round(subclass_hits * 100 / | |
220 variable_nucs) | |
221 else: | |
222 subclass_percentage = 100 | |
223 class_percentage = round(class_hits * 100 / class_chunks[best_match]) | |
224 output.write(f"{id}\t{best_match}\t{subclass_percentage}\t" | |
225 f"{class_percentage}\t{start_locations}\n") | |
226 | |
227 | |
228 if __name__ == "__main__": | |
229 main() |