annotate SNP_Mapping.py @ 1:619b943864a9

Uploaded
author gregory-minevich
date Tue, 20 Mar 2012 10:40:15 -0400
parents 30fa4d84e84c
children 7d6bccd5f88c
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
1 #!/usr/bin/python
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
2
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
3 import re
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
4 import sys
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
5 import optparse
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
6 import csv
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
7 import re
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
8 from decimal import *
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
9 from rpy import *
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
10
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
11 def main():
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
12 csv.field_size_limit(1000000000)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
13
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
14 parser = optparse.OptionParser()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
15 parser.add_option('-p', '--sample_pileup', dest = 'sample_pileup', action = 'store', type = 'string', default = None, help = "Sample pileup from mpileup")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
16 parser.add_option('-v', '--haw_vcf', dest = 'haw_vcf', action = 'store', type = 'string', default = None, help = "vcf file of Hawaiian SNPs")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
17 parser.add_option('-l', '--loess_span', dest = 'loess_span', action = 'store', type = 'float', default = .01, help = "Loess span")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
18 parser.add_option('-d', '--d_yaxis', dest = 'd_yaxis', action = 'store', type = 'float', default = .7, help = "y-axis upper limit for dot plot")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
19 parser.add_option('-y', '--h_yaxis', dest = 'h_yaxis', action = 'store', type = 'int', default = 500, help = "y-axis upper limit for dot plot")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
20 parser.add_option('-c', '--points_color', dest = 'points_color', action = 'store', type = 'string', default = "gray27", help = "Color for data points")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
21 parser.add_option('-k', '--loess_color', dest = 'loess_color', action = 'store', type = 'string', default = "red", help = "Color for loess regression line")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
22 parser.add_option('-z', '--standardize', dest = 'standardize', default= 'false', help = "Standardize X-axis")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
23
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
24 parser.add_option('-o', '--output', dest = 'output', action = 'store', type = 'string', default = None, help = "Output file name")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
25 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")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
26
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
27 #For plotting with map units on the X-axis instead of physical distance
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
28 #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")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
29 (options, args) = parser.parse_args()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
30
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
31 haw_snps = build_haw_snp_dictionary(haw_vcf = options.haw_vcf)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
32 pileup_info = parse_pileup(sample_pileup = options.sample_pileup, haw_snps = haw_snps)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
33 output_pileup_info(output = options.output, pileup_info = pileup_info)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
34
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
35 #output plot with all ratios
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
36 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)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
37
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
38 #For plotting with map units on the X-axis instead of physical distance)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
39 #output_scatter_plots_by_mapping_units(mpu_plot_output = options.mpu_plot_output, pileup_info = pileup_info)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
40
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
41 def skip_headers(reader = None, i_file = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
42 # count headers
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
43 comment = 0
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
44 while reader.next()[0].startswith('#'):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
45 comment = comment + 1
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
46
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
47 # skip headers
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
48 i_file.seek(0)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
49 for i in range(0, comment):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
50 reader.next()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
51
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
52 def location_comparer(location_1, location_2):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
53 chr_loc_1 = location_1.split(':')[0]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
54 pos_loc_1 = int(location_1.split(':')[1])
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
55
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
56 chr_loc_2 = location_2.split(':')[0]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
57 pos_loc_2 = int(location_2.split(':')[1])
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
58
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
59 if chr_loc_1 == chr_loc_2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
60 if pos_loc_1 < pos_loc_2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
61 return -1
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
62 elif pos_loc_1 == pos_loc_1:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
63 return 0
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
64 elif pos_loc_1 > pos_loc_2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
65 return 1
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
66 elif chr_loc_1 < chr_loc_2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
67 return -1
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
68 elif chr_loc_1 > chr_loc_2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
69 return 1
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
70
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
71 def output_pileup_info(output = None, pileup_info = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
72 o_file = open(output, 'wb')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
73 writer = csv.writer(o_file, delimiter = '\t')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
74
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
75 writer.writerow(["#Chr\t", "Pos\t", "ID\t", "Alt Count\t", "Ref Count\t", "Read Depth\t", "Ratio\t", "Mapping Unit"])
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
76
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
77 location_sorted_pileup_info_keys = sorted(pileup_info.keys(), cmp=location_comparer)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
78
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
79 for location in location_sorted_pileup_info_keys:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
80 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
81
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
82 location_info = location.split(':')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
83 chromosome = location_info[0]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
84 position = location_info[1]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
85
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
86 writer.writerow([chromosome, position, haw_snp_id, alt_allele_count, ref_allele_count, read_depth, ratio, mapping_unit])
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
87
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
88 o_file.close()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
89
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
90 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):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
91 i = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
92 ii = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
93 iii = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
94 iv = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
95 v = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
96 x = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
97
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
98 breaks = { 'I' : 16 , 'II' : 16, 'III' : 14, 'IV' : 18, 'V' : 21, 'X' : 18 }
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
99
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
100 for location in pileup_info:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
101 chromosome = location.split(':')[0]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
102 position = location.split(':')[1]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
103
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
104 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
105
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
106 if chromosome == "I":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
107 i[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
108 elif chromosome == "II":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
109 ii[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
110 elif chromosome == "III":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
111 iii[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
112 elif chromosome == "IV":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
113 iv[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
114 elif chromosome == "V":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
115 v[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
116 elif chromosome == "X":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
117 x[position] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
118
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
119 x_label = "Location (Mb)"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
120 filtered_label = ''
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
121
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
122
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
123 try:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
124 r.pdf(location_plot_output, 8, 8)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
125 if i:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
126 plot_data(chr_dict = i, chr = "I" + 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["I"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
127 if ii:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
128 plot_data(chr_dict = ii, chr = "II" + 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["II"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
129 if iii:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
130 plot_data(chr_dict = iii, chr = "III" + 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["III"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
131 if iv:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
132 plot_data(chr_dict = iv, chr = "IV" + 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["IV"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
133 if v:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
134 plot_data(chr_dict = v, chr = "V" + 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["V"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
135 if x:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
136 plot_data(chr_dict = x, chr = "X" + 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["X"], standardize=standardize)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
137
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
138 r.dev_off()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
139 except Exception as inst:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
140 print inst
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
141 print "There was an error creating the location plot pdf... Please try again"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
142
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
143 def output_scatter_plots_by_mapping_units(mpu_plot_output = None, pileup_info = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
144 i = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
145 ii = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
146 iii = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
147 iv = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
148 v = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
149 x = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
150
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
151 for location in pileup_info:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
152 chromosome = location.split(':')[0]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
153 position = location.split(':')[1]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
154
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
155 alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit = pileup_info[location]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
156
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
157 if chromosome == "I":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
158 i[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
159 elif chromosome == "II":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
160 ii[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
161 elif chromosome == "III":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
162 iii[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
163 elif chromosome == "IV":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
164 iv[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
165 elif chromosome == "V":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
166 v[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
167 elif chromosome == "X":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
168 x[mapping_unit] = ratio
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
169
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
170 x_label = "Map Units"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
171
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
172 try:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
173 r.pdf(mpu_plot_output, 8, 8)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
174 plot_data(chr_dict = i, chr = "I", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
175 plot_data(chr_dict = ii, chr = "II", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
176 plot_data(chr_dict = iii, chr = "III", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
177 plot_data(chr_dict = iv, chr = "IV", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
178 plot_data(chr_dict = v, chr = "V", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
179 plot_data(chr_dict = x, chr = "X", x_label = "Map Units")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
180 r.dev_off()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
181 except Exception as inst:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
182 print inst
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
183 print "There was an error creating the map unit plot pdf... Please try again"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
184
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
185
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
186 def plot_data(chr_dict = 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):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
187 ratios = "c("
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
188 positions = "c("
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
189 z_ratios = "c("
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
190 z_positions = "c("
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
191
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
192 for position in chr_dict:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
193 ratio = chr_dict[position]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
194 if divide_position:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
195 position = float(position) / 1000000.0
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
196 positions = positions + str(position) + ", "
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
197 ratios = ratios + str(ratio) + ", "
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
198 if ratio == 0:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
199 if divide_position:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
200 z_position = float(position) / 1000000.0
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
201 z_positions = z_positions + str(position) + ", "
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
202 z_ratios = z_ratios + str(ratio) + ", "
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
203
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
204
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
205 if len(ratios) == 2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
206 ratios = ratios + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
207 else:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
208 ratios = ratios[0:len(ratios) - 2] + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
209
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
210 if len(z_ratios) == 2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
211 z_ratios = z_ratios + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
212 else:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
213 z_ratios = z_ratios[0:len(z_ratios) - 2] + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
214
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
215
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
216 if len(positions) == 2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
217 positions = positions + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
218 else:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
219 positions = positions[0:len(positions) - 2] + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
220
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
221 if len(z_positions) == 2:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
222 z_positions = z_positions + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
223 else:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
224 z_positions = z_positions[0:len(z_positions) - 2] + ")"
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
225
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
226 r("x <- " + positions)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
227 r("y <- " + ratios)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
228
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
229 r("xz <- " + z_positions)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
230 r("yz <- " + z_ratios)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
231
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
232 if (standardize=='true'):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
233 r("plot(x, y, cex=0.60, xlim=c(0,21), main='LG " + chr + "', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='HA Ratio [Hawaiian Reads / Total Read Depth]', pch=18, col='"+ points_color +"')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
234 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
235 r("axis(1, at=seq(0, 21, by=1), labels=FALSE, tcl=-0.5)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
236 r("axis(1, at=seq(0, 21, by=0.5), labels=FALSE, tcl=-0.25)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
237 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
238 elif (standardize=='false'):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
239 r("plot(x, y, cex=0.60, main='LG " + chr + "', xlab= '" + x_label + "', ylim = c(0, %f " %d_yaxis + "), ylab='HA Ratio [Hawaiian Reads / Total Read Depth]', pch=18, col='"+ points_color +"')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
240 r("lines(loess.smooth(x, y, span = %f "%loess_span + "), lwd=5, col='"+ loess_color +"')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
241 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=1), labels=FALSE, tcl=-0.5)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
242 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=0.5), labels=FALSE, tcl=-0.25)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
243 r("axis(2, at=seq(floor(min(y)), 1, by=0.1), labels=FALSE, tcl=-0.2)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
244
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
245
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
246 if draw_secondary_grid_lines:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
247 r("abline(h = seq(floor(min(y)), 1, by=0.1), v = seq(floor(min(x)), length(x), by= 1), col='gray')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
248 else:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
249 r("grid(lty = 1, col = 'gray')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
250
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
251 if (standardize=='true'):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
252 r("hist(xz, col='darkgray', xlim=c(0,21), xlab='Location (Mb)', ylab='Frequency SNP Positions Where HA Ratio=0 ', ylim=c(0, %f " %h_yaxis + "), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=1), main='LG " + chr + "')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
253 r("hist(xz, add=TRUE, col=rgb(1, 0, 0, 1), xlim=c(0,21), xlab='Location (Mb)', ylab='Frequency SNP Positions Where HA Ratio=0 ', ylim=c(0, %f " %h_yaxis + "), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=.5), main='Chr " + chr + "')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
254 r("axis(1, at=seq(0, 21, by=1), labels=FALSE, tcl=-0.5)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
255 r("axis(1, at=seq(0, 21, by=0.5), labels=FALSE, tcl=-0.25)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
256 elif (standardize=='false'):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
257 r("hist(xz, col='darkgray', xlab='Location (Mb)', ylab='Frequency SNP Positions Where HA Ratio=0 ', ylim=c(0, %f " %h_yaxis + "), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=1), main='LG " + chr + "')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
258 r("hist(xz, add=TRUE, col=rgb(1, 0, 0, 1), xlab='Location (Mb)', ylab='Frequency SNP Positions Where HA Ratio=0 ', ylim=c(0, %f " %h_yaxis + "), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=.5), main='Chr " + chr + "')")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
259 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=1), labels=FALSE, tcl=-0.5)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
260 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=0.5), labels=FALSE, tcl=-0.25)")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
261
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
262
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
263 def build_haw_snp_dictionary(haw_vcf = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
264 haw_snps = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
265
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
266 i_file = open(haw_vcf, 'rU')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
267 reader = csv.reader(i_file, delimiter = '\t')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
268
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
269 skip_headers(reader = reader, i_file = i_file)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
270
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
271 for row in reader:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
272 #print row
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
273 chromosome = row[0].upper()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
274 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
275 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
276
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
277 position = row[1]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
278 haw_snp_id = row[2]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
279 ref_allele = row[3]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
280 alt_allele = row[4]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
281 info = row[7]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
282
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
283 mapping_unit = info.replace("MPU=", "")
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
284
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
285 location = chromosome + ":" + position
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
286 haw_snps[location] = (alt_allele, haw_snp_id, mapping_unit)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
287
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
288 i_file.close()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
289
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
290 return haw_snps
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
291
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
292 def parse_pileup(sample_pileup = None, haw_snps = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
293 i_file = open(sample_pileup, 'rU')
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
294 reader = csv.reader(i_file, delimiter = '\t', quoting = csv.QUOTE_NONE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
295
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
296 pileup_info = {}
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
297
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
298 for row in reader:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
299 chromosome = row[0].upper()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
300 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
301 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
302
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
303 position = row[1]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
304 ref_allele = row[2]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
305 read_depth = row[3]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
306 read_bases = row[4]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
307
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
308 location = chromosome + ":" + position
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
309 if location in haw_snps:
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
310 alt_allele, haw_snp_id, mapping_unit = haw_snps[location]
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
311 ref_allele_count, alt_allele_count = parse_read_bases(read_bases = read_bases, alt_allele = alt_allele)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
312
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
313 getcontext().prec = 6
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
314 ratio = Decimal(alt_allele_count) / Decimal(read_depth)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
315
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
316 pileup_info[location] = (alt_allele_count, ref_allele_count, read_depth, ratio, haw_snp_id, mapping_unit)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
317
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
318 #debug line
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
319 #print chromosome, position, read_depth, ref_allele_count, alt_allele_count, ratio, id
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
320
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
321 i_file.close()
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
322
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
323 return pileup_info
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
324
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
325 def parse_read_bases(read_bases = None, alt_allele = None):
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
326 read_bases = re.sub('\$', '', read_bases)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
327 read_bases = re.sub('\^[^\s]', '', read_bases)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
328
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
329 ref_allele_matches = re.findall("\.|\,", read_bases)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
330 ref_allele_count = len(ref_allele_matches)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
331
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
332 alt_allele_matches = re.findall(alt_allele, read_bases, flags = re.IGNORECASE)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
333 alt_allele_count = len(alt_allele_matches)
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
334
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
335 #debug line
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
336 #print read_bases, alt_allele, alt_allele_count, ref_allele_count
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
337
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
338 return ref_allele_count, alt_allele_count
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
339
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
340 if __name__ == "__main__":
30fa4d84e84c Uploaded
gregory-minevich
parents:
diff changeset
341 main()