Mercurial > repos > iuc > medaka_consensus_pipeline

comparison convert_VCF_info_fields.py @ 6:d915879a88bc 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 9b7d28ac59ad082874670ee989836631ba8d7fb4"
author iuc
date Wed, 10 Feb 2021 08:28:46 +0000
parents
children 57beb8ea55b4
comparison
equal deleted inserted replaced
5:b35de691c42c 6:d915879a88bc
1 #!/usr/bin/env python3
2
3 # Takes in VCF file annotated with medaka tools annotate and converts
4 #
5 # Usage statement:
6 # python convert_VCF_info_fields.py in_vcf.vcf out_vcf.vcf
7
8 # 10/21/2020 - Nathan P. Roach, natproach@gmail.com
9
10 import sys
11 from collections import OrderedDict
12 from math import log10
13
14 from scipy.stats import fisher_exact
15
16
17 def pval_to_phredqual(pval):
18 return round(-10 * log10(pval))
19
20
21 def parseInfoField(info):
22 info_fields = info.split(';')
23 info_dict = OrderedDict()
24 for info_field in info_fields:
25 code, val = info_field.split('=')
26 info_dict[code] = val
27 return info_dict
28
29
30 def annotateVCF(in_vcf_filepath, out_vcf_filepath):
31 in_vcf = open(in_vcf_filepath, 'r')
32 out_vcf = open(out_vcf_filepath, 'w')
33 to_skip = set(['SC', 'SR'])
34 for i, line in enumerate(in_vcf):
35 if i == 1:
36 out_vcf.write("##convert_VCF_info_fields=0.1\n")
37 if line[0:2] == "##":
38 if line[0:11] == "##INFO=<ID=":
39 id_ = line[11:].split(',')[0]
40 if id_ in to_skip:
41 continue
42 out_vcf.write(line)
43 elif line[0] == "#":
44 out_vcf.write('##INFO=<ID=DPSPS,Number=2,Type=Integer,Description="Spanning Reads Allele Frequency By Strand">\n')
45 out_vcf.write('##INFO=<ID=AF,Number=1,Type=Float,Description="Spanning Reads Allele Frequency">\n')
46 out_vcf.write('##INFO=<ID=FAF,Number=1,Type=Float,Description="Forward Spanning Reads Allele Frequency">\n')
47 out_vcf.write('##INFO=<ID=RAF,Number=1,Type=Float,Description="Reverse Spanning Reads Allele Frequency">\n')
48 out_vcf.write('##INFO=<ID=SB,Number=1,Type=Integer,Description="Phred-scaled strand bias of spanning reads at this position">\n')
49 out_vcf.write('##INFO=<ID=DP4,Number=4,Type=Integer,Description="Counts for ref-forward bases, ref-reverse, alt-forward and alt-reverse bases in spanning reads">\n')
50 out_vcf.write('##INFO=<ID=AS,Number=4,Type=Integer,Description="Total alignment score to ref and alt allele of spanning reads by strand (ref fwd, ref rev, alt fwd, alt rev) aligned with parasail match 5, mismatch -4, open 5, extend 3">\n')
51 out_vcf.write(line)
52 else:
53 fields = line.split('\t')
54 info_dict = parseInfoField(fields[7])
55 sr_list = [int(x) for x in info_dict["SR"].split(',')]
56 sc_list = [int(x) for x in info_dict["SC"].split(',')]
57 if len(sr_list) == len(sc_list):
58 variant_list = fields[4].split(',')
59 dpsp = int(info_dict["DPSP"])
60 ref_fwd, ref_rev = 0, 1
61 dpspf, dpspr = (int(x) for x in info_dict["AR"].split(','))
62 for i in range(0, len(sr_list), 2):
63 dpspf += sr_list[i]
64 dpspr += sr_list[i + 1]
65 for j, i in enumerate(range(2, len(sr_list), 2)):
66 dp4 = (sr_list[ref_fwd], sr_list[ref_rev], sr_list[i], sr_list[i + 1])
67 dp2x2 = [[dp4[0], dp4[1]], [dp4[2], dp4[3]]]
68 _, p_val = fisher_exact(dp2x2)
69 sb = pval_to_phredqual(p_val)
70
71 as_ = (sc_list[ref_fwd], sc_list[ref_rev], sc_list[i], sc_list[i + 1])
72
73 info = []
74 for code in info_dict:
75 if code in to_skip:
76 continue
77 val = info_dict[code]
78 info.append("%s=%s" % (code, val))
79
80 info.append("DPSPS=%d,%d" % (dpspf, dpspr))
81
82 if dpsp == 0:
83 info.append("AF=NaN")
84 else:
85 af = dp4[2] + dp4[3] / dpsp
86 info.append("AF=%.6f" % (af))
87 if dpspf == 0:
88 info.append("FAF=NaN")
89 else:
90 faf = dp4[2] / dpspf
91 info.append("FAF=%.6f" % (faf))
92 if dpspr == 0:
93 info.append("RAF=NaN")
94 else:
95 raf = dp4[3] / dpspr
96 info.append("RAF=%.6f" % (raf))
97 info.append("SB=%d" % (sb))
98 info.append("DP4=%d,%d,%d,%d" % (dp4))
99 info.append("AS=%d,%d,%d,%d" % (as_))
100 new_info = ';'.join(info)
101 fields[4] = variant_list[j]
102 fields[7] = new_info
103 out_vcf.write("%s" % ("\t".join(fields)))
104 else:
105 print("WARNING - SR and SC are different lengths, skipping variant")
106 print(line.strip()) # Print the line for debugging purposes
107 in_vcf.close()
108 out_vcf.close()
109
110
111 if __name__ == "__main__":
112 annotateVCF(sys.argv[1], sys.argv[2])