annotate EMS_VariantDensityMapping.py @ 9:58a3878549ef draft

Uploaded
author gregory-minevich
date Mon, 25 Jun 2012 16:09:09 -0400
parents 6d2dbdfa11e3
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
6
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
1 #!/usr/bin/python
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
2
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
3 import re
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
4 import sys
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
5 import optparse
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
6 import csv
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
7 from rpy import *
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
8
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
9 def main():
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
10 parser = optparse.OptionParser()
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
11 parser.add_option('-s', '--snp_vcf', dest = 'snp_vcf', action = 'store', type = 'string', default = None, help = "VCF of SNPs")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
12 parser.add_option('-c', '--hist_color', dest = 'hist_color', action = 'store', type = 'string', default = "darkgray", help = "Color for 1Mb histograms")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
13 parser.add_option('-y', '--ylim', dest = 'ylim', action = 'store', type = 'int', default= 100, help = "Upper limit of Y axis")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
14 parser.add_option('-z', '--standardize', dest = 'standardize', default= 'false', help = "Standardize X-axis")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
15 parser.add_option('-e', '--ems', dest = 'ems', default= 'false', help = "Whether EMS variants should be filtered for")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
16 parser.add_option('-o', '--output', dest = 'plot_output', action = 'store', type = 'string', default = 'EMS_Variant_Density_Plot.pdf', help = "Output file name of plot")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
17 (options, args) = parser.parse_args()
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
18
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
19
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
20 i, ii, iii, iv, v, x = parse_snp_vcf(snp_vcf = options.snp_vcf, ems=options.ems)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
21 create_histograms(plot_output = options.plot_output, hist_color=options.hist_color, ylim=options.ylim, ems=options.ems, standardize=options.standardize, i = i, ii = ii, iii = iii, iv = iv, v = v, x = x)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
22
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
23 def create_histograms(plot_output = None, hist_color=None, ylim=None, ems=None, standardize=None , i = None, ii = None, iii = None, iv = None, v = None, x = None):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
24 breaks = { 'I' : 16 , 'II' : 16, 'III' : 14, 'IV' : 18, 'V' : 21, 'X' : 18 }
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
25
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
26 try:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
27 r.pdf(plot_output, 8, 8)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
28 if len(i) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
29 plot_data(position_list = i, chr = "I", breaks = breaks["I"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
30 if len(ii) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
31 plot_data(position_list = ii, chr = "II", breaks = breaks["II"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
32 if len(iii) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
33 plot_data(position_list = iii, chr = "III", breaks = breaks["III"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
34 if len(iv) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
35 plot_data(position_list = iv, chr = "IV", breaks = breaks["IV"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
36 if len(v) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
37 plot_data(position_list = v, chr = "V", breaks = breaks["V"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
38 if len(x) > 0:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
39 plot_data(position_list = x, chr = "X", breaks = breaks["X"], hist_color=hist_color, ylim=ylim, ems=ems, standardize=standardize)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
40 r.dev_off()
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
41 except Exception as inst:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
42 print inst
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
43 print "There was an error creating the plot pdf... Please try again"
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
44
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
45 def parse_snp_vcf(snp_vcf = None, ems=None):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
46 i_file = open(snp_vcf, 'rU')
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
47 reader = csv.reader(i_file, delimiter = '\t', quoting = csv.QUOTE_NONE)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
48
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
49 skip_headers(reader = reader, i_file = i_file)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
50
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
51 i_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
52 ii_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
53 iii_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
54 iv_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
55 v_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
56 x_position_list = []
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
57
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
58 for row in reader:
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
59 chromosome = row[0].upper()
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
60 chromosome = re.sub("chr", "", chromosome, flags = re.IGNORECASE)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
61 chromosome = re.sub("CHROMOSOME_", "", chromosome, flags = re.IGNORECASE)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
62
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
63 position = row[1]
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
64 ref_allele = row[3]
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
65 alt_allele = row[4]
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
66
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
67 if (ems=='true'):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
68 if (ref_allele =="G" or ref_allele =="C") and (alt_allele =="A" or alt_allele =="T"):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
69 if chromosome == "I":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
70 i_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
71 elif chromosome == "II":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
72 ii_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
73 elif chromosome == "III":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
74 iii_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
75 elif chromosome == "IV":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
76 iv_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
77 elif chromosome == "V":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
78 v_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
79 elif chromosome == "X":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
80 x_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
81 elif (ems=='false'):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
82 if chromosome == "I":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
83 i_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
84 elif chromosome == "II":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
85 ii_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
86 elif chromosome == "III":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
87 iii_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
88 elif chromosome == "IV":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
89 iv_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
90 elif chromosome == "V":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
91 v_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
92 elif chromosome == "X":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
93 x_position_list.append(position)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
94
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
95 return i_position_list, ii_position_list, iii_position_list, iv_position_list, v_position_list, x_position_list
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
96
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
97 def skip_headers(reader = None, i_file = None):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
98 # count headers
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
99 comment = 0
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
100 while reader.next()[0].startswith('#'):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
101 comment = comment + 1
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
102
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
103 # skip headers
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
104 i_file.seek(0)
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
105 for i in range(0, comment):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
106 reader.next()
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
107
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
108 def plot_data(position_list = None, chr = None, breaks = None, hist_color=None, ylim = None, ems=None, standardize=None):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
109 positions = ",".join(map(str, map(lambda x: float(x) / 1000000, position_list)))
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
110 positions = "c(" + positions + ")"
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
111
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
112 if (standardize=='true'):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
113 r("hist(" + positions + ", xlim=c(0,21), ylim=c(0, %d "%ylim +"),col='"+ hist_color + "', breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=1), main = 'LG " + chr + "', ylab = 'Frequency Of SNPs', xlab = 'Location (Mb)')")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
114 r("hist(" + positions + ", xlim=c(0,21), add=TRUE, ylim=c(0, %d "%ylim +"), col=rgb(1, 0, 0, 1), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=.5), main = 'Chr " + chr + "', ylab = 'Number Of SNPs', xlab = 'Location (Mb)')")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
115 r("axis(1, at=seq(0, 21, by=1), labels=FALSE, tcl=-0.5)")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
116 r("axis(1, at=seq(0, 21, by=0.5), labels=FALSE, tcl=-0.25)")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
117 elif (standardize=='false'):
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
118 r("hist(" + positions + ", xlim=c(0,as.integer( ' " + str(breaks) + " ')), ylim=c(0, %d "%ylim +"),col='"+ hist_color + "', breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=1), main = 'LG " + chr + "', ylab = 'Frequency Of SNPs', xlab = 'Location (Mb)')")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
119 r("hist(" + positions + ", xlim=c(0,as.integer( ' " + str(breaks) + " ')), add=TRUE, ylim=c(0, %d "%ylim +"), col=rgb(1, 0, 0, 1), breaks = seq(0, as.integer( ' " + str(breaks) + " '), by=.5), main = 'Chr " + chr + "', ylab = 'Number Of SNPs', xlab = 'Location (Mb)')")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
120 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=1), labels=FALSE, tcl=-0.5)")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
121 r("axis(1, at=seq(0, as.integer( ' " + str(breaks) + " '), by=0.5), labels=FALSE, tcl=-0.25)")
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
122
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
123
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
124
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
125 if __name__ == "__main__":
6d2dbdfa11e3 Uploaded
gregory-minevich
parents:
diff changeset
126 main()