annotate vcf_gff.py @ 1:a0689dc29b7f draft

Updated vcf to gff conversion tool
author john-mccallum
date Tue, 31 Jul 2012 00:33:11 -0400
parents
children 402c3f0fe807
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
1 #!/usr/local/bin/python2.6
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
2 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
3 Usage vcf_gff.py <vcf_file input> <gff_file output>
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
4
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
5 Input vcf file format:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
6 CHROM POS ID REF ALT QUAL FILTER INFO FORMAT
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
7
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
8 Note: Generating vcf from a single merged bam file, using multiple bam files results in multiple FORMAT columns!!!
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
9
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
10 Output gff format:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
11 SEQID SOURCE TYPE START END SCORE STRAND PHASE ATTRIBUTES
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
12
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
13 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
14 #Copyright 2012 Susan Thomson
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
15 #New Zealand Institute for Plant and Food Research
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
16 #This program is free software: you can redistribute it and/or modify
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
17 # it under the terms of the GNU General Public License as published by
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
18 # the Free Software Foundation, either version 3 of the License, or
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
19 # (at your option) any later version.
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
20 #
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
21 # This program is distributed in the hope that it will be useful,
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
22 # but WITHOUT ANY WARRANTY; without even the implied warranty of
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
23 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
24 # GNU General Public License for more details.
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
25 #
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
26 # You should have received a copy of the GNU General Public License
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
27 # along with this program. If not, see <http://www.gnu.org/licenses/>.
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
28
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
29
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
30 import sys
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
31
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
32 in_vcf_file = open(sys.argv[1], 'r')
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
33 out_gff_file = open(sys.argv[2], 'w')
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
34
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
35 def get_info(attribute_input):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
36 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
37 Get the record_type by determining if INDEL is stated in column INFO; get
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
38 raw read count
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
39 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
40 INFO = {}
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
41 rec = attribute_input.split(";")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
42 if "INDEL" in rec:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
43 record_type = "INDEL"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
44 rec = rec[1:]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
45 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
46 record_type = "SNP"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
47 for entry in rec:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
48 detail = entry.split("=")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
49 INFO[detail[0]] = detail[1]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
50 if INFO.has_key("DP"):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
51 reads = INFO.get("DP")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
52 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
53 reads = "NA"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
54 data = (record_type, reads)
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
55 return data
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
56
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
57 def get_gen(formatcols, ref):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
58 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
59 Get info on heterozyosity or homozygosity (could be useful later),
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
60 by estimating genotype(GT) calling ref allele=0, variant allele=1
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
61
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
62 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
63 formats = []
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
64 sample_dict = {}
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
65 genos = ""
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
66 format_types = formatcols[0].split(":")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
67 samples = formatcols[1:]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
68 for entry in format_types:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
69 formats.append(entry)
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
70 for sample in samples:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
71 var = ""
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
72 data = sample.split(":")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
73 sample_dict = dict(zip(formats, data))
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
74 if sample_dict.has_key("DP"):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
75 reads = sample_dict.get("DP")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
76 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
77 reads = "NA"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
78 if sample_dict.has_key("NF"):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
79 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
80 polysnp tool output
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
81 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
82 freq = sample_dict.get("NF")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
83 variants = freq.split("|")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
84 if int(variants[0]) >= 1:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
85 var += "A"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
86 if int(variants[1]) >= 1:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
87 var += "C"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
88 if int(variants[2]) >= 1:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
89 var += "G"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
90 if int(variants[3]) >= 1:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
91 var += "T"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
92 if var == "":
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
93 gen = "NA"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
94 elif var == ref:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
95 gen = "HOM_ref"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
96 elif len(var) >= 2:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
97 gen = "HET"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
98 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
99 gen = "HOM_mut"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
100 elif sample_dict.has_key("GT"):
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
101 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
102 mpileup output, recommend ALWAYS have GT, note this only good scoring for diploids too!
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
103 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
104 genotypes = sample_dict.get("GT")
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
105 if genotypes == "1/1":
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
106 gen = "HOM_mut"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
107 if genotypes == "0/1":
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
108 gen = "HET"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
109 if genotypes == "0/0":
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
110 gen = "HOM_ref"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
111 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
112 gen = "NA"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
113 geno = ("%s:%s;" % (reads, gen))
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
114 genos += geno
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
115 sample_dict = {}
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
116 return genos
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
117
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
118 attributes = {}
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
119 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
120 Get relevant info from vcf file and put to proper gff columns
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
121 """
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
122
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
123 for line in in_vcf_file:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
124 if line.startswith("#") == False:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
125 info = line.split()
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
126 seqid = info[0].strip()
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
127 source = "SAMTOOLS"
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
128 start = int(info[1])
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
129 score = info[5]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
130 strand = "."
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
131 phase = "."
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
132 reference = info[3].strip()
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
133 variant = info[4].strip()
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
134 attr = get_info(info[7])
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
135 record_type = attr[0]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
136 reads = attr[1]
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
137 if record_type == "INDEL":
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
138 end = start + len(reference)
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
139 else:
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
140 end = start
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
141 gen = get_gen(info[8:], reference)
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
142 out_gff_file.write(
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
143 ("%s\t%s\t%s\t%d\t%d\t%s\t%s\t%s\tID=%s:%s:%d;Variant" +
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
144 "_seq=%s;Reference_seq=%s;Total_reads=%s:Zygosity=%s\n") %
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
145 ( seqid, source,record_type, start, end, score, strand, phase,seqid,
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
146 record_type, start, variant, reference, reads, gen))
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
147
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
148 out_gff_file.close()
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
149
a0689dc29b7f Updated vcf to gff conversion tool
john-mccallum
parents:
diff changeset
150