Mercurial > repos > gregory-minevich > snp_mapping_using_wgs

comparison SNP_Mapping.py @ 12:9c28b8aebe84 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author gregory-minevich
date Thu, 14 Jun 2012 20:34:56 -0400
parents
children
comparison
equal deleted inserted replaced
11:cb4c388fb155 12:9c28b8aebe84
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('-p', '--sample_pileup', dest = 'sample_pileup', action = 'store', type = 'string', default = None, help = "Sample pileup from mpileup")
17 parser.add_option('-v', '--haw_vcf', dest = 'haw_vcf', action = 'store', type = 'string', default = None, help = "vcf file of Hawaiian SNPs")
18 parser.add_option('-l', '--loess_span', dest = 'loess_span', action = 'store', type = 'float', default = .01, help = "Loess span")
19 parser.add_option('-d', '--d_yaxis', dest = 'd_yaxis', action = 'store', type = 'float', default = .7, help = "y-axis upper limit for dot plot")
20 parser.add_option('-y', '--h_yaxis', dest = 'h_yaxis', action = 'store', type = 'int', default = 5, help = "y-axis upper limit for histogram plot")
21 parser.add_option('-c', '--points_color', dest = 'points_color', action = 'store', type = 'string', default = "gray27", help = "Color for data points")
22 parser.add_option('-k', '--loess_color', dest = 'loess_color', action = 'store', type = 'string', default = "red", help = "Color for loess regression line")
23 parser.add_option('-z', '--standardize', dest = 'standardize', default= 'true', help = "Standardize X-axis")
24 parser.add_option('-b', '--break_file', dest = 'break_file', action = 'store', type = 'string', default = 'C.elegans', help = "File defining the breaks per chromosome")
25 parser.add_option('-x', '--bin_size', dest = 'bin_size', action = 'store', type = 'int', default = 1000000, help = "Size of histogram bins, default is 1mb")
26 parser.add_option('-n', '--do_not_normalize_bin', dest = 'do_not_normalize_bin', action = 'store_true', help = "Do not Normalize histograms")
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 #For plotting with map units on the X-axis instead of physical distance
32 #parser.add_option('-u', '--mpu_plot_output', dest = 'mpu_plot_output', action = 'store', type = 'string', default = None, help = "Output file name of SNP plots by map unit location")
33 (options, args) = parser.parse_args()
34
35 haw_snps = build_haw_snp_dictionary(haw_vcf = options.haw_vcf)
36 pileup_info = parse_pileup(sample_pileup = options.sample_pileup, haw_snps = haw_snps)
37
38 output_pileup_info(output = options.output, pileup_info = pileup_info)
39
40 #output plot with all ratios
41 rounded_bin_size = int(round((float(options.bin_size) / 1000000), 1) * 1000000)
42
43 normalized_histogram_bins_per_mb = calculate_normalized_histogram_bins_per_xbase(pileup_info = pileup_info, xbase = rounded_bin_size, do_not_normalize = options.do_not_normalize_bin)
44 normalized_histogram_bins_per_5kb = calculate_normalized_histogram_bins_per_xbase(pileup_info = pileup_info, xbase = (rounded_bin_size / 2), do_not_normalize = options.do_not_normalize_bin)
45
46 break_dict = parse_breaks(break_file = options.break_file)
47
48 output_scatter_plots_by_location(location_plot_output = options.location_plot_output, pileup_info = pileup_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)
49
50 #For plotting with map units on the X-axis instead of physical distance)
51 #output_scatter_plots_by_mapping_units(mpu_plot_output = options.mpu_plot_output, pileup_info = pileup_info)
52
53 def skip_headers(reader = None, i_file = None):
54 # count headers
55 comment = 0
56 while reader.next()[0].startswith('#'):
57 comment = comment + 1
58
59 # skip headers
60 i_file.seek(0)
61 for i in range(0, comment):
62 reader.next()
63
64 def parse_breaks(break_file = None):
65 if break_file == 'C.elegans':
66 break_dict = { 'I' : 16 , 'II' : 16, 'III' : 14, 'IV' : 18, 'V' : 21, 'X' : 18 }
67 return break_dict
68 elif break_file == 'Arabadopsis':
69 break_dict = { '1' : 16 , '2' : 16, '3' : 21, '4' : 18, '5' : 21 }
70 return break_dict
71 else:
72 i_file = open(break_file, 'rU')
73 break_dict = {}
74 reader = csv.reader(i_file, delimiter = '\t')
75 for row in reader:
76 chromosome = row[0].upper()
77 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
78 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
79 break_count = row[1]
80 break_dict[chromosome] = int(break_count)
81 return break_dict
82
83
84 def location_comparer(location_1, location_2):
85 chr_loc_1 = location_1.split(':')[0]
86 pos_loc_1 = int(location_1.split(':')[1])
87
88 chr_loc_2 = location_2.split(':')[0]
89 pos_loc_2 = int(location_2.split(':')[1])
90
91 if chr_loc_1 == chr_loc_2:
92 if pos_loc_1 < pos_loc_2:
93 return -1
94 elif pos_loc_1 == pos_loc_1:
95 return 0
96 elif pos_loc_1 > pos_loc_2:
97 return 1
98 elif chr_loc_1 < chr_loc_2:
99 return -1
100 elif chr_loc_1 > chr_loc_2:
101 return 1
102
103 def output_pileup_info(output = None, pileup_info = None):
104 o_file = open(output, 'wb')
105 writer = csv.writer(o_file, delimiter = '\t')
106
107 writer.writerow(["#Chr\t", "Pos\t", "ID\t", "Alt Count\t", "Ref Count\t", "Read Depth\t", "Ratio\t", "Mapping Unit"])
108
109 location_sorted_pileup_info_keys = sorted(pileup_info.keys(), cmp=location_comparer)
110
111 for location in location_sorted_pileup_info_keys:
112 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
113
114 location_info = location.split(':')
115 chromosome = location_info[0]
116 position = location_info[1]
117
118 writer.writerow([chromosome, position, haw_snp_id, alt_allele_count, ref_allele_count, read_depth, ratio, mapping_unit])
119
120 o_file.close()
121
122 def output_scatter_plots_by_location(location_plot_output = None, pileup_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):
123 positions = {}
124 current_chr = ""
125 prev_chr = ""
126
127 x_label = "Location (Mb)"
128 filtered_label = ''
129
130 location_sorted_pileup_info_keys = sorted(pileup_info.keys(), cmp=location_comparer)
131
132 break_unit = Decimal(rounded_bin_size) / Decimal(1000000)
133 max_breaks = max(breaks.values())
134
135 try:
136 r.pdf(location_plot_output, 8, 8)
137
138 for location in location_sorted_pileup_info_keys:
139 current_chr = location.split(':')[0]
140 position = location.split(':')[1]
141
142 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
143
144 if prev_chr != current_chr:
145 if prev_chr != "":
146 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = prev_chr)
147 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = prev_chr)
148
149 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)
150
151 prev_chr = current_chr
152 positions = {}
153
154 positions[position] = ratio
155
156 hist_dict_mb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_mb, chr = current_chr)
157 hist_dict_5kb = get_hist_dict_by_chr(normalized_hist_per_xbase = normalized_hist_per_5kb, chr = current_chr)
158
159 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)
160
161 r.dev_off()
162
163 except Exception as inst:
164 print inst
165 print "There was an error creating the location plot pdf... Please try again"
166
167 def get_hist_dict_by_chr(normalized_hist_per_xbase = None, chr = ''):
168 hist_dict = {}
169
170 for location in normalized_hist_per_xbase:
171 chromosome = location.split(':')[0]
172 if chromosome == chr:
173 position = int(location.split(':')[1])
174 hist_dict[position] = normalized_hist_per_xbase[location]
175
176 return hist_dict
177
178 '''
179 def output_scatter_plots_by_mapping_units(mpu_plot_output = None, pileup_info = None):
180 i = {}
181 ii = {}
182 iii = {}
183 iv = {}
184 v = {}
185 x = {}
186
187 for location in pileup_info:
188 chromosome = location.split(':')[0]
189 position = location.split(':')[1]
190
191 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
192
193 if chromosome == "I":
194 i[mapping_unit] = ratio
195 elif chromosome == "II":
196 ii[mapping_unit] = ratio
197 elif chromosome == "III":
198 iii[mapping_unit] = ratio
199 elif chromosome == "IV":
200 iv[mapping_unit] = ratio
201 elif chromosome == "V":
202 v[mapping_unit] = ratio
203 elif chromosome == "X":
204 x[mapping_unit] = ratio
205
206 x_label = "Map Units"
207
208 try:
209 r.pdf(mpu_plot_output, 8, 8)
210 plot_data(chr_dict = i, chr = "I", x_label = "Map Units")
211 plot_data(chr_dict = ii, chr = "II", x_label = "Map Units")
212 plot_data(chr_dict = iii, chr = "III", x_label = "Map Units")
213 plot_data(chr_dict = iv, chr = "IV", x_label = "Map Units")
214 plot_data(chr_dict = v, chr = "V", x_label = "Map Units")
215 plot_data(chr_dict = x, chr = "X", x_label = "Map Units")
216 r.dev_off()
217 except Exception as inst:
218 print inst
219 print "There was an error creating the map unit plot pdf... Please try again"
220 '''
221
222 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):
223 ratios = "c("
224 positions = "c("
225
226 for position in chr_dict:
227 ratio = chr_dict[position]
228 if divide_position:
229 position = float(position) / 1000000.0
230 positions = positions + str(position) + ", "
231 ratios = ratios + str(ratio) + ", "
232
233 if len(ratios) == 2:
234 ratios = ratios + ")"
235 else:
236 ratios = ratios[0:len(ratios) - 2] + ")"
237
238 if len(positions) == 2:
239 positions = positions + ")"
240 else:
241 positions = positions[0:len(positions) - 2] + ")"
242
243 r("x <- " + positions)
244 r("y <- " + ratios)
245
246 hist_mb_values = "c("
247 for position in sorted(hist_dict_mb):
248 hist_mb_values = hist_mb_values + str(hist_dict_mb[position]) + ", "
249
250 if len(hist_mb_values) == 2:
251 hist_mb_values = hist_mb_values + ")"
252 else:
253 hist_mb_values = hist_mb_values[0:len(hist_mb_values) - 2] + ")"
254
255 hist_5kb_values = "c("
256 for position in sorted(hist_dict_5kb):
257 hist_5kb_values = hist_5kb_values + str(hist_dict_5kb[position]) + ", "
258
259 if len(hist_5kb_values) == 2:
260 hist_5kb_values = hist_5kb_values + ")"
261 else:
262 hist_5kb_values = hist_5kb_values[0:len(hist_5kb_values) - 2] + ")"
263
264 r("xz <- " + hist_mb_values)
265 r("yz <- " + hist_5kb_values)
266
267 max_break_str = str(max_breaks)
268 break_unit_str = str(Decimal(break_unit))
269 half_break_unit_str = str(Decimal(break_unit) / Decimal(2))
270 break_penta_unit_str = str(Decimal(break_unit) * Decimal(5))
271
272 if (standardize=='true'):
273 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='Ratios of mapping strain alleles/total reads (at SNP positions)', pch=18, col='"+ points_color +"')")
274 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
275 r("axis(1, at=seq(0, " + max_break_str + ", by=" + break_unit_str + "), labels=FALSE, tcl=-0.5)")
276 r("axis(1, at=seq(0, " + max_break_str + ", by=" + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
277 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
278 elif (standardize=='false'):
279 r("plot(x, y, cex=0.60, main='LG " + chr + "', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='Ratios of mapping strain alleles/total reads (at SNP positions)', pch=18, col='"+ points_color +"')")
280 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
281 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
282 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
283 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
284
285 if draw_secondary_grid_lines:
286 r("abline(h = seq(floor(min(y)), 1, by=0.1), v = seq(floor(min(x)), length(x), by= 1), col='gray')")
287 else:
288 r("grid(lty = 1, col = 'gray')")
289
290 if (standardize=='true'):
291 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 + "')")
292 r("barplot(yz, space = 0, add=TRUE, width = " + half_break_unit_str + ", col=rgb(1, 0, 0, 1))")
293 r("axis(1, hadj = 1, at=seq(0, " + max_break_str + ", by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
294 r("axis(1, at=seq(0, " + max_break_str + ", by= " + break_penta_unit_str + "), labels=TRUE, tcl=-0.5)")
295 r("axis(1, at=seq(0, " + max_break_str + ", by= " + half_break_unit_str + "), labels=FALSE, tcl=-0.25)")
296 elif (standardize=='false'):
297 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 + "')")
298 r("barplot(yz, space = 0, add=TRUE, width = 0.5, col=rgb(1, 0, 0, 1))")
299 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_unit_str + "), labels=FALSE, tcl=-0.5)")
300 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_penta_unit_str + ", labels=TRUE, tcl=-0.5)")
301 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by= " + break_half_unit_str + "), labels=FALSE, tcl=-0.25)")
302
303
304 def build_haw_snp_dictionary(haw_vcf = None):
305 haw_snps = {}
306
307 i_file = open(haw_vcf, 'rU')
308 reader = csv.reader(i_file, delimiter = '\t')
309
310 skip_headers(reader = reader, i_file = i_file)
311
312 for row in reader:
313 #print row
314 chromosome = row[0].upper()
315 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
316 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
317
318 position = row[1]
319 haw_snp_id = row[2]
320 ref_allele = row[3]
321 alt_allele = row[4]
322
323 info = row[7]
324
325 mapping_unit = info.replace("MPU=", "")
326
327 location = chromosome + ":" + position
328 haw_snps[location] = (alt_allele, haw_snp_id, mapping_unit)
329
330 i_file.close()
331
332 return haw_snps
333
334 def calculate_normalized_histogram_bins_per_xbase(pileup_info = None, xbase = 1000000, do_not_normalize = False):
335 normalized_histogram_bins_per_xbase = {}
336
337 ref_snp_count_per_xbase = get_ref_snp_count_per_xbase(pileup_info = pileup_info, xbase = xbase)
338 mean_zero_snp_count_per_chromosome = get_mean_zero_snp_count_per_chromosome(pileup_info = pileup_info, xbase = xbase)
339 zero_snp_count_per_xbase = get_zero_snp_count_per_xbase(pileup_info = pileup_info, xbase = xbase)
340
341 for location in ref_snp_count_per_xbase:
342 chromosome = location.split(':')[0]
343 mean_zero_snp_count = mean_zero_snp_count_per_chromosome[chromosome]
344 ref_snp_count = ref_snp_count_per_xbase[location]
345
346 zero_snp_count = 0
347 if location in zero_snp_count_per_xbase:
348 zero_snp_count = zero_snp_count_per_xbase[location]
349
350 if do_not_normalize == True:
351 normalized_histogram_bins_per_xbase[location] = zero_snp_count
352 else:
353 if zero_snp_count == 0 or ref_snp_count == 0:
354 normalized_histogram_bins_per_xbase[location] = 0
355 elif zero_snp_count == ref_snp_count:
356 normalized_histogram_bins_per_xbase[location] = 0
357 else:
358 normalized_histogram_bins_per_xbase[location] = (Decimal(zero_snp_count) / (Decimal(ref_snp_count)-Decimal(zero_snp_count))) * Decimal(mean_zero_snp_count)
359
360 return normalized_histogram_bins_per_xbase
361
362 def get_ref_snp_count_per_xbase(pileup_info = None, xbase = 1000000):
363 ref_snps_per_xbase = {}
364
365 for location in pileup_info:
366 location_info = location.split(':')
367
368 chromosome = location_info[0].upper()
369 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
370 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
371
372 position = location_info[1]
373 xbase_position = (int(position) / xbase) + 1
374
375 location = chromosome + ":" + str(xbase_position)
376 if location in ref_snps_per_xbase:
377 ref_snps_per_xbase[location] = ref_snps_per_xbase[location] + 1
378 else:
379 ref_snps_per_xbase[location] = 1
380
381 return ref_snps_per_xbase
382
383 def get_mean_zero_snp_count_per_chromosome(pileup_info, xbase = 1000000):
384 sample_snp_count_per_xbase = {}
385
386 for location in pileup_info:
387 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_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 alt_allele_count == 0:
396 if xbase_location in sample_snp_count_per_xbase:
397 sample_snp_count_per_xbase[xbase_location] = sample_snp_count_per_xbase[xbase_location] + 1
398 else:
399 sample_snp_count_per_xbase[xbase_location] = 1
400
401 elif alt_allele_count != 0 and xbase_location not in sample_snp_count_per_xbase:
402 sample_snp_count_per_xbase[xbase_location] = 0
403
404 mean_zero_snp_count_per_chromosome = {}
405 for location in sample_snp_count_per_xbase:
406 chromosome = location.split(':')[0]
407 sample_count = sample_snp_count_per_xbase[location]
408 if chromosome in mean_zero_snp_count_per_chromosome:
409 mean_zero_snp_count_per_chromosome[chromosome].append(sample_count)
410 else:
411 mean_zero_snp_count_per_chromosome[chromosome] = [sample_count]
412
413 for chromosome in mean_zero_snp_count_per_chromosome:
414 summa = sum(mean_zero_snp_count_per_chromosome[chromosome])
415 count = len(mean_zero_snp_count_per_chromosome[chromosome])
416
417 mean_zero_snp_count_per_chromosome[chromosome] = Decimal(summa) / Decimal(count)
418
419 return mean_zero_snp_count_per_chromosome
420
421 def get_zero_snp_count_per_xbase(pileup_info = None, xbase = 1000000):
422 zero_snp_count_per_xbase = {}
423
424 for location in pileup_info:
425 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
426
427 location_info = location.split(':')
428 chromosome = location_info[0]
429 position = location_info[1]
430 xbase_position = (int(position) / xbase) + 1
431 xbase_location = chromosome + ":" + str(xbase_position)
432
433 if alt_allele_count == 0:
434 if xbase_location in zero_snp_count_per_xbase:
435 zero_snp_count_per_xbase[xbase_location] = zero_snp_count_per_xbase[xbase_location] + 1
436 else:
437 zero_snp_count_per_xbase[xbase_location] = 1
438
439 elif alt_allele_count != 0 and xbase_location not in zero_snp_count_per_xbase:
440 zero_snp_count_per_xbase[xbase_location] = 0
441
442 return zero_snp_count_per_xbase
443
444 def parse_pileup(sample_pileup = None, haw_snps = None):
445 i_file = open(sample_pileup, 'rU')
446 reader = csv.reader(i_file, delimiter = '\t', quoting = csv.QUOTE_NONE)
447
448 pileup_info = {}
449
450 for row in reader:
451 chromosome = row[0].upper()
452 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
453 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
454
455 position = row[1]
456 ref_allele = row[2]
457 read_depth = row[3]
458 read_bases = row[4]
459
460 location = chromosome + ":" + position
461 if location in haw_snps:
462 alt_allele, haw_snp_id, mapping_unit = haw_snps[location]
463 ref_allele_count, alt_allele_count = parse_read_bases(read_bases = read_bases, alt_allele = alt_allele)
464
465 if Decimal(read_depth!=0):
466 getcontext().prec = 6
467 ratio = Decimal(alt_allele_count) / Decimal(read_depth)
468 #ratio = Decimal(alt_allele_count) / Decimal(ref_allele_count)
469
470 pileup_info[location] = (alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit)
471
472 #debug line
473 #print chromosome, position, read_depth, ref_allele_count, alt_allele_count, ratio, id
474
475 i_file.close()
476
477 return pileup_info
478
479 def parse_read_bases(read_bases = None, alt_allele = None):
480 read_bases = re.sub('\$', '', read_bases)
481 read_bases = re.sub('\^[^\s]', '', read_bases)
482
483 ref_allele_matches = re.findall("\.|\,", read_bases)
484 ref_allele_count = len(ref_allele_matches)
485
486 alt_allele_matches = re.findall(alt_allele, read_bases, flags = re.IGNORECASE)
487 alt_allele_count = len(alt_allele_matches)
488
489 #debug line
490 #print read_bases, alt_allele, alt_allele_count, ref_allele_count
491
492 return ref_allele_count, alt_allele_count
493
494 if __name__ == "__main__":
495 main()