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