Mercurial > repos > in_silico > cravat_score_and_annotate
comparison cravat_convert/vcf_converter.py @ 21:67a13940d0bf draft default tip
Uploaded
author | in_silico |
---|---|
date | Thu, 16 Aug 2018 15:10:43 -0400 |
parents | 6da9a12f04b5 |
children |
comparison
equal
deleted
inserted
replaced
20:6da9a12f04b5 | 21:67a13940d0bf |
---|---|
1 """ | |
2 A module originally obtained from the cravat package. Modified to use in the vcf | |
3 converter galaxy tool. | |
4 | |
5 | |
6 Register of changes made (Chris Jacoby): | |
7 1) Changed imports as galaxy tool won't have access to complete cravat python package | |
8 2) Defined BadFormatError in BaseConverted file, as I didn't have the BadFormatError module | |
9 """ | |
10 | |
11 from base_converter import BaseConverter, BadFormatError | |
12 import re | |
13 | |
14 class CravatConverter(BaseConverter): | |
15 | |
16 def __init__(self): | |
17 self.format_name = 'vcf' | |
18 self.samples = [] | |
19 self.var_counter = 0 | |
20 self.addl_cols = [{'name':'phred', | |
21 'title':'Phred', | |
22 'type':'string'}, | |
23 {'name':'filter', | |
24 'title':'VCF filter', | |
25 'type':'string'}, | |
26 {'name':'zygosity', | |
27 'title':'Zygosity', | |
28 'type':'string'}, | |
29 {'name':'alt_reads', | |
30 'title':'Alternate reads', | |
31 'type':'int'}, | |
32 {'name':'tot_reads', | |
33 'title':'Total reads', | |
34 'type':'int'}, | |
35 {'name':'af', | |
36 'title':'Variant allele frequency', | |
37 'type':'float'}] | |
38 | |
39 def check_format(self, f): | |
40 return f.readline().startswith('##fileformat=VCF') | |
41 | |
42 def setup(self, f): | |
43 | |
44 vcf_line_no = 0 | |
45 for line in f: | |
46 vcf_line_no += 1 | |
47 if len(line) < 6: | |
48 continue | |
49 if line[:6] == '#CHROM': | |
50 toks = re.split('\s+', line.rstrip()) | |
51 if len(toks) > 8: | |
52 self.samples = toks[9:] | |
53 break | |
54 | |
55 def convert_line(self, l): | |
56 if l.startswith('#'): return None | |
57 self.var_counter += 1 | |
58 toks = l.strip('\r\n').split('\t') | |
59 all_wdicts = [] | |
60 if len(toks) < 8: | |
61 raise BadFormatError('Wrong VCF format') | |
62 [chrom, pos, tag, ref, alts, qual, filter, info] = toks[:8] | |
63 if tag == '': | |
64 raise BadFormatError('ID column is blank') | |
65 elif tag == '.': | |
66 tag = 'VAR' + str(self.var_counter) | |
67 if chrom[:3] != 'chr': | |
68 chrom = 'chr' + chrom | |
69 alts = alts.split(',') | |
70 len_alts = len(alts) | |
71 if len(toks) == 8: | |
72 for altno in range(len_alts): | |
73 wdict = None | |
74 alt = alts[altno] | |
75 newpos, newref, newalt = self.extract_vcf_variant('+', pos, ref, alt) | |
76 wdict = {'tags':tag, | |
77 'chrom':chrom, | |
78 'pos':newpos, | |
79 'ref_base':newref, | |
80 'alt_base':newalt, | |
81 'sample_id':'no_sample', | |
82 'phred': qual, | |
83 'filter': filter} | |
84 all_wdicts.append(wdict) | |
85 elif len(toks) > 8: | |
86 sample_datas = toks[9:] | |
87 genotype_fields = {} | |
88 genotype_field_no = 0 | |
89 for genotype_field in toks[8].split(':'): | |
90 genotype_fields[genotype_field] = genotype_field_no | |
91 genotype_field_no += 1 | |
92 if not ('GT' in genotype_fields): | |
93 raise BadFormatError('No GT Field') | |
94 gt_field_no = genotype_fields['GT'] | |
95 for sample_no in range(len(sample_datas)): | |
96 sample = self.samples[sample_no] | |
97 sample_data = sample_datas[sample_no].split(':') | |
98 gts = {} | |
99 for gt in sample_data[gt_field_no].replace('/', '|').split('|'): | |
100 if gt == '.': | |
101 continue | |
102 else: | |
103 gts[int(gt)] = True | |
104 for gt in sorted(gts.keys()): | |
105 wdict = None | |
106 if gt == 0: | |
107 continue | |
108 else: | |
109 alt = alts[gt - 1] | |
110 newpos, newref, newalt = self.extract_vcf_variant('+', pos, ref, alt) | |
111 zyg = self.homo_hetro(sample_data[gt_field_no]) | |
112 depth, alt_reads, af = self.extract_read_info(sample_data, gt, gts, genotype_fields) | |
113 | |
114 wdict = {'tags':tag, | |
115 'chrom':chrom, | |
116 'pos':newpos, | |
117 'ref_base':newref, | |
118 'alt_base':newalt, | |
119 'sample_id':sample, | |
120 'phred': qual, | |
121 'filter': filter, | |
122 'zygosity': zyg, | |
123 'tot_reads': depth, | |
124 'alt_reads': alt_reads, | |
125 'af': af, | |
126 } | |
127 all_wdicts.append(wdict) | |
128 return all_wdicts | |
129 | |
130 #The vcf genotype string has a call for each allele separated by '\' or '/' | |
131 #If the call is the same for all allels, return 'hom' otherwise 'het' | |
132 def homo_hetro(self, gt_str): | |
133 if '.' in gt_str: | |
134 return ''; | |
135 | |
136 gts = gt_str.strip().replace('/', '|').split('|') | |
137 for gt in gts: | |
138 if gt != gts[0]: | |
139 return 'het' | |
140 return 'hom' | |
141 | |
142 #Extract read depth, allele count, and allele frequency from optional VCR information | |
143 def extract_read_info (self, sample_data, gt, gts, genotype_fields): | |
144 depth = '' | |
145 alt_reads = '' | |
146 ref_reads = '' | |
147 af = '' | |
148 | |
149 #AD contains 2 values usually ref count and alt count unless there are | |
150 #multiple alts then it will have alt 1 then alt 2. | |
151 if 'AD' in genotype_fields and genotype_fields['AD'] <= len(sample_data): | |
152 if 0 in gts.keys(): | |
153 #if part of the genotype is reference, then AD will have #ref reads, #alt reads | |
154 ref_reads = sample_data[genotype_fields['AD']].split(',')[0] | |
155 alt_reads = sample_data[genotype_fields['AD']].split(',')[1] | |
156 elif gt == max(gts.keys()): | |
157 #if geontype has multiple alt bases, then AD will have #alt1 reads, #alt2 reads | |
158 alt_reads = sample_data[genotype_fields['AD']].split(',')[1] | |
159 else: | |
160 alt_reads = sample_data[genotype_fields['AD']].split(',')[0] | |
161 | |
162 if 'DP' in genotype_fields and genotype_fields['DP'] <= len(sample_data): | |
163 depth = sample_data[genotype_fields['DP']] | |
164 elif alt_reads != '' and ref_reads != '': | |
165 #if DP is not present but we have alt and ref reads count, dp = ref+alt | |
166 depth = int(alt_reads) + int(ref_reads) | |
167 | |
168 if 'AF' in genotype_fields and genotype_fields['AF'] <= len(sample_data): | |
169 af = float(sample_data[genotype_fields['AF']] ) | |
170 elif depth != '' and alt_reads != '': | |
171 #if AF not specified, calc it from alt and ref reads | |
172 af = float(alt_reads) / float(depth) | |
173 | |
174 return depth, alt_reads, af | |
175 | |
176 def extract_vcf_variant (self, strand, pos, ref, alt): | |
177 | |
178 reflen = len(ref) | |
179 altlen = len(alt) | |
180 | |
181 # Returns without change if same single nucleotide for ref and alt. | |
182 if reflen == 1 and altlen == 1 and ref == alt: | |
183 return pos, ref, alt | |
184 | |
185 # Trimming from the start and then the end of the sequence | |
186 # where the sequences overlap with the same nucleotides | |
187 new_ref2, new_alt2, new_pos = \ | |
188 self.trimming_vcf_input(ref, alt, pos, strand) | |
189 | |
190 if new_ref2 == '': | |
191 new_ref2 = '-' | |
192 if new_alt2 == '': | |
193 new_alt2 = '-' | |
194 | |
195 return new_pos, new_ref2, new_alt2 | |
196 | |
197 # This function looks at the ref and alt sequences and removes | |
198 # where the overlapping sequences contain the same nucleotide. | |
199 # This trims from the end first but does not remove the first nucleotide | |
200 # because based on the format of VCF input the | |
201 # first nucleotide of the ref and alt sequence occur | |
202 # at the position specified. | |
203 # End removed first, not the first nucleotide | |
204 # Front removed and position changed | |
205 def trimming_vcf_input(self, ref, alt, pos, strand): | |
206 pos = int(pos) | |
207 reflen = len(ref) | |
208 altlen = len(alt) | |
209 minlen = min(reflen, altlen) | |
210 new_ref = ref | |
211 new_alt = alt | |
212 new_pos = pos | |
213 # Trims from the end. Except don't remove the first nucleotide. | |
214 # 1:6530968 CTCA -> GTCTCA becomes C -> GTC. | |
215 for nt_pos in range(0, minlen - 1): | |
216 if ref[reflen - nt_pos - 1] == alt[altlen - nt_pos - 1]: | |
217 new_ref = ref[:reflen - nt_pos - 1] | |
218 new_alt = alt[:altlen - nt_pos - 1] | |
219 else: | |
220 break | |
221 new_ref_len = len(new_ref) | |
222 new_alt_len = len(new_alt) | |
223 minlen = min(new_ref_len, new_alt_len) | |
224 new_ref2 = new_ref | |
225 new_alt2 = new_alt | |
226 # Trims from the start. 1:6530968 G -> GT becomes 1:6530969 - -> T. | |
227 for nt_pos in range(0, minlen): | |
228 if new_ref[nt_pos] == new_alt[nt_pos]: | |
229 if strand == '+': | |
230 new_pos += 1 | |
231 elif strand == '-': | |
232 new_pos -= 1 | |
233 new_ref2 = new_ref[nt_pos + 1:] | |
234 new_alt2 = new_alt[nt_pos + 1:] | |
235 else: | |
236 new_ref2 = new_ref[nt_pos:] | |
237 new_alt2 = new_alt[nt_pos:] | |
238 break | |
239 return new_ref2, new_alt2, new_pos | |
240 | |
241 | |
242 if __name__ == "__main__": | |
243 c = CravatConverter() |