annotate VDM_Mapping.py @ 6:d1273c7210da draft

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