|
38
|
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 = "red", 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
|
|
|
28 parser.add_option('-o', '--output', dest = 'output', action = 'store', type = 'string', default = None, help = "Output file name")
|
|
|
29 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")
|
|
|
30
|
|
|
31 (options, args) = parser.parse_args()
|
|
|
32
|
|
|
33 vcf_info = parse_vcf(sample_vcf = options.sample_vcf)
|
|
|
34
|
|
|
35 output_vcf_info(output = options.output, vcf_info = vcf_info)
|
|
|
36
|
|
|
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 break_dict = parse_breaks(break_file = options.break_file)
|
|
|
48
|
|
|
49 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)
|
|
|
50
|
|
|
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 == 'Arabidopsis':
|
|
|
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 #Brachy
|
|
|
85 chromosome = re.sub("Bd", "", chromosome, flags = re.IGNORECASE)
|
|
|
86 chromosome = re.sub("bd", "", chromosome, flags = re.IGNORECASE)
|
|
|
87
|
|
|
88 break_count = row[1]
|
|
|
89 break_dict[chromosome] = int(break_count)
|
|
|
90 return break_dict
|
|
|
91
|
|
|
92
|
|
|
93 def location_comparer(location_1, location_2):
|
|
|
94 chr_loc_1 = location_1.split(':')[0]
|
|
|
95 pos_loc_1 = int(location_1.split(':')[1])
|
|
|
96
|
|
|
97 chr_loc_2 = location_2.split(':')[0]
|
|
|
98 pos_loc_2 = int(location_2.split(':')[1])
|
|
|
99
|
|
|
100 if chr_loc_1 == chr_loc_2:
|
|
|
101 if pos_loc_1 < pos_loc_2:
|
|
|
102 return -1
|
|
|
103 elif pos_loc_1 == pos_loc_1:
|
|
|
104 return 0
|
|
|
105 elif pos_loc_1 > pos_loc_2:
|
|
|
106 return 1
|
|
|
107 elif chr_loc_1 < chr_loc_2:
|
|
|
108 return -1
|
|
|
109 elif chr_loc_1 > chr_loc_2:
|
|
|
110 return 1
|
|
|
111
|
|
|
112 def output_vcf_info(output = None, vcf_info = None):
|
|
|
113 o_file = open(output, 'wb')
|
|
|
114 writer = csv.writer(o_file, delimiter = '\t')
|
|
|
115
|
|
|
116 writer.writerow(["#Chr\t", "Pos\t", "Alt Count\t", "Ref Count\t", "Read Depth\t", "Ratio\t"])
|
|
|
117
|
|
|
118 location_sorted_vcf_info_keys = sorted(vcf_info.keys(), cmp=location_comparer)
|
|
|
119
|
|
|
120 for location in location_sorted_vcf_info_keys:
|
|
|
121 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location]
|
|
|
122
|
|
|
123 location_info = location.split(':')
|
|
|
124 chromosome = location_info[0]
|
|
|
125 position = location_info[1]
|
|
|
126
|
|
|
127 writer.writerow([chromosome, position, alt_allele_count, ref_allele_count, read_depth, ratio])
|
|
|
128
|
|
|
129 o_file.close()
|
|
|
130
|
|
|
131 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 = ""):
|
|
|
132 positions = {}
|
|
|
133 current_chr = ""
|
|
|
134 prev_chr = ""
|
|
|
135
|
|
|
136 x_label = "Location (Mb)"
|
|
|
137 filtered_label = ''
|
|
|
138
|
|
|
139 location_sorted_vcf_info_keys = sorted(vcf_info.keys(), cmp=location_comparer)
|
|
|
140
|
|
|
141 break_unit = Decimal(rounded_bin_size) / Decimal(1000000)
|
|
|
142 max_breaks = max(breaks.values())
|
|
|
143
|
|
|
144 try:
|
|
|
145 r.pdf(location_plot_output, 8, 8)
|
|
|
146
|
|
|
147 for location in location_sorted_vcf_info_keys:
|
|
|
148 current_chr = location.split(':')[0]
|
|
|
149 position = location.split(':')[1]
|
|
|
150
|
|
|
151 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location]
|
|
|
152
|
|
|
153 if prev_chr != current_chr:
|
|
|
154 if prev_chr != "":
|
|
|
155 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = prev_chr)
|
|
|
156 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = prev_chr)
|
|
|
157
|
|
|
158 if h_yaxis == 0:
|
|
|
159 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)
|
|
|
160 else:
|
|
|
161 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)
|
|
|
162
|
|
|
163 prev_chr = current_chr
|
|
|
164 positions = {}
|
|
|
165
|
|
|
166 positions[position] = ratio
|
|
|
167
|
|
|
168 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = current_chr)
|
|
|
169 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = current_chr)
|
|
|
170
|
|
|
171 if h_yaxis == 0:
|
|
|
172 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)
|
|
|
173 else:
|
|
|
174 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)
|
|
|
175
|
|
|
176 r.dev_off()
|
|
|
177
|
|
|
178 except Exception as inst:
|
|
|
179 print inst
|
|
|
180 print "There was an error creating the location plot pdf... Please try again"
|
|
|
181
|
|
|
182 def get_hist_dict_by_chr(normalized_hist_per_xbase = None, chr = ''):
|
|
|
183 hist_dict = {}
|
|
|
184
|
|
|
185 for location in normalized_hist_per_xbase:
|
|
|
186 chromosome = location.split(':')[0]
|
|
|
187 if chromosome == chr:
|
|
|
188 position = int(location.split(':')[1])
|
|
|
189 hist_dict[position] = normalized_hist_per_xbase[location]
|
|
|
190
|
|
|
191 max_location = max(hist_dict.keys(), key=int)
|
|
|
192 for i in range(1, max_location):
|
|
|
193 if i not in hist_dict:
|
|
|
194 hist_dict[i] = 0
|
|
|
195
|
|
|
196 return hist_dict
|
|
|
197
|
|
|
198
|
|
|
199 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):
|
|
|
200 ratios = "c("
|
|
|
201 positions = "c("
|
|
|
202
|
|
|
203 for position in chr_dict:
|
|
|
204 ratio = chr_dict[position]
|
|
|
205 if divide_position:
|
|
|
206 position = float(position) / 1000000.0
|
|
|
207 positions = positions + str(position) + ", "
|
|
|
208 ratios = ratios + str(ratio) + ", "
|
|
|
209
|
|
|
210 if len(ratios) == 2:
|
|
|
211 ratios = ratios + ")"
|
|
|
212 else:
|
|
|
213 ratios = ratios[0:len(ratios) - 2] + ")"
|
|
|
214
|
|
|
215 if len(positions) == 2:
|
|
|
216 positions = positions + ")"
|
|
|
217 else:
|
|
|
218 positions = positions[0:len(positions) - 2] + ")"
|
|
|
219
|
|
|
220 r("x <- " + positions)
|
|
|
221 r("y <- " + ratios)
|
|
|
222
|
|
|
223 hist_mb_values = "c("
|
|
|
224 for position in sorted(hist_dict_mb):
|
|
|
225 hist_mb_values = hist_mb_values + str(hist_dict_mb[position]) + ", "
|
|
|
226
|
|
|
227 if len(hist_mb_values) == 2:
|
|
|
228 hist_mb_values = hist_mb_values + ")"
|
|
|
229 else:
|
|
|
230 hist_mb_values = hist_mb_values[0:len(hist_mb_values) - 2] + ")"
|
|
|
231
|
|
|
232 hist_5kb_values = "c("
|
|
|
233 for position in sorted(hist_dict_5kb):
|
|
|
234 hist_5kb_values = hist_5kb_values + str(hist_dict_5kb[position]) + ", "
|
|
|
235
|
|
|
236 if len(hist_5kb_values) == 2:
|
|
|
237 hist_5kb_values = hist_5kb_values + ")"
|
|
|
238 else:
|
|
|
239 hist_5kb_values = hist_5kb_values[0:len(hist_5kb_values) - 2] + ")"
|
|
|
240
|
|
|
241 r("xz <- " + hist_mb_values)
|
|
|
242 r("yz <- " + hist_5kb_values)
|
|
|
243
|
|
|
244
|
|
|
245 max_break_str = str(max_breaks)
|
|
|
246 break_unit_str = str(Decimal(break_unit))
|
|
|
247 half_break_unit_str = str(Decimal(break_unit) / Decimal(2))
|
|
|
248 break_penta_unit_str = str(Decimal(break_unit) * Decimal(5))
|
|
|
249
|
|
|
250 if (standardize=='true'):
|
|
|
251 r("plot(x, y, cex=0.60, xlim=c(0," + max_break_str + "), main='LG " + chr + "', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='Percentage of mapping strain alleles/total reads (at SNP positions)', pch=10, col='"+ points_color +"')")
|
|
|
252 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
|
|
|
253 r("axis(1, at=seq(0, " + max_break_str + ", by=" + break_unit_str + "), labels=FALSE, tcl=-0.5)")
|
|
|
254 r("axis(1, at=seq(0, " + max_break_str + ", by=" + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
|
|
|
255 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
|
|
|
256 elif (standardize=='false'):
|
|
|
257 r("plot(x, y, cex=0.60, main='LG " + chr + " ', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='Percentage of mapping strain alleles/total reads (at SNP positions)', pch=10, col='"+ points_color +"')")
|
|
|
258 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
|
|
|
259 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
|
|
|
260 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
|
|
|
261 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
|
|
|
262
|
|
|
263 if draw_secondary_grid_lines:
|
|
|
264 r("abline(h = seq(floor(min(y)), 1, by=0.1), v = seq(floor(min(x)), length(x), by= 1), col='gray')")
|
|
|
265 else:
|
|
|
266 r("grid(lty = 1, col = 'gray')")
|
|
|
267
|
|
|
268 if (standardize=='true'):
|
|
|
269 #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 + " (Hawaiian Variant Mapping)')")
|
|
|
270 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 + "')")
|
|
|
271 r("barplot(yz, space = 0, add=TRUE, width = " + half_break_unit_str + ", col=rgb(1, 0, 0, 1))")
|
|
|
272 r("axis(1, hadj = 1, at=seq(0, " + max_break_str + ", by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
|
|
|
273 r("axis(1, at=seq(0, " + max_break_str + ", by= " + break_penta_unit_str + "), labels=TRUE, tcl=-0.5)")
|
|
|
274 r("axis(1, at=seq(0, " + max_break_str + ", by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
|
|
|
275 elif (standardize=='false'):
|
|
|
276 #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 + " (Hawaiian Variant Mapping)')")
|
|
|
277 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 + " ')")
|
|
|
278 r("barplot(yz, space = 0, add=TRUE, width = 0.5, col=rgb(1, 0, 0, 1))")
|
|
|
279 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
|
|
|
280 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_penta_unit_str + "), labels=TRUE, tcl=-0.5)")
|
|
|
281 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
|
|
|
282
|
|
|
283
|
|
|
284
|
|
|
285 def calculate_normalized_histogram_bins_per_xbase(vcf_info = None, xbase = 1000000, normalize_bins = None):
|
|
|
286 normalized_histogram_bins_per_xbase = {}
|
|
|
287
|
|
|
288 ref_snp_count_per_xbase = get_ref_snp_count_per_xbase(vcf_info = vcf_info, xbase = xbase)
|
|
|
289
|
|
|
290 mean_zero_snp_count_per_chromosome = get_mean_zero_snp_count_per_chromosome(vcf_info = vcf_info, xbase = xbase)
|
|
|
291
|
|
|
292 zero_snp_count_per_xbase = get_zero_snp_count_per_xbase(vcf_info = vcf_info, xbase = xbase)
|
|
|
293
|
|
|
294
|
|
|
295 for location in ref_snp_count_per_xbase:
|
|
|
296 chromosome = location.split(':')[0]
|
|
|
297 mean_zero_snp_count = mean_zero_snp_count_per_chromosome[chromosome]
|
|
|
298 ref_snp_count = ref_snp_count_per_xbase[location]
|
|
|
299
|
|
|
300 zero_snp_count = 0
|
|
|
301 if location in zero_snp_count_per_xbase:
|
|
|
302 zero_snp_count = zero_snp_count_per_xbase[location]
|
|
|
303
|
|
|
304 if normalize_bins == 'true':
|
|
|
305 if zero_snp_count == 0 or ref_snp_count == 0:
|
|
|
306 normalized_histogram_bins_per_xbase[location] = 0
|
|
|
307 elif zero_snp_count == ref_snp_count:
|
|
|
308 normalized_histogram_bins_per_xbase[location] = 0
|
|
|
309 else:
|
|
|
310 normalized_histogram_bins_per_xbase[location] = (Decimal(zero_snp_count) / (Decimal(ref_snp_count)-Decimal(zero_snp_count))) * Decimal(mean_zero_snp_count)
|
|
|
311 else:
|
|
|
312 normalized_histogram_bins_per_xbase[location] = zero_snp_count
|
|
|
313
|
|
|
314 return normalized_histogram_bins_per_xbase
|
|
|
315
|
|
|
316
|
|
|
317 def get_ref_snp_count_per_xbase(vcf_info = None, xbase = 1000000):
|
|
|
318 ref_snps_per_xbase = {}
|
|
|
319
|
|
|
320 for location in vcf_info:
|
|
|
321 location_info = location.split(':')
|
|
|
322
|
|
|
323 chromosome = location_info[0].upper()
|
|
|
324 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
|
|
|
325 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
|
|
|
326
|
|
|
327 #Brachy
|
|
|
328 chromosome = re.sub("Bd", "", chromosome, flags = re.IGNORECASE)
|
|
|
329 chromosome = re.sub("bd", "", chromosome, flags = re.IGNORECASE)
|
|
|
330
|
|
|
331 position = location_info[1]
|
|
|
332 xbase_position = (int(position) / xbase) + 1
|
|
|
333
|
|
|
334 location = chromosome + ":" + str(xbase_position)
|
|
|
335 if location in ref_snps_per_xbase:
|
|
|
336 ref_snps_per_xbase[location] = ref_snps_per_xbase[location] + 1
|
|
|
337 else:
|
|
|
338 ref_snps_per_xbase[location] = 1
|
|
|
339
|
|
|
340 return ref_snps_per_xbase
|
|
|
341
|
|
|
342
|
|
|
343
|
|
|
344 def get_mean_zero_snp_count_per_chromosome(vcf_info, xbase = 1000000):
|
|
|
345 sample_snp_count_per_xbase = {}
|
|
|
346
|
|
|
347 for location in vcf_info:
|
|
|
348 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location]
|
|
|
349
|
|
|
350 location_info = location.split(':')
|
|
|
351 chromosome = location_info[0]
|
|
|
352 position = location_info[1]
|
|
|
353 xbase_position = (int(position) / xbase) + 1
|
|
|
354 xbase_location = chromosome + ":" + str(xbase_position)
|
|
|
355
|
|
|
356 if int(alt_allele_count) == 0:
|
|
|
357 if xbase_location in sample_snp_count_per_xbase:
|
|
|
358 sample_snp_count_per_xbase[xbase_location] = sample_snp_count_per_xbase[xbase_location] + 1
|
|
|
359 else:
|
|
|
360 sample_snp_count_per_xbase[xbase_location] = 1
|
|
|
361
|
|
|
362 elif int(alt_allele_count) != 0 and xbase_location not in sample_snp_count_per_xbase:
|
|
|
363 sample_snp_count_per_xbase[xbase_location] = 0
|
|
|
364
|
|
|
365 mean_zero_snp_count_per_chromosome = {}
|
|
|
366 for location in sample_snp_count_per_xbase:
|
|
|
367 chromosome = location.split(':')[0]
|
|
|
368 sample_count = sample_snp_count_per_xbase[location]
|
|
|
369 if chromosome in mean_zero_snp_count_per_chromosome:
|
|
|
370 mean_zero_snp_count_per_chromosome[chromosome].append(sample_count)
|
|
|
371 else:
|
|
|
372 mean_zero_snp_count_per_chromosome[chromosome] = [sample_count]
|
|
|
373
|
|
|
374 for chromosome in mean_zero_snp_count_per_chromosome:
|
|
|
375 summa = sum(mean_zero_snp_count_per_chromosome[chromosome])
|
|
|
376 count = len(mean_zero_snp_count_per_chromosome[chromosome])
|
|
|
377
|
|
|
378 mean_zero_snp_count_per_chromosome[chromosome] = Decimal(summa) / Decimal(count)
|
|
|
379
|
|
|
380 return mean_zero_snp_count_per_chromosome
|
|
|
381
|
|
|
382
|
|
|
383 def get_zero_snp_count_per_xbase(vcf_info = None, xbase = 1000000):
|
|
|
384 zero_snp_count_per_xbase = {}
|
|
|
385
|
|
|
386 for location in vcf_info:
|
|
|
387 alt_allele_count, ref_allele_count, read_depth, ratio = vcf_info[location]
|
|
|
388
|
|
|
389 location_info = location.split(':')
|
|
|
390 chromosome = location_info[0]
|
|
|
391 position = location_info[1]
|
|
|
392 xbase_position = (int(position) / xbase) + 1
|
|
|
393 xbase_location = chromosome + ":" + str(xbase_position)
|
|
|
394
|
|
|
395 if int(alt_allele_count) == 0:
|
|
|
396 if xbase_location in zero_snp_count_per_xbase:
|
|
|
397 zero_snp_count_per_xbase[xbase_location] = zero_snp_count_per_xbase[xbase_location] + 1
|
|
|
398 else:
|
|
|
399 zero_snp_count_per_xbase[xbase_location] = 1
|
|
|
400
|
|
|
401 elif int(alt_allele_count) != 0 and xbase_location not in zero_snp_count_per_xbase:
|
|
|
402 zero_snp_count_per_xbase[xbase_location] = 0
|
|
|
403
|
|
|
404 return zero_snp_count_per_xbase
|
|
|
405
|
|
|
406
|
|
|
407 def parse_vcf(sample_vcf = None):
|
|
|
408 i_file = open(sample_vcf, 'rU')
|
|
|
409 reader = csv.reader(i_file, delimiter = '\t', quoting = csv.QUOTE_NONE)
|
|
|
410
|
|
|
411 skip_headers(reader = reader, i_file = i_file)
|
|
|
412 vcf_info = {}
|
|
|
413
|
|
|
414 for row in reader:
|
|
|
415 chromosome = row[0].upper()
|
|
|
416 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
|
|
|
417 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
|
|
|
418
|
|
|
419 #Brachy
|
|
|
420 chromosome = re.sub("Bd", "", chromosome, flags = re.IGNORECASE)
|
|
|
421 chromosome = re.sub("bd_", "", chromosome, flags = re.IGNORECASE)
|
|
|
422
|
|
|
423 if chromosome != 'MTDNA':
|
|
|
424 position = row[1]
|
|
|
425 #ref_allele = row[2]
|
|
|
426 #read_depth = row[3]
|
|
|
427 #read_bases = row[4]
|
|
|
428
|
|
|
429 vcf_format_info = row[8].split(":")
|
|
|
430 vcf_allele_freq_data = row[9]
|
|
|
431
|
|
|
432 read_depth_data_index = vcf_format_info.index("DP")
|
|
|
433 read_depth = vcf_allele_freq_data.split(":")[read_depth_data_index]
|
|
|
434
|
|
|
435 ref_and_alt_counts_data_index = vcf_format_info.index("AD")
|
|
|
436 ref_and_alt_counts = vcf_allele_freq_data.split(":")[ref_and_alt_counts_data_index]
|
|
|
437 ref_allele_count = ref_and_alt_counts.split(",")[0]
|
|
|
438 alt_allele_count = ref_and_alt_counts.split(",")[1]
|
|
|
439
|
|
|
440 location = chromosome + ":" + position
|
|
|
441
|
|
|
442 if (Decimal(read_depth)!=0):
|
|
|
443 getcontext().prec = 6
|
|
|
444 ratio = Decimal(alt_allele_count) / Decimal(read_depth)
|
|
|
445
|
|
|
446 vcf_info[location] = (alt_allele_count, ref_allele_count, read_depth, ratio)
|
|
|
447
|
|
|
448 #debug line
|
|
|
449 #print chromosome, position, read_depth, ref_allele_count, alt_allele_count, ratio, id
|
|
|
450
|
|
|
451 i_file.close()
|
|
|
452
|
|
|
453 return vcf_info
|
|
|
454
|
|
|
455 def parse_read_bases(read_bases = None, alt_allele = None):
|
|
|
456 read_bases = re.sub('\$', '', read_bases)
|
|
|
457 read_bases = re.sub('\^[^\s]', '', read_bases)
|
|
|
458
|
|
|
459 ref_allele_matches = re.findall("\.|\,", read_bases)
|
|
|
460 ref_allele_count = len(ref_allele_matches)
|
|
|
461
|
|
|
462 alt_allele_matches = re.findall(alt_allele, read_bases, flags = re.IGNORECASE)
|
|
|
463 alt_allele_count = len(alt_allele_matches)
|
|
|
464
|
|
|
465 #debug line
|
|
|
466 #print read_bases, alt_allele, alt_allele_count, ref_allele_count
|
|
|
467
|
|
|
468 return ref_allele_count, alt_allele_count
|
|
|
469
|
|
|
470 if __name__ == "__main__":
|
|
|
471 main()
|
