|
12
|
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() |
