Mercurial > repos > iuc > medaka_consensus_pipeline

comparison annotateVCF.py @ 5:b35de691c42c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/medaka commit 13769e7d51b30a1d15eb62a9ba89ee2064f3ddc3"
author iuc
date Fri, 16 Oct 2020 17:04:09 +0000
parents 35666d44fe7a
children
comparison
equal deleted inserted replaced
4:a1b70f038b4a 5:b35de691c42c
73 if ins_flag: 73 if ins_flag:
74 if base.isdigit(): 74 if base.isdigit():
75 ins_str += base 75 ins_str += base
76 else: 76 else:
77 ins_len = int(ins_str) - 1 77 ins_len = int(ins_str) - 1
78 ins_str = ""
78 insertion = base 79 insertion = base
79 ins_flag = False 80 ins_flag = False
80 elif del_flag: 81 elif del_flag:
81 if base.isdigit(): 82 if base.isdigit():
82 del_str += base 83 del_str += base
83 else: 84 else:
84 del_len = int(del_str) - 1 85 del_len = int(del_str) - 1
86 del_str = ""
85 deletion = base 87 deletion = base
86 del_flag = False 88 del_flag = False
87 else: 89 else:
88 if base == '^': 90 if base == '^':
89 carrot_flag = True 91 carrot_flag = True
105 elif base == '*': 107 elif base == '*':
106 continue 108 continue
107 else: 109 else:
108 counts[base_to_idx[base]] += 1 110 counts[base_to_idx[base]] += 1
109 stranded_counts[base_to_idx_stranded[base]] += 1 111 stranded_counts[base_to_idx_stranded[base]] += 1
110 af = float(counts[base_to_idx[alt_base]]) / float(sum(counts)) 112 if sum(counts) == 0:
113 af = float("nan")
114 else:
115 af = float(counts[base_to_idx[alt_base]]) / float(sum(counts))
111 if float(sum(stranded_counts[0:4])) == 0: 116 if float(sum(stranded_counts[0:4])) == 0:
112 faf = float("nan") 117 faf = float("nan")
113 else: 118 else:
114 faf = float(stranded_counts[base_to_idx_stranded[alt_base]]) / float(sum(stranded_counts[0:4])) 119 faf = float(stranded_counts[base_to_idx_stranded[alt_base]]) / float(sum(stranded_counts[0:4]))
115 if float(sum(stranded_counts[4:])) == 0: 120 if float(sum(stranded_counts[4:])) == 0:
291 if forward_flag: 296 if forward_flag:
292 counts[7] += 1 297 counts[7] += 1
293 else: 298 else:
294 counts[8] += 1 299 counts[8] += 1
295 dp = int(fields[3]) 300 dp = int(fields[3])
296 af = float(counts[3]) / float(sum([counts[0], counts[3], counts[6]])) 301 if sum([counts[0], counts[3], counts[6]]) == 0:
302 af = float("nan")
303 else:
304 af = float(counts[3]) / float(sum([counts[0], counts[3], counts[6]]))
297 if sum([counts[1], counts[4], counts[7]]) == 0: 305 if sum([counts[1], counts[4], counts[7]]) == 0:
298 faf = float("nan") 306 faf = float("nan")
299 else: 307 else:
300 faf = float(counts[4]) / float(sum([counts[1], counts[4], counts[7]])) 308 faf = float(counts[4]) / float(sum([counts[1], counts[4], counts[7]]))
301 if sum([counts[2], counts[5], counts[8]]) == 0: 309 if sum([counts[2], counts[5], counts[8]]) == 0:
315 to_examine = {} 323 to_examine = {}
316 for line in in_vcf: 324 for line in in_vcf:
317 if line[0:2] == "##": 325 if line[0:2] == "##":
318 out_vcf.write(line) 326 out_vcf.write(line)
319 elif line[0] == "#": 327 elif line[0] == "#":
320 out_vcf.write("##annotateVCFVersion=0.1\n") 328 out_vcf.write("##annotateVCFVersion=0.2\n")
321 out_vcf.write("##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Raw Depth\">\n") 329 out_vcf.write("##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Raw Depth\">\n")
322 out_vcf.write("##INFO=<ID=AF,Number=1,Type=Float,Description=\"Allele Frequency\">\n") 330 out_vcf.write("##INFO=<ID=AF,Number=1,Type=Float,Description=\"Allele Frequency\">\n")
323 out_vcf.write("##INFO=<ID=FAF,Number=1,Type=Float,Description=\"Forward Allele Frequency\">\n") 331 out_vcf.write("##INFO=<ID=FAF,Number=1,Type=Float,Description=\"Forward Allele Frequency\">\n")
324 out_vcf.write("##INFO=<ID=RAF,Number=1,Type=Float,Description=\"Reverse Allele Frequency\">\n") 332 out_vcf.write("##INFO=<ID=RAF,Number=1,Type=Float,Description=\"Reverse Allele Frequency\">\n")
325 out_vcf.write("##INFO=<ID=SB,Number=1,Type=Integer,Description=\"Phred-scaled strand bias at this position\">\n") 333 out_vcf.write("##INFO=<ID=SB,Number=1,Type=Integer,Description=\"Phred-scaled strand bias at this position\">\n")
374 if fields[7] == "": 382 if fields[7] == "":
375 info = [] 383 info = []
376 else: 384 else:
377 info = fields[7].split(';') 385 info = fields[7].split(';')
378 info.append("DP=%d" % (dp)) 386 info.append("DP=%d" % (dp))
379 info.append("AF=%.6f" % (af)) 387 if isnan(af):
388 info.append("AF=NaN")
389 else:
390 info.append("AF=%.6f" % (af))
380 if isnan(faf): 391 if isnan(faf):
381 info.append("FAF=NaN") 392 info.append("FAF=NaN")
382 else: 393 else:
383 info.append("FAF=%.6f" % (faf)) 394 info.append("FAF=%.6f" % (faf))
384 if isnan(raf): 395 if isnan(raf):