annotate convert_VCF_info_fields.py @ 7:28c13c42de01 draft

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