annotate SNP_Mapping.py @ 15:6dc4cc8da6b4 draft

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