Mercurial > repos > gregory-minevich > cloudmap_variant_discovery_mapping
comparison VDM_Mapping.py @ 9:3dba6603108e draft
Uploaded
author | gregory-minevich |
---|---|
date | Fri, 09 May 2014 16:55:03 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
8:985aacf4d157 | 9:3dba6603108e |
---|---|
1 #!/usr/bin/python | |
2 | |
3 import re | |
4 import sys | |
5 import optparse | |
6 import csv | |
7 import re | |
8 import pprint | |
9 from decimal import * | |
10 from rpy import * | |
11 | |
12 def main(): | |
13 csv.field_size_limit(1000000000) | |
14 | |
15 parser = optparse.OptionParser() | |
16 parser.add_option('-v', '--sample_vcf', dest = 'sample_vcf', action = 'store', type = 'string', default = None, help = "Sample VCF from GATK Unified Genotyper") | |
17 parser.add_option('-l', '--loess_span', dest = 'loess_span', action = 'store', type = 'float', default = .1, help = "Loess span") | |
18 parser.add_option('-d', '--d_yaxis', dest = 'd_yaxis', action = 'store', type = 'float', default = 1, help = "y-axis upper limit for dot plot") | |
19 parser.add_option('-y', '--h_yaxis', dest = 'h_yaxis', action = 'store', type = 'int', default = 0, help = "y-axis upper limit for histogram plot") | |
20 parser.add_option('-c', '--points_color', dest = 'points_color', action = 'store', type = 'string', default = "gray27", help = "Color for data points") | |
21 parser.add_option('-k', '--loess_color', dest = 'loess_color', action = 'store', type = 'string', default = "green2", help = "Color for loess regression line") | |
22 parser.add_option('-z', '--standardize', dest = 'standardize', default= 'true', help = "Standardize X-axis") | |
23 parser.add_option('-b', '--break_file', dest = 'break_file', action = 'store', type = 'string', default = 'C.elegans', help = "File defining the breaks per chromosome") | |
24 parser.add_option('-x', '--bin_size', dest = 'bin_size', action = 'store', type = 'int', default = 1000000, help = "Size of histogram bins, default is 1mb") | |
25 parser.add_option('-n', '--normalize_bins', dest = 'normalize_bins', default= 'true', help = "Normalize histograms") | |
26 | |
27 parser.add_option('-o', '--output', dest = 'output', action = 'store', type = 'string', default = None, help = "Output file name") | |
28 parser.add_option('-s', '--location_plot_output', dest = 'location_plot_output', action = 'store', type = 'string', default = "SNP_Mapping_Plot.pdf", help = "Output file name of SNP plots by chromosomal location") | |
29 | |
30 (options, args) = parser.parse_args() | |
31 | |
32 vcf_info = parse_vcf(sample_vcf = options.sample_vcf) | |
33 | |
34 output_vcf_info(output = options.output, vcf_info = vcf_info) | |
35 | |
36 #output plot with all ratios | |
37 rounded_bin_size = int(round((float(options.bin_size) / 1000000), 1) * 1000000) | |
38 | |
39 normalized_histogram_bins_per_mb = calculate_normalized_histogram_bins_per_xbase(vcf_info = vcf_info, xbase = rounded_bin_size, normalize_bins = options.normalize_bins) | |
40 max_y_hist_mb = normalized_histogram_bins_per_mb[max(normalized_histogram_bins_per_mb, key = lambda x: normalized_histogram_bins_per_mb.get(x) )] | |
41 | |
42 normalized_histogram_bins_per_5kb = calculate_normalized_histogram_bins_per_xbase(vcf_info = vcf_info, xbase = (rounded_bin_size / 2), normalize_bins = options.normalize_bins) | |
43 max_y_hist_5kb = normalized_histogram_bins_per_5kb[max(normalized_histogram_bins_per_5kb, key = lambda x: normalized_histogram_bins_per_5kb.get(x) )] | |
44 | |
45 max_y_hist_overall = myround(max(max_y_hist_mb, max_y_hist_5kb) + int(round(round(max(max_y_hist_mb, max_y_hist_5kb)) * .1))) | |
46 | |
47 | |
48 break_dict = parse_breaks(break_file = options.break_file) | |
49 | |
50 output_scatter_plots_by_location(location_plot_output = options.location_plot_output, vcf_info = vcf_info, loess_span=options.loess_span, d_yaxis=options.d_yaxis, h_yaxis=options.h_yaxis, points_color=options.points_color, loess_color=options.loess_color, standardize =options.standardize, normalized_hist_per_mb = normalized_histogram_bins_per_mb, normalized_hist_per_5kb = normalized_histogram_bins_per_5kb, breaks = break_dict, rounded_bin_size = rounded_bin_size, max_y_hist_overall = max_y_hist_overall) | |
51 | |
52 def myround(x, base=10): | |
53 return int(base * round(float(x)/base)) | |
54 | |
55 def skip_headers(reader = None, i_file = None): | |
56 # count headers | |
57 comment = 0 | |
58 while reader.next()[0].startswith('#'): | |
59 comment = comment + 1 | |
60 | |
61 # skip headers | |
62 i_file.seek(0) | |
63 for i in range(0, comment): | |
64 reader.next() | |
65 | |
66 def parse_breaks(break_file = None): | |
67 if break_file == 'C.elegans': | |
68 break_dict = { 'I' : 16 , 'II' : 16, 'III' : 14, 'IV' : 18, 'V' : 21, 'X' : 18 } | |
69 return break_dict | |
70 elif break_file == 'Brachypodium': | |
71 break_dict = { '1' : 75 , '2' : 60, '3' : 60, '4' : 50, '5' : 30 } | |
72 return break_dict | |
73 elif break_file == 'Arabadopsis': | |
74 break_dict = { '1' : 31 , '2' : 20, '3' : 24, '4' : 19, '5' : 27 } | |
75 return break_dict | |
76 else: | |
77 i_file = open(break_file, 'rU') | |
78 break_dict = {} | |
79 reader = csv.reader(i_file, delimiter = '\t') | |
80 for row in reader: | |
81 chromosome = row[0].upper() | |
82 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE) | |
83 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE) | |
84 break_count = row[1] | |
85 break_dict[chromosome] = int(break_count) | |
86 return break_dict | |
87 | |
88 | |
89 def location_comparer(location_1, location_2): | |
90 chr_loc_1 = location_1.split(':')[0] | |
91 pos_loc_1 = int(location_1.split(':')[1]) | |
92 | |
93 chr_loc_2 = location_2.split(':')[0] | |
94 pos_loc_2 = int(location_2.split(':')[1]) | |
95 | |
96 if chr_loc_1 == chr_loc_2: | |
97 if pos_loc_1 < pos_loc_2: | |
98 return -1 | |
99 elif pos_loc_1 == pos_loc_1: | |
100 return 0 | |
101 elif pos_loc_1 > pos_loc_2: | |
102 return 1 | |
103 elif chr_loc_1 < chr_loc_2: | |
104 return -1 | |
105 elif chr_loc_1 > chr_loc_2: | |
106 return 1 | |
107 | |
108 def output_vcf_info(output = None, vcf_info = None): | |
109 o_file = open(output, 'wb') | |
110 writer = csv.writer(o_file, delimiter = '\t') | |
111 | |
112 writer.writerow(["#Chr\t", "Pos\t", "Alt Count\t", "Ref Count\t", "Read Depth\t", "Ratio\t"]) | |
113 | |
114 location_sorted_vcf_info_keys = sorted(vcf_info.keys(), cmp=location_comparer) | |
115 | |
116 for location in location_sorted_vcf_info_keys: | |
117 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location] | |
118 | |
119 location_info = location.split(':') | |
120 chromosome = location_info[0] | |
121 position = location_info[1] | |
122 | |
123 writer.writerow([chromosome, position, alt_allele_count, ref_allele_count, read_depth, ratio]) | |
124 | |
125 o_file.close() | |
126 | |
127 def output_scatter_plots_by_location(location_plot_output = None, vcf_info = None, loess_span="", d_yaxis="", h_yaxis="", points_color="", loess_color="", standardize=None, normalized_hist_per_mb = None, normalized_hist_per_5kb = None, breaks = None, rounded_bin_size = 1000000, max_y_hist_overall = ""): | |
128 positions = {} | |
129 current_chr = "" | |
130 prev_chr = "" | |
131 | |
132 x_label = "Location (Mb)" | |
133 filtered_label = '' | |
134 | |
135 location_sorted_vcf_info_keys = sorted(vcf_info.keys(), cmp=location_comparer) | |
136 | |
137 break_unit = Decimal(rounded_bin_size) / Decimal(1000000) | |
138 max_breaks = max(breaks.values()) | |
139 | |
140 try: | |
141 r.pdf(location_plot_output, 8, 8) | |
142 | |
143 for location in location_sorted_vcf_info_keys: | |
144 current_chr = location.split(':')[0] | |
145 position = location.split(':')[1] | |
146 | |
147 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location] | |
148 | |
149 if prev_chr != current_chr: | |
150 if prev_chr != "": | |
151 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = prev_chr) | |
152 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = prev_chr) | |
153 | |
154 if h_yaxis == 0: | |
155 plot_data(chr_dict = positions, hist_dict_mb = hist_dict_mb, hist_dict_5kb = hist_dict_5kb, chr = prev_chr + filtered_label, x_label = "Location (Mb)", divide_position = True, draw_secondary_grid_lines = True, loess_span=loess_span, d_yaxis=d_yaxis, h_yaxis=max_y_hist_overall, points_color=points_color, loess_color=loess_color, breaks = breaks[prev_chr], standardize=standardize, max_breaks = max_breaks, break_unit = break_unit) | |
156 else: | |
157 plot_data(chr_dict = positions, hist_dict_mb = hist_dict_mb, hist_dict_5kb = hist_dict_5kb, chr = prev_chr + filtered_label, x_label = "Location (Mb)", divide_position = True, draw_secondary_grid_lines = True, loess_span=loess_span, d_yaxis=d_yaxis, h_yaxis=h_yaxis, points_color=points_color, loess_color=loess_color, breaks = breaks[prev_chr], standardize=standardize, max_breaks = max_breaks, break_unit = break_unit) | |
158 | |
159 prev_chr = current_chr | |
160 positions = {} | |
161 | |
162 positions[position] = ratio | |
163 | |
164 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = current_chr) | |
165 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = current_chr) | |
166 | |
167 if h_yaxis == 0: | |
168 plot_data(chr_dict = positions, hist_dict_mb = hist_dict_mb, hist_dict_5kb = hist_dict_5kb, chr = current_chr + filtered_label, x_label = "Location (Mb)", divide_position = True, draw_secondary_grid_lines = True, loess_span=loess_span, d_yaxis=d_yaxis, h_yaxis=max_y_hist_overall, points_color=points_color, loess_color=loess_color, breaks = breaks[current_chr], standardize=standardize, max_breaks = max_breaks, break_unit = break_unit) | |
169 else: | |
170 plot_data(chr_dict = positions, hist_dict_mb = hist_dict_mb, hist_dict_5kb = hist_dict_5kb, chr = current_chr + filtered_label, x_label = "Location (Mb)", divide_position = True, draw_secondary_grid_lines = True, loess_span=loess_span, d_yaxis=d_yaxis, h_yaxis=h_yaxis, points_color=points_color, loess_color=loess_color, breaks = breaks[current_chr], standardize=standardize, max_breaks = max_breaks, break_unit = break_unit) | |
171 | |
172 r.dev_off() | |
173 | |
174 except Exception as inst: | |
175 print inst | |
176 print "There was an error creating the location plot pdf... Please try again" | |
177 | |
178 def get_hist_dict_by_chr(normalized_hist_per_xbase = None, chr = ''): | |
179 hist_dict = {} | |
180 | |
181 for location in normalized_hist_per_xbase: | |
182 chromosome = location.split(':')[0] | |
183 if chromosome == chr: | |
184 position = int(location.split(':')[1]) | |
185 hist_dict[position] = normalized_hist_per_xbase[location] | |
186 | |
187 max_location = max(hist_dict.keys(), key=int) | |
188 for i in range(1, max_location): | |
189 if i not in hist_dict: | |
190 hist_dict[i] = 0 | |
191 | |
192 return hist_dict | |
193 | |
194 def plot_data(chr_dict = None, hist_dict_mb = None, hist_dict_5kb = None, chr = "", x_label = "", divide_position = False, draw_secondary_grid_lines = False, loess_span=None, d_yaxis=None, h_yaxis=None, points_color="", loess_color="", breaks = None, standardize= None, max_breaks = 1, break_unit = 1): | |
195 ratios = "c(" | |
196 positions = "c(" | |
197 | |
198 for position in chr_dict: | |
199 ratio = chr_dict[position] | |
200 if divide_position: | |
201 position = float(position) / 1000000.0 | |
202 positions = positions + str(position) + ", " | |
203 ratios = ratios + str(ratio) + ", " | |
204 | |
205 if len(ratios) == 2: | |
206 ratios = ratios + ")" | |
207 else: | |
208 ratios = ratios[0:len(ratios) - 2] + ")" | |
209 | |
210 if len(positions) == 2: | |
211 positions = positions + ")" | |
212 else: | |
213 positions = positions[0:len(positions) - 2] + ")" | |
214 | |
215 r("x <- " + positions) | |
216 r("y <- " + ratios) | |
217 | |
218 hist_mb_values = "c(" | |
219 for position in sorted(hist_dict_mb): | |
220 hist_mb_values = hist_mb_values + str(hist_dict_mb[position]) + ", " | |
221 | |
222 if len(hist_mb_values) == 2: | |
223 hist_mb_values = hist_mb_values + ")" | |
224 else: | |
225 hist_mb_values = hist_mb_values[0:len(hist_mb_values) - 2] + ")" | |
226 | |
227 hist_5kb_values = "c(" | |
228 for position in sorted(hist_dict_5kb): | |
229 hist_5kb_values = hist_5kb_values + str(hist_dict_5kb[position]) + ", " | |
230 | |
231 if len(hist_5kb_values) == 2: | |
232 hist_5kb_values = hist_5kb_values + ")" | |
233 else: | |
234 hist_5kb_values = hist_5kb_values[0:len(hist_5kb_values) - 2] + ")" | |
235 | |
236 r("xz <- " + hist_mb_values) | |
237 r("yz <- " + hist_5kb_values) | |
238 | |
239 | |
240 max_break_str = str(max_breaks) | |
241 break_unit_str = str(Decimal(break_unit)) | |
242 half_break_unit_str = str(Decimal(break_unit) / Decimal(2)) | |
243 break_penta_unit_str = str(Decimal(break_unit) * Decimal(5)) | |
244 | |
245 if (standardize=='true'): | |
246 r("plot(x, y, ,cex=0.60, xlim=c(0," + max_break_str + "), main='LG " + chr + " (Variant Discovery Mapping)', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='Ratios of variant reads/total reads (at variant positions)', pch=10, col='"+ points_color +"')") | |
247 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')") | |
248 r("axis(1, at=seq(0, " + max_break_str + ", by=" + break_unit_str + "), labels=FALSE, tcl=-0.5)") | |
249 r("axis(1, at=seq(0, " + max_break_str + ", by=" + half_break_unit_str + "), labels=FALSE, tcl=-0.25)") | |
250 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)") | |
251 elif (standardize=='false'): | |
252 r("plot(x, y, cex=0.60, main='LG " + chr + " (Variant Discovery Mapping)', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='Ratios of variant reads/total reads (at variant positions)', pch=10, col='"+ points_color +"')") | |
253 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')") | |
254 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)") | |
255 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)") | |
256 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)") | |
257 | |
258 if draw_secondary_grid_lines: | |
259 r("abline(h = seq(floor(min(y)), 1, by=0.1), v = seq(floor(min(x)), length(x), by= 1), col='gray')") | |
260 else: | |
261 r("grid(lty = 1, col = 'gray')") | |
262 | |
263 if (standardize=='true'): | |
264 r("barplot(xz, xlim=c(0, " + max_break_str + "), ylim = c(0, " + str(h_yaxis) + "), yaxp=c(0, " + str(h_yaxis) + ", 1), space = 0, col='darkgray', width = " + break_unit_str + ", xlab='Location (Mb)', ylab='Normalized frequency of pure parental alleles ', main='LG " + chr + " (Variant Discovery Mapping)')") | |
265 r("barplot(yz, space = 0, add=TRUE, width = " + half_break_unit_str + ", col='green2')") | |
266 r("axis(1, hadj = 1, at=seq(0, " + max_break_str + ", by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)") | |
267 r("axis(1, at=seq(0, " + max_break_str + ", by= " + break_penta_unit_str + "), labels=TRUE, tcl=-0.5)") | |
268 r("axis(1, at=seq(0, " + max_break_str + ", by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)") | |
269 elif (standardize=='false'): | |
270 r("barplot(xz, ylim = c(0, " + str(h_yaxis) + "), yaxp=c(0, " + str(h_yaxis) + ", 1), space = 0, col='darkgray', width = 1, xlab='Location (Mb)', ylab='Normalized frequency of pure parental alleles ', main='LG " + chr + " (Variant Discovery Mapping)')") | |
271 r("barplot(yz, space = 0, add=TRUE, width = 0.5, col='green2')") | |
272 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)") | |
273 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_penta_unit_str + "), labels=TRUE, tcl=-0.5)") | |
274 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)") | |
275 | |
276 | |
277 | |
278 def calculate_normalized_histogram_bins_per_xbase(vcf_info = None, xbase = 1000000, normalize_bins = None): | |
279 normalized_histogram_bins_per_xbase = {} | |
280 | |
281 ref_snp_count_per_xbase = get_ref_snp_count_per_xbase(vcf_info = vcf_info, xbase = xbase) | |
282 mean_one_ratio_snp_count_per_chromosome = get_mean_one_ratio_snp_count_per_chromosome(vcf_info = vcf_info, xbase = xbase) | |
283 one_ratio_snp_count_per_xbase = get_one_ratio_snp_count_per_xbase(vcf_info = vcf_info, xbase = xbase) | |
284 | |
285 for location in ref_snp_count_per_xbase: | |
286 chromosome = location.split(':')[0] | |
287 mean_one_ratio_snp_count = mean_one_ratio_snp_count_per_chromosome[chromosome] | |
288 ref_snp_count = ref_snp_count_per_xbase[location] | |
289 | |
290 one_ratio_snp_count = 0 | |
291 if location in one_ratio_snp_count_per_xbase: | |
292 one_ratio_snp_count = one_ratio_snp_count_per_xbase[location] | |
293 | |
294 if normalize_bins == 'true': | |
295 if one_ratio_snp_count == 0: | |
296 normalized_histogram_bins_per_xbase[location] = 0 | |
297 elif one_ratio_snp_count == ref_snp_count: | |
298 normalized_histogram_bins_per_xbase[location] = (Decimal(one_ratio_snp_count**2) * Decimal(mean_one_ratio_snp_count)) | |
299 else: | |
300 normalized_histogram_bins_per_xbase[location] = (Decimal(one_ratio_snp_count**2) / (Decimal(ref_snp_count)-Decimal(one_ratio_snp_count))) * Decimal(mean_one_ratio_snp_count) | |
301 else: | |
302 normalized_histogram_bins_per_xbase[location] = one_ratio_snp_count | |
303 | |
304 return normalized_histogram_bins_per_xbase | |
305 | |
306 | |
307 def get_ref_snp_count_per_xbase(vcf_info = None, xbase = 1000000): | |
308 ref_snps_per_xbase = {} | |
309 | |
310 for location in vcf_info: | |
311 location_info = location.split(':') | |
312 | |
313 chromosome = location_info[0].upper() | |
314 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE) | |
315 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE) | |
316 | |
317 position = location_info[1] | |
318 xbase_position = (int(position) / xbase) + 1 | |
319 | |
320 location = chromosome + ":" + str(xbase_position) | |
321 if location in ref_snps_per_xbase: | |
322 ref_snps_per_xbase[location] = ref_snps_per_xbase[location] + 1 | |
323 else: | |
324 ref_snps_per_xbase[location] = 1 | |
325 | |
326 return ref_snps_per_xbase | |
327 | |
328 | |
329 | |
330 def get_mean_one_ratio_snp_count_per_chromosome(vcf_info, xbase = 1000000): | |
331 sample_snp_count_per_xbase = {} | |
332 | |
333 for location in vcf_info: | |
334 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location] | |
335 | |
336 location_info = location.split(':') | |
337 chromosome = location_info[0] | |
338 position = location_info[1] | |
339 xbase_position = (int(position) / xbase) + 1 | |
340 xbase_location = chromosome + ":" + str(xbase_position) | |
341 | |
342 if ratio == 1: | |
343 if xbase_location in sample_snp_count_per_xbase: | |
344 sample_snp_count_per_xbase[xbase_location] = sample_snp_count_per_xbase[xbase_location] + 1 | |
345 else: | |
346 sample_snp_count_per_xbase[xbase_location] = 1 | |
347 | |
348 elif ratio != 1 and xbase_location not in sample_snp_count_per_xbase: | |
349 sample_snp_count_per_xbase[xbase_location] = 0 | |
350 | |
351 mean_one_ratio_snp_count_per_chromosome = {} | |
352 for location in sample_snp_count_per_xbase: | |
353 chromosome = location.split(':')[0] | |
354 sample_count = sample_snp_count_per_xbase[location] | |
355 if chromosome in mean_one_ratio_snp_count_per_chromosome: | |
356 mean_one_ratio_snp_count_per_chromosome[chromosome].append(sample_count) | |
357 else: | |
358 mean_one_ratio_snp_count_per_chromosome[chromosome] = [sample_count] | |
359 | |
360 for chromosome in mean_one_ratio_snp_count_per_chromosome: | |
361 summa = sum(mean_one_ratio_snp_count_per_chromosome[chromosome]) | |
362 count = len(mean_one_ratio_snp_count_per_chromosome[chromosome]) | |
363 | |
364 mean_one_ratio_snp_count_per_chromosome[chromosome] = Decimal(summa) / Decimal(count) | |
365 | |
366 return mean_one_ratio_snp_count_per_chromosome | |
367 | |
368 | |
369 def get_one_ratio_snp_count_per_xbase(vcf_info = None, xbase = 1000000): | |
370 one_ratio_snp_count_per_xbase = {} | |
371 | |
372 for location in vcf_info: | |
373 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location] | |
374 | |
375 location_info = location.split(':') | |
376 chromosome = location_info[0] | |
377 position = location_info[1] | |
378 xbase_position = (int(position) / xbase) + 1 | |
379 xbase_location = chromosome + ":" + str(xbase_position) | |
380 | |
381 if ratio == 1: | |
382 if xbase_location in one_ratio_snp_count_per_xbase: | |
383 one_ratio_snp_count_per_xbase[xbase_location] = one_ratio_snp_count_per_xbase[xbase_location] + 1 | |
384 else: | |
385 one_ratio_snp_count_per_xbase[xbase_location] = 1 | |
386 | |
387 elif ratio != 1 and xbase_location not in one_ratio_snp_count_per_xbase: | |
388 one_ratio_snp_count_per_xbase[xbase_location] = 0 | |
389 | |
390 return one_ratio_snp_count_per_xbase | |
391 | |
392 | |
393 def parse_vcf(sample_vcf = None): | |
394 i_file = open(sample_vcf, 'rU') | |
395 reader = csv.reader(i_file, delimiter = '\t', quoting = csv.QUOTE_NONE) | |
396 | |
397 skip_headers(reader = reader, i_file = i_file) | |
398 vcf_info = {} | |
399 | |
400 for row in reader: | |
401 chromosome = row[0].upper() | |
402 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE) | |
403 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE) | |
404 | |
405 | |
406 if chromosome != 'MTDNA': | |
407 position = row[1] | |
408 #ref_allele = row[2] | |
409 #read_depth = row[3] | |
410 #read_bases = row[4] | |
411 | |
412 vcf_format_info = row[8].split(":") | |
413 vcf_allele_freq_data = row[9] | |
414 | |
415 read_depth_data_index = vcf_format_info.index("DP") | |
416 read_depth = vcf_allele_freq_data.split(":")[read_depth_data_index] | |
417 | |
418 ref_and_alt_counts_data_index = vcf_format_info.index("AD") | |
419 ref_and_alt_counts = vcf_allele_freq_data.split(":")[ref_and_alt_counts_data_index] | |
420 ref_allele_count = ref_and_alt_counts.split(",")[0] | |
421 alt_allele_count = ref_and_alt_counts.split(",")[1] | |
422 | |
423 location = chromosome + ":" + position | |
424 | |
425 if Decimal(read_depth!=0): | |
426 getcontext().prec = 6 | |
427 ratio = Decimal(alt_allele_count) / Decimal(read_depth) | |
428 | |
429 vcf_info[location] = (alt_allele_count, ref_allele_count, read_depth, ratio) | |
430 | |
431 #debug line | |
432 #print chromosome, position, read_depth, ref_allele_count, alt_allele_count, ratio, id | |
433 | |
434 i_file.close() | |
435 | |
436 return vcf_info | |
437 | |
438 def parse_read_bases(read_bases = None, alt_allele = None): | |
439 read_bases = re.sub('\$', '', read_bases) | |
440 read_bases = re.sub('\^[^\s]', '', read_bases) | |
441 | |
442 ref_allele_matches = re.findall("\.|\,", read_bases) | |
443 ref_allele_count = len(ref_allele_matches) | |
444 | |
445 alt_allele_matches = re.findall(alt_allele, read_bases, flags = re.IGNORECASE) | |
446 alt_allele_count = len(alt_allele_matches) | |
447 | |
448 #debug line | |
449 #print read_bases, alt_allele, alt_allele_count, ref_allele_count | |
450 | |
451 return ref_allele_count, alt_allele_count | |
452 | |
453 if __name__ == "__main__": | |
454 main() |