Mercurial > repos > davidvanzessen > shm_csr
annotate shm_csr.py @ 85:a539d14d66b6 draft
"planemo upload commit 95767c430ea956636f95ec28bb749a84eabe0156"
| author | rhpvorderman |
|---|---|
| date | Wed, 27 Oct 2021 09:41:34 +0000 |
| parents | 729738462297 |
| children | 6809c63d9161 |
| rev | line source |
|---|---|
| 81 | 1 import argparse |
| 2 import logging | |
| 3 import sys | |
| 4 import os | |
| 5 import re | |
| 6 | |
| 7 from collections import defaultdict | |
| 8 | |
| 9 def main(): | |
| 10 parser = argparse.ArgumentParser() | |
| 11 parser.add_argument("--input", help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation") | |
| 12 parser.add_argument("--genes", help="The genes available in the 'best_match' column") | |
| 13 parser.add_argument("--empty_region_filter", help="Where does the sequence start?", choices=['leader', 'FR1', 'CDR1', 'FR2']) | |
| 14 parser.add_argument("--output", help="Output file") | |
| 15 | |
| 16 args = parser.parse_args() | |
| 17 | |
| 18 infile = args.input | |
| 19 genes = str(args.genes).split(",") | |
| 20 empty_region_filter = args.empty_region_filter | |
| 21 outfile = args.output | |
| 22 | |
| 23 genedic = dict() | |
| 24 | |
| 25 mutationdic = dict() | |
| 26 mutationMatcher = re.compile("^(.)(\d+).(.),?[ ]?(.)?(\d+)?.?(.)?(.?.?.?.?.?)?") | |
| 27 mutationMatcher = re.compile("^([actg])(\d+).([actg]),?[ ]?([A-Z])?(\d+)?.?([A-Z])?(.*)?") | |
| 28 mutationMatcher = re.compile("^([actg])(\d+).([actg]),?[ ]?([A-Z])?(\d+)?[>]?([A-Z;])?(.*)?") | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
29 mutationMatcher = re.compile(r"^([nactg])(\d+).([nactg]),?[ ]?([A-Z*])?(\d+)?[>]?([A-Z*;])?(.*)?") |
| 81 | 30 NAMatchResult = (None, None, None, None, None, None, '') |
| 31 geneMatchers = {gene: re.compile("^" + gene + ".*") for gene in genes} | |
| 32 linecount = 0 | |
| 33 | |
| 34 IDIndex = 0 | |
| 35 best_matchIndex = 0 | |
| 36 fr1Index = 0 | |
| 37 cdr1Index = 0 | |
| 38 fr2Index = 0 | |
| 39 cdr2Index = 0 | |
| 40 fr3Index = 0 | |
| 41 first = True | |
| 42 IDlist = [] | |
| 43 mutationList = [] | |
| 44 mutationListByID = {} | |
| 45 cdr1LengthDic = {} | |
| 46 cdr2LengthDic = {} | |
| 47 | |
| 48 fr1LengthDict = {} | |
| 49 fr2LengthDict = {} | |
| 50 fr3LengthDict = {} | |
| 51 | |
| 52 cdr1LengthIndex = 0 | |
| 53 cdr2LengthIndex = 0 | |
| 54 | |
| 55 fr1SeqIndex = 0 | |
| 56 fr2SeqIndex = 0 | |
| 57 fr3SeqIndex = 0 | |
| 58 | |
| 59 tandem_sum_by_class = defaultdict(int) | |
| 60 expected_tandem_sum_by_class = defaultdict(float) | |
| 61 | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
62 with open(infile, 'r') as i: |
| 81 | 63 for line in i: |
| 64 if first: | |
| 65 linesplt = line.split("\t") | |
| 66 IDIndex = linesplt.index("Sequence.ID") | |
| 67 best_matchIndex = linesplt.index("best_match") | |
| 68 fr1Index = linesplt.index("FR1.IMGT") | |
| 69 cdr1Index = linesplt.index("CDR1.IMGT") | |
| 70 fr2Index = linesplt.index("FR2.IMGT") | |
| 71 cdr2Index = linesplt.index("CDR2.IMGT") | |
| 72 fr3Index = linesplt.index("FR3.IMGT") | |
| 73 cdr1LengthIndex = linesplt.index("CDR1.IMGT.length") | |
| 74 cdr2LengthIndex = linesplt.index("CDR2.IMGT.length") | |
| 75 fr1SeqIndex = linesplt.index("FR1.IMGT.seq") | |
| 76 fr2SeqIndex = linesplt.index("FR2.IMGT.seq") | |
| 77 fr3SeqIndex = linesplt.index("FR3.IMGT.seq") | |
| 78 first = False | |
| 79 continue | |
| 80 linecount += 1 | |
| 81 linesplt = line.split("\t") | |
| 82 ID = linesplt[IDIndex] | |
| 83 genedic[ID] = linesplt[best_matchIndex] | |
| 84 | |
| 85 mutationdic[ID + "_FR1"] = [] | |
| 86 if len(linesplt[fr1Index]) > 5 and empty_region_filter == "leader": | |
| 87 mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] | |
| 88 | |
| 89 mutationdic[ID + "_CDR1"] = [] | |
| 90 if len(linesplt[cdr1Index]) > 5 and empty_region_filter in ["leader", "FR1"]: | |
| 91 mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x] | |
| 92 | |
| 93 mutationdic[ID + "_FR2"] = [] | |
| 94 if len(linesplt[fr2Index]) > 5 and empty_region_filter in ["leader", "FR1", "CDR1"]: | |
| 95 mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x] | |
| 96 | |
| 97 mutationdic[ID + "_CDR2"] = [] | |
| 98 if len(linesplt[cdr2Index]) > 5: | |
| 99 mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x] | |
| 100 | |
| 101 mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] | |
| 102 | |
| 103 mutationdic[ID + "_FR3"] = [] | |
| 104 if len(linesplt[fr3Index]) > 5: | |
| 105 mutationdic[ID + "_FR3"] = [mutationMatcher.match(x).groups() for x in linesplt[fr3Index].split("|") if x] | |
| 106 | |
| 107 mutationList += mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"] | |
| 108 mutationListByID[ID] = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"] | |
| 109 | |
| 110 try: | |
| 111 cdr1Length = int(linesplt[cdr1LengthIndex]) | |
| 112 except: | |
| 113 cdr1Length = 0 | |
| 114 | |
| 115 try: | |
| 116 cdr2Length = int(linesplt[cdr2LengthIndex]) | |
| 117 except: | |
| 118 cdr2Length = 0 | |
| 119 | |
| 120 #print linesplt[fr2SeqIndex] | |
| 121 fr1Length = len(linesplt[fr1SeqIndex]) if empty_region_filter == "leader" else 0 | |
| 122 fr2Length = len(linesplt[fr2SeqIndex]) if empty_region_filter in ["leader", "FR1", "CDR1"] else 0 | |
| 123 fr3Length = len(linesplt[fr3SeqIndex]) | |
| 124 | |
| 125 cdr1LengthDic[ID] = cdr1Length | |
| 126 cdr2LengthDic[ID] = cdr2Length | |
| 127 | |
| 128 fr1LengthDict[ID] = fr1Length | |
| 129 fr2LengthDict[ID] = fr2Length | |
| 130 fr3LengthDict[ID] = fr3Length | |
| 131 | |
| 132 IDlist += [ID] | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
133 print("len(mutationdic) =", len(mutationdic)) |
| 81 | 134 |
| 135 with open(os.path.join(os.path.dirname(os.path.abspath(infile)), "mutationdict.txt"), 'w') as out_handle: | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
136 for ID, lst in mutationdic.items(): |
| 81 | 137 for mut in lst: |
| 138 out_handle.write("{0}\t{1}\n".format(ID, "\t".join([str(x) for x in mut]))) | |
| 139 | |
| 140 #tandem mutation stuff | |
| 141 tandem_frequency = defaultdict(int) | |
| 142 mutation_frequency = defaultdict(int) | |
| 143 | |
| 144 mutations_by_id_dic = {} | |
| 145 first = True | |
| 146 mutation_by_id_file = os.path.join(os.path.dirname(outfile), "mutation_by_id.txt") | |
| 147 with open(mutation_by_id_file, 'r') as mutation_by_id: | |
| 148 for l in mutation_by_id: | |
| 149 if first: | |
| 150 first = False | |
| 151 continue | |
| 152 splt = l.split("\t") | |
| 153 mutations_by_id_dic[splt[0]] = int(splt[1]) | |
| 154 | |
| 155 tandem_file = os.path.join(os.path.dirname(outfile), "tandems_by_id.txt") | |
| 156 with open(tandem_file, 'w') as o: | |
| 157 highest_tandem_length = 0 | |
| 158 | |
| 159 o.write("Sequence.ID\tnumber_of_mutations\tnumber_of_tandems\tregion_length\texpected_tandems\tlongest_tandem\ttandems\n") | |
| 160 for ID in IDlist: | |
| 161 mutations = mutationListByID[ID] | |
| 162 if len(mutations) == 0: | |
| 163 continue | |
| 164 last_mut = max(mutations, key=lambda x: int(x[1])) | |
| 165 | |
| 166 last_mut_pos = int(last_mut[1]) | |
| 167 | |
| 168 mut_positions = [False] * (last_mut_pos + 1) | |
| 169 | |
| 170 for mutation in mutations: | |
| 171 frm, where, to, frmAA, whereAA, toAA, thing = mutation | |
| 172 where = int(where) | |
| 173 mut_positions[where] = True | |
| 174 | |
| 175 tandem_muts = [] | |
| 176 tandem_start = -1 | |
| 177 tandem_length = 0 | |
| 178 for i in range(len(mut_positions)): | |
| 179 if mut_positions[i]: | |
| 180 if tandem_start == -1: | |
| 181 tandem_start = i | |
| 182 tandem_length += 1 | |
| 183 #print "".join(["1" if x else "0" for x in mut_positions[:i+1]]) | |
| 184 else: | |
| 185 if tandem_length > 1: | |
| 186 tandem_muts.append((tandem_start, tandem_length)) | |
| 187 #print "{0}{1} {2}:{3}".format(" " * (i - tandem_length), "^" * tandem_length, tandem_start, tandem_length) | |
| 188 tandem_start = -1 | |
| 189 tandem_length = 0 | |
| 190 if tandem_length > 1: # if the sequence ends with a tandem mutation | |
| 191 tandem_muts.append((tandem_start, tandem_length)) | |
| 192 | |
| 193 if len(tandem_muts) > 0: | |
| 194 if highest_tandem_length < len(tandem_muts): | |
| 195 highest_tandem_length = len(tandem_muts) | |
| 196 | |
| 197 region_length = fr1LengthDict[ID] + cdr1LengthDic[ID] + fr2LengthDict[ID] + cdr2LengthDic[ID] + fr3LengthDict[ID] | |
| 198 longest_tandem = max(tandem_muts, key=lambda x: x[1]) if len(tandem_muts) else (0, 0) | |
| 199 num_mutations = mutations_by_id_dic[ID] # len(mutations) | |
| 200 f_num_mutations = float(num_mutations) | |
| 201 num_tandem_muts = len(tandem_muts) | |
| 202 expected_tandem_muts = f_num_mutations * (f_num_mutations - 1.0) / float(region_length) | |
| 203 o.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\n".format(ID, | |
| 204 str(num_mutations), | |
| 205 str(num_tandem_muts), | |
| 206 str(region_length), | |
| 207 str(round(expected_tandem_muts, 2)), | |
| 208 str(longest_tandem[1]), | |
| 209 str(tandem_muts))) | |
| 210 gene = genedic[ID] | |
| 211 if gene.find("unmatched") == -1: | |
| 212 tandem_sum_by_class[gene] += num_tandem_muts | |
| 213 expected_tandem_sum_by_class[gene] += expected_tandem_muts | |
| 214 | |
| 215 tandem_sum_by_class["all"] += num_tandem_muts | |
| 216 expected_tandem_sum_by_class["all"] += expected_tandem_muts | |
| 217 | |
| 218 gene = gene[:3] | |
| 219 if gene in ["IGA", "IGG"]: | |
| 220 tandem_sum_by_class[gene] += num_tandem_muts | |
| 221 expected_tandem_sum_by_class[gene] += expected_tandem_muts | |
| 222 else: | |
| 223 tandem_sum_by_class["unmatched"] += num_tandem_muts | |
| 224 expected_tandem_sum_by_class["unmatched"] += expected_tandem_muts | |
| 225 | |
| 226 | |
| 227 for tandem_mut in tandem_muts: | |
| 228 tandem_frequency[str(tandem_mut[1])] += 1 | |
| 229 #print "\t".join([ID, str(len(tandem_muts)), str(longest_tandem[1]) , str(tandem_muts)]) | |
| 230 | |
| 231 tandem_freq_file = os.path.join(os.path.dirname(outfile), "tandem_frequency.txt") | |
| 232 with open(tandem_freq_file, 'w') as o: | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
233 for frq in sorted([int(x) for x in list(tandem_frequency.keys())]): |
| 81 | 234 o.write("{0}\t{1}\n".format(frq, tandem_frequency[str(frq)])) |
| 235 | |
| 236 tandem_row = [] | |
| 237 genes_extra = list(genes) | |
| 238 genes_extra.append("all") | |
| 239 for x, y, in zip([tandem_sum_by_class[x] for x in genes_extra], [expected_tandem_sum_by_class[x] for x in genes_extra]): | |
| 240 if y != 0: | |
| 241 tandem_row += [x, round(y, 2), round(x / y, 2)] | |
| 242 else: | |
| 243 tandem_row += [x, round(y, 2), 0] | |
| 244 | |
| 245 tandem_freq_file = os.path.join(os.path.dirname(outfile), "shm_overview_tandem_row.txt") | |
| 246 with open(tandem_freq_file, 'w') as o: | |
| 247 o.write("Tandems/Expected (ratio),{0}\n".format(",".join([str(x) for x in tandem_row]))) | |
| 248 | |
| 249 #print mutationList, linecount | |
| 250 | |
| 251 AALength = (int(max(mutationList, key=lambda i: int(i[4]) if i[4] and i[5] != ";" else 0)[4]) + 1) # [4] is the position of the AA mutation, None if silent | |
| 252 if AALength < 60: | |
| 253 AALength = 64 | |
| 254 | |
| 255 AA_mutation = [0] * AALength | |
| 256 AA_mutation_dic = {"IGA": AA_mutation[:], "IGG": AA_mutation[:], "IGM": AA_mutation[:], "IGE": AA_mutation[:], "unm": AA_mutation[:], "all": AA_mutation[:]} | |
| 257 AA_mutation_empty = AA_mutation[:] | |
| 258 | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
259 print("AALength:", AALength) |
| 81 | 260 aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/aa_id_mutations.txt" |
| 261 with open(aa_mutations_by_id_file, 'w') as o: | |
| 262 o.write("ID\tbest_match\t" + "\t".join([str(x) for x in range(1,AALength)]) + "\n") | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
263 for ID in list(mutationListByID.keys()): |
| 81 | 264 AA_mutation_for_ID = AA_mutation_empty[:] |
| 265 for mutation in mutationListByID[ID]: | |
| 266 if mutation[4] and mutation[5] != ";": | |
| 267 AA_mutation_position = int(mutation[4]) | |
| 268 try: | |
| 269 AA_mutation[AA_mutation_position] += 1 | |
| 270 AA_mutation_for_ID[AA_mutation_position] += 1 | |
| 271 except Exception as e: | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
272 print(e) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
273 print(mutation) |
| 81 | 274 sys.exit() |
| 275 clss = genedic[ID][:3] | |
| 276 AA_mutation_dic[clss][AA_mutation_position] += 1 | |
| 277 o.write(ID + "\t" + genedic[ID] + "\t" + "\t".join([str(x) for x in AA_mutation_for_ID[1:]]) + "\n") | |
| 278 | |
| 279 | |
| 280 | |
| 281 #absent AA stuff | |
| 282 absentAACDR1Dic = defaultdict(list) | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
283 absentAACDR1Dic[5] = list(range(29,36)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
284 absentAACDR1Dic[6] = list(range(29,35)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
285 absentAACDR1Dic[7] = list(range(30,35)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
286 absentAACDR1Dic[8] = list(range(30,34)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
287 absentAACDR1Dic[9] = list(range(31,34)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
288 absentAACDR1Dic[10] = list(range(31,33)) |
| 81 | 289 absentAACDR1Dic[11] = [32] |
| 290 | |
| 291 absentAACDR2Dic = defaultdict(list) | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
292 absentAACDR2Dic[0] = list(range(55,65)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
293 absentAACDR2Dic[1] = list(range(56,65)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
294 absentAACDR2Dic[2] = list(range(56,64)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
295 absentAACDR2Dic[3] = list(range(57,64)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
296 absentAACDR2Dic[4] = list(range(57,63)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
297 absentAACDR2Dic[5] = list(range(58,63)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
298 absentAACDR2Dic[6] = list(range(58,62)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
299 absentAACDR2Dic[7] = list(range(59,62)) |
|
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
300 absentAACDR2Dic[8] = list(range(59,61)) |
| 81 | 301 absentAACDR2Dic[9] = [60] |
| 302 | |
| 303 absentAA = [len(IDlist)] * (AALength-1) | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
304 for k, cdr1Length in cdr1LengthDic.items(): |
| 81 | 305 for c in absentAACDR1Dic[cdr1Length]: |
| 306 absentAA[c] -= 1 | |
| 307 | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
308 for k, cdr2Length in cdr2LengthDic.items(): |
| 81 | 309 for c in absentAACDR2Dic[cdr2Length]: |
| 310 absentAA[c] -= 1 | |
| 311 | |
| 312 | |
| 313 aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/absent_aa_id.txt" | |
| 314 with open(aa_mutations_by_id_file, 'w') as o: | |
| 315 o.write("ID\tcdr1length\tcdr2length\tbest_match\t" + "\t".join([str(x) for x in range(1,AALength)]) + "\n") | |
| 316 for ID in IDlist: | |
| 317 absentAAbyID = [1] * (AALength-1) | |
| 318 cdr1Length = cdr1LengthDic[ID] | |
| 319 for c in absentAACDR1Dic[cdr1Length]: | |
| 320 absentAAbyID[c] -= 1 | |
| 321 | |
| 322 cdr2Length = cdr2LengthDic[ID] | |
| 323 for c in absentAACDR2Dic[cdr2Length]: | |
| 324 absentAAbyID[c] -= 1 | |
| 325 o.write(ID + "\t" + str(cdr1Length) + "\t" + str(cdr2Length) + "\t" + genedic[ID] + "\t" + "\t".join([str(x) for x in absentAAbyID]) + "\n") | |
| 326 | |
| 327 if linecount == 0: | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
328 print("No data, exiting") |
| 81 | 329 with open(outfile, 'w') as o: |
| 330 o.write("RGYW (%)," + ("0,0,0\n" * len(genes))) | |
| 331 o.write("WRCY (%)," + ("0,0,0\n" * len(genes))) | |
| 332 o.write("WA (%)," + ("0,0,0\n" * len(genes))) | |
| 333 o.write("TW (%)," + ("0,0,0\n" * len(genes))) | |
| 334 sys.exit() | |
| 335 | |
| 336 hotspotMatcher = re.compile("[actg]+,(\d+)-(\d+)\((.*)\)") | |
| 337 RGYWCount = {} | |
| 338 WRCYCount = {} | |
| 339 WACount = {} | |
| 340 TWCount = {} | |
| 341 | |
| 342 #IDIndex = 0 | |
| 343 ataIndex = 0 | |
| 344 tatIndex = 0 | |
| 345 aggctatIndex = 0 | |
| 346 atagcctIndex = 0 | |
| 347 first = True | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
348 with open(infile, 'r') as i: |
| 81 | 349 for line in i: |
| 350 if first: | |
| 351 linesplt = line.split("\t") | |
| 352 ataIndex = linesplt.index("X.a.t.a") | |
| 353 tatIndex = linesplt.index("t.a.t.") | |
| 354 aggctatIndex = linesplt.index("X.a.g.g.c.t..a.t.") | |
| 355 atagcctIndex = linesplt.index("X.a.t..a.g.c.c.t.") | |
| 356 first = False | |
| 357 continue | |
| 358 linesplt = line.split("\t") | |
| 359 gene = linesplt[best_matchIndex] | |
| 360 ID = linesplt[IDIndex] | |
| 361 RGYW = [(int(x), int(y), z) for (x, y, z) in | |
| 362 [hotspotMatcher.match(x).groups() for x in linesplt[aggctatIndex].split("|") if x]] | |
| 363 WRCY = [(int(x), int(y), z) for (x, y, z) in | |
| 364 [hotspotMatcher.match(x).groups() for x in linesplt[atagcctIndex].split("|") if x]] | |
| 365 WA = [(int(x), int(y), z) for (x, y, z) in | |
| 366 [hotspotMatcher.match(x).groups() for x in linesplt[ataIndex].split("|") if x]] | |
| 367 TW = [(int(x), int(y), z) for (x, y, z) in | |
| 368 [hotspotMatcher.match(x).groups() for x in linesplt[tatIndex].split("|") if x]] | |
| 369 RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0, 0, 0, 0 | |
| 370 | |
| 371 with open(os.path.join(os.path.dirname(os.path.abspath(infile)), "RGYW.txt"), 'a') as out_handle: | |
| 372 for hotspot in RGYW: | |
| 373 out_handle.write("{0}\t{1}\n".format(ID, "\t".join([str(x) for x in hotspot]))) | |
| 374 | |
| 375 mutationList = mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"] | |
| 376 for mutation in mutationList: | |
| 377 frm, where, to, AAfrm, AAwhere, AAto, junk = mutation | |
| 378 mutation_in_RGYW = any(((start <= int(where) <= end) for (start, end, region) in RGYW)) | |
| 379 mutation_in_WRCY = any(((start <= int(where) <= end) for (start, end, region) in WRCY)) | |
| 380 mutation_in_WA = any(((start <= int(where) <= end) for (start, end, region) in WA)) | |
| 381 mutation_in_TW = any(((start <= int(where) <= end) for (start, end, region) in TW)) | |
| 382 | |
| 383 in_how_many_motifs = sum([mutation_in_RGYW, mutation_in_WRCY, mutation_in_WA, mutation_in_TW]) | |
| 384 | |
| 385 if in_how_many_motifs > 0: | |
| 386 RGYWCount[ID] += (1.0 * int(mutation_in_RGYW)) / in_how_many_motifs | |
| 387 WRCYCount[ID] += (1.0 * int(mutation_in_WRCY)) / in_how_many_motifs | |
| 388 WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs | |
| 389 TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs | |
| 390 | |
| 391 mutations_in_motifs_file = os.path.join(os.path.dirname(os.path.abspath(infile)), "mutation_in_motifs.txt") | |
| 392 if not os.path.exists(mutation_by_id_file): | |
| 393 with open(mutations_in_motifs_file, 'w') as out_handle: | |
| 394 out_handle.write("{0}\n".format("\t".join([ | |
| 395 "Sequence.ID", | |
| 396 "mutation_position", | |
| 397 "region", | |
| 398 "from_nt", | |
| 399 "to_nt", | |
| 400 "mutation_position_AA", | |
| 401 "from_AA", | |
| 402 "to_AA", | |
| 403 "motif", | |
| 404 "motif_start_nt", | |
| 405 "motif_end_nt", | |
| 406 "rest" | |
| 407 ]))) | |
| 408 | |
| 409 with open(mutations_in_motifs_file, 'a') as out_handle: | |
| 410 motif_dic = {"RGYW": RGYW, "WRCY": WRCY, "WA": WA, "TW": TW} | |
| 411 for mutation in mutationList: | |
| 412 frm, where, to, AAfrm, AAwhere, AAto, junk = mutation | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
413 for motif in list(motif_dic.keys()): |
| 81 | 414 |
| 415 for start, end, region in motif_dic[motif]: | |
| 416 if start <= int(where) <= end: | |
| 417 out_handle.write("{0}\n".format( | |
| 418 "\t".join([ | |
| 419 ID, | |
| 420 where, | |
| 421 region, | |
| 422 frm, | |
| 423 to, | |
| 424 str(AAwhere), | |
| 425 str(AAfrm), | |
| 426 str(AAto), | |
| 427 motif, | |
| 428 str(start), | |
| 429 str(end), | |
| 430 str(junk) | |
| 431 ]) | |
| 432 )) | |
| 433 | |
| 434 | |
| 435 | |
| 436 def mean(lst): | |
| 437 return (float(sum(lst)) / len(lst)) if len(lst) > 0 else 0.0 | |
| 438 | |
| 439 | |
| 440 def median(lst): | |
| 441 lst = sorted(lst) | |
| 442 l = len(lst) | |
| 443 if l == 0: | |
| 444 return 0 | |
| 445 if l == 1: | |
| 446 return lst[0] | |
| 447 | |
| 448 l = int(l / 2) | |
| 449 | |
| 450 if len(lst) % 2 == 0: | |
| 451 return float(lst[l] + lst[(l - 1)]) / 2.0 | |
| 452 else: | |
| 453 return lst[l] | |
| 454 | |
| 455 funcs = {"mean": mean, "median": median, "sum": sum} | |
| 456 | |
| 457 directory = outfile[:outfile.rfind("/") + 1] | |
| 458 value = 0 | |
| 459 valuedic = dict() | |
| 460 | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
461 for fname in list(funcs.keys()): |
| 81 | 462 for gene in genes: |
| 463 with open(directory + gene + "_" + fname + "_value.txt", 'r') as v: | |
| 464 valuedic[gene + "_" + fname] = float(v.readlines()[0].rstrip()) | |
| 465 with open(directory + "all_" + fname + "_value.txt", 'r') as v: | |
| 466 valuedic["total_" + fname] = float(v.readlines()[0].rstrip()) | |
| 467 | |
| 468 | |
| 469 def get_xyz(lst, gene, f, fname): | |
| 470 x = round(round(f(lst), 1)) | |
| 471 y = valuedic[gene + "_" + fname] | |
| 472 z = str(round(x / float(y) * 100, 1)) if y != 0 else "0" | |
| 473 return (str(x), str(y), z) | |
| 474 | |
| 475 dic = {"RGYW": RGYWCount, "WRCY": WRCYCount, "WA": WACount, "TW": TWCount} | |
| 476 arr = ["RGYW", "WRCY", "WA", "TW"] | |
| 477 | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
478 for fname in list(funcs.keys()): |
| 81 | 479 func = funcs[fname] |
| 480 foutfile = outfile[:outfile.rindex("/")] + "/hotspot_analysis_" + fname + ".txt" | |
| 481 with open(foutfile, 'w') as o: | |
| 482 for typ in arr: | |
| 483 o.write(typ + " (%)") | |
| 484 curr = dic[typ] | |
| 485 for gene in genes: | |
| 486 geneMatcher = geneMatchers[gene] | |
| 487 if valuedic[gene + "_" + fname] is 0: | |
| 488 o.write(",0,0,0") | |
| 489 else: | |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
490 x, y, z = get_xyz([curr[x] for x in [y for y, z in genedic.items() if geneMatcher.match(z)]], gene, func, fname) |
| 81 | 491 o.write("," + x + "," + y + "," + z) |
|
83
729738462297
"planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents:
81
diff
changeset
|
492 x, y, z = get_xyz([y for x, y in curr.items() if not genedic[x].startswith("unmatched")], "total", func, fname) |
| 81 | 493 #x, y, z = get_xyz([y for x, y in curr.iteritems()], "total", func, fname) |
| 494 o.write("," + x + "," + y + "," + z + "\n") | |
| 495 | |
| 496 | |
| 497 # for testing | |
| 498 seq_motif_file = outfile[:outfile.rindex("/")] + "/motif_per_seq.txt" | |
| 499 with open(seq_motif_file, 'w') as o: | |
| 500 o.write("ID\tRGYW\tWRCY\tWA\tTW\n") | |
| 501 for ID in IDlist: | |
| 502 #o.write(ID + "\t" + str(round(RGYWCount[ID], 2)) + "\t" + str(round(WRCYCount[ID], 2)) + "\t" + str(round(WACount[ID], 2)) + "\t" + str(round(TWCount[ID], 2)) + "\n") | |
| 503 o.write(ID + "\t" + str(RGYWCount[ID]) + "\t" + str(WRCYCount[ID]) + "\t" + str(WACount[ID]) + "\t" + str(TWCount[ID]) + "\n") | |
| 504 | |
| 505 if __name__ == "__main__": | |
| 506 main() |