annotate hetbox.py @ 4:dfa2e75da6aa draft

fixed script typo
author nick
date Thu, 30 May 2013 15:41:27 -0400
parents 3a1ce69571e5
children 99cda2ff12b8
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
1 #!/usr/bin/env python
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
2 # New in this version:
4
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
3 # - option to generate report is triggered simply by including report filename
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
4 # as third argument
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
5
128db16c9399 Uploaded script
nick
parents:
diff changeset
6 import os,sys,numpy
128db16c9399 Uploaded script
nick
parents:
diff changeset
7 from rpy2.robjects import Formula
128db16c9399 Uploaded script
nick
parents:
diff changeset
8 from rpy2.robjects.packages import importr
128db16c9399 Uploaded script
nick
parents:
diff changeset
9 from rpy2 import robjects
128db16c9399 Uploaded script
nick
parents:
diff changeset
10
128db16c9399 Uploaded script
nick
parents:
diff changeset
11 def fail(message):
128db16c9399 Uploaded script
nick
parents:
diff changeset
12 sys.stderr.write(message+'\n')
128db16c9399 Uploaded script
nick
parents:
diff changeset
13 sys.exit(1)
128db16c9399 Uploaded script
nick
parents:
diff changeset
14
128db16c9399 Uploaded script
nick
parents:
diff changeset
15 args = sys.argv[1:]
128db16c9399 Uploaded script
nick
parents:
diff changeset
16 if len(args) >= 1:
128db16c9399 Uploaded script
nick
parents:
diff changeset
17 infile = args[0]
128db16c9399 Uploaded script
nick
parents:
diff changeset
18 else:
128db16c9399 Uploaded script
nick
parents:
diff changeset
19 fail('Error: No input filename provided (as argument 1).')
128db16c9399 Uploaded script
nick
parents:
diff changeset
20 if len(args) >= 2:
128db16c9399 Uploaded script
nick
parents:
diff changeset
21 outfile = args[1]
128db16c9399 Uploaded script
nick
parents:
diff changeset
22 else:
128db16c9399 Uploaded script
nick
parents:
diff changeset
23 fail('Error: No output filename provided (as argument 2).')
128db16c9399 Uploaded script
nick
parents:
diff changeset
24 if len(args) >= 3:
128db16c9399 Uploaded script
nick
parents:
diff changeset
25 report = args[2]
4
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
26 else:
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
27 report = ''
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
28
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
29 # Check input file
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
30 if not os.path.exists(infile):
128db16c9399 Uploaded script
nick
parents:
diff changeset
31 fail('Error: Input file '+infile+' could not be found.')
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
32 with open(infile, 'r') as lines:
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
33 line = lines.readline()
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
34 if not line:
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
35 fail('Error: Input file seems to be empty')
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
36 line = line.strip().lstrip('#') # rm whitespace, comment chars
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
37 labels = line.split("\t")
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
38 if 'SAMPLE' not in labels or labels[11] != 'MINOR.FREQ.PERC.':
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
39 fail('Error: Input file does not seem to have a proper header line.')
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
40
128db16c9399 Uploaded script
nick
parents:
diff changeset
41
128db16c9399 Uploaded script
nick
parents:
diff changeset
42 utils = importr('utils')
128db16c9399 Uploaded script
nick
parents:
diff changeset
43 graphics = importr('graphics')
128db16c9399 Uploaded script
nick
parents:
diff changeset
44 base = importr('base')
128db16c9399 Uploaded script
nick
parents:
diff changeset
45 stats = importr('stats')
128db16c9399 Uploaded script
nick
parents:
diff changeset
46 rprint = robjects.globalenv.get("print")
128db16c9399 Uploaded script
nick
parents:
diff changeset
47 grdevices = importr('grDevices')
128db16c9399 Uploaded script
nick
parents:
diff changeset
48 grdevices.png(file=outfile, width=1024, height=768)
128db16c9399 Uploaded script
nick
parents:
diff changeset
49
128db16c9399 Uploaded script
nick
parents:
diff changeset
50 # Read file into a data frame
128db16c9399 Uploaded script
nick
parents:
diff changeset
51 DATA = utils.read_delim(infile)
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
52 # Remove comment from header, if
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
53 labels = robjects.r.names(DATA)
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
54 if labels[0][0:2] == 'X.':
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
55 labels[0] = labels[0][2:]
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
56 #robjects.r.assign('data', DATA)
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
57 #robjects.r('data$MINOR.FREQ.PERC. = data$MINOR.FREQ.PERC. * 100')
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
58 #DATA = robjects.r('data')
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
59 #index = data.names.index('MINOR.FREQ.PERC.')
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
60 # Multiply minor allele frequencies by 100 to get percentage
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
61 # .rx2() looks up a column by its label and returns it as a vector
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
62 # .ro turns the returned object into one that can be operated on per-element
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
63 minor_freq = DATA.rx2('MINOR.FREQ.PERC.').ro * 100
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
64 samples = DATA.rx2('SAMPLE')
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
65 # Formula() creates a Python object representing the R object returned by x ~ y
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
66 formula = Formula('minor_freq ~ samples')
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
67 # The "environment" in .getenvironment() is the entire R workspace in which the
128db16c9399 Uploaded script
nick
parents:
diff changeset
68 # Formula object exists. The R workspace meaning all the defined variables.
128db16c9399 Uploaded script
nick
parents:
diff changeset
69 # Here, the .getenvironment() method is being used to set some variables in the
128db16c9399 Uploaded script
nick
parents:
diff changeset
70 # R workspace
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
71 formula.getenvironment()['minor_freq'] = minor_freq
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
72 formula.getenvironment()['samples'] = samples
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
73
128db16c9399 Uploaded script
nick
parents:
diff changeset
74 # create boxplot - fill kwargs1 with the options for the boxplot function
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
75 kwargs1 = {'ylab':"Minor allele frequency (%)", 'col':"gray", 'xaxt':"n", 'outpch':"*",'main':"Distribution of minor allele frequencies >= 2%", 'cex.lab':"1.5"}
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
76 p = graphics.boxplot(formula, **kwargs1)
128db16c9399 Uploaded script
nick
parents:
diff changeset
77
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
78 table = base.table(DATA.rx2('SAMPLE'))
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
79 graphics.text(0.5, 1, 'N:', font=2)
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
80 for i in range(1, base.length(table)[0]+1, 1):
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
81 graphics.text(i, 1, table[i-1], font=2)
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
82
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
83 graphlabels = base.names(table)
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
84 kwargs3 = {'pos':"0", 'las':"2", 'cex.axis':"1"}
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
85 graphics.axis(1, at=range(1, len(graphlabels)+1, 1), labels=graphlabels, **kwargs3)
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
86 grdevices.dev_off()
128db16c9399 Uploaded script
nick
parents:
diff changeset
87
4
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
88 if not report:
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
89 sys.exit(0)
128db16c9399 Uploaded script
nick
parents:
diff changeset
90
2
3a1ce69571e5 Allow commented lines, de-kludge data frame handling
nick
parents: 0
diff changeset
91
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
92 ####################################
128db16c9399 Uploaded script
nick
parents:
diff changeset
93 # GENERATE REPORT
128db16c9399 Uploaded script
nick
parents:
diff changeset
94 # report should be something like:
128db16c9399 Uploaded script
nick
parents:
diff changeset
95 # SAMPLE NoHET MEDIAN MAD TEST
128db16c9399 Uploaded script
nick
parents:
diff changeset
96 # s1 7 10% n p/w/f
128db16c9399 Uploaded script
nick
parents:
diff changeset
97 # n <= 5 pass
128db16c9399 Uploaded script
nick
parents:
diff changeset
98 # 6 <= n <=10 warn
128db16c9399 Uploaded script
nick
parents:
diff changeset
99 # n >= 11 fail
128db16c9399 Uploaded script
nick
parents:
diff changeset
100 # MAD <= 2.0 fail
128db16c9399 Uploaded script
nick
parents:
diff changeset
101 # MAD > 2.0 pass
128db16c9399 Uploaded script
nick
parents:
diff changeset
102 ###################################
128db16c9399 Uploaded script
nick
parents:
diff changeset
103
128db16c9399 Uploaded script
nick
parents:
diff changeset
104 SAMPLES=[]
4
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
105 for i in range(len(table)):
dfa2e75da6aa fixed script typo
nick
parents: 2
diff changeset
106 SAMPLES.append(base.names(table)[i])
0
128db16c9399 Uploaded script
nick
parents:
diff changeset
107
128db16c9399 Uploaded script
nick
parents:
diff changeset
108 def boxstats(data,sample):
128db16c9399 Uploaded script
nick
parents:
diff changeset
109 VALUES = [100*float(x.strip().split('\t')[11]) for x in list(open(data)) if x.strip().split('\t')[0]==sample]
128db16c9399 Uploaded script
nick
parents:
diff changeset
110 NoHET = len(VALUES)
128db16c9399 Uploaded script
nick
parents:
diff changeset
111 MEDIAN = numpy.median(VALUES)
128db16c9399 Uploaded script
nick
parents:
diff changeset
112 MAD = numpy.median([abs(i - MEDIAN) for i in VALUES]) # Median absolute distance (robust spread statistic)
128db16c9399 Uploaded script
nick
parents:
diff changeset
113 return [NoHET,MEDIAN, MAD]
128db16c9399 Uploaded script
nick
parents:
diff changeset
114
128db16c9399 Uploaded script
nick
parents:
diff changeset
115 boxreport = open(report, "w+")
128db16c9399 Uploaded script
nick
parents:
diff changeset
116 boxreport.write("SAMPLE\tTOTAL.SITES\tMEDIAN.FREQ.\tMAD.FREQ\tEVAL\n")
128db16c9399 Uploaded script
nick
parents:
diff changeset
117 for sample in SAMPLES:
128db16c9399 Uploaded script
nick
parents:
diff changeset
118 ENTRY = [sample] + boxstats(infile,sample)
128db16c9399 Uploaded script
nick
parents:
diff changeset
119 if ENTRY[1] <= 5:
128db16c9399 Uploaded script
nick
parents:
diff changeset
120 ENTRY.append('pass')
128db16c9399 Uploaded script
nick
parents:
diff changeset
121 elif 6 <= ENTRY[1] <=10:
128db16c9399 Uploaded script
nick
parents:
diff changeset
122 ENTRY.append('warn')
128db16c9399 Uploaded script
nick
parents:
diff changeset
123 elif ENTRY[1] >= 11:
128db16c9399 Uploaded script
nick
parents:
diff changeset
124 ENTRY.append('fail')
128db16c9399 Uploaded script
nick
parents:
diff changeset
125 if ENTRY[3] <=2.0:
128db16c9399 Uploaded script
nick
parents:
diff changeset
126 ENTRY.append('fail')
128db16c9399 Uploaded script
nick
parents:
diff changeset
127 elif ENTRY[3] >2.0:
128db16c9399 Uploaded script
nick
parents:
diff changeset
128 ENTRY.append('pass')
128db16c9399 Uploaded script
nick
parents:
diff changeset
129 if len(set(ENTRY[4:6])) == 2:
128db16c9399 Uploaded script
nick
parents:
diff changeset
130 ENTRY.append('warn')
128db16c9399 Uploaded script
nick
parents:
diff changeset
131 else:
128db16c9399 Uploaded script
nick
parents:
diff changeset
132 ENTRY.append(list(set(ENTRY[4:6]))[0])
128db16c9399 Uploaded script
nick
parents:
diff changeset
133 boxreport.write ('%s\t%d\t%.1f\t%.1f\t%s\n' % tuple([ENTRY[i] for i in [0,1,2,3,6]]))
128db16c9399 Uploaded script
nick
parents:
diff changeset
134
128db16c9399 Uploaded script
nick
parents:
diff changeset
135 boxreport.close()
128db16c9399 Uploaded script
nick
parents:
diff changeset
136
128db16c9399 Uploaded script
nick
parents:
diff changeset
137
128db16c9399 Uploaded script
nick
parents:
diff changeset
138
128db16c9399 Uploaded script
nick
parents:
diff changeset
139
128db16c9399 Uploaded script
nick
parents:
diff changeset
140
128db16c9399 Uploaded script
nick
parents:
diff changeset
141
128db16c9399 Uploaded script
nick
parents:
diff changeset
142 #####################################
128db16c9399 Uploaded script
nick
parents:
diff changeset
143 #STUFF
128db16c9399 Uploaded script
nick
parents:
diff changeset
144 #
128db16c9399 Uploaded script
nick
parents:
diff changeset
145 #kwargs = {'ylab':"Minor allele frequency (sites >= 2%)",'col':"blue", 'axes':"FALSE", 'outpch':"*", 'ylim':"c(-0.03,0.5)", 'main':"Minor allele frequencies run020 at 2%"}
128db16c9399 Uploaded script
nick
parents:
diff changeset
146 #boxplot(freq~id,data=run020[run020$freq>=0.02,], col="blue", axes=FALSE, outpch="*", ylab="Minor allele frequency (sites >= 2%)", ylim=c(-0.03,0.5), main="Minor allele frequencies run020 at 2%")
128db16c9399 Uploaded script
nick
parents:
diff changeset
147 #biglabels2=sort(as.vector(unique(run020$id)))
128db16c9399 Uploaded script
nick
parents:
diff changeset
148 #biglabels=apply(as.matrix(biglabels2), 1,function(x) substr(x,1,5))
128db16c9399 Uploaded script
nick
parents:
diff changeset
149 #axis(1, at=c(1:length(biglabels)),labels=biglabels, pos=-0.02, las=2, cex.axis=0.9)
128db16c9399 Uploaded script
nick
parents:
diff changeset
150 #axis(2, at=c(seq(0,50,5)/100), pos=0,cex.axis=0.9)
128db16c9399 Uploaded script
nick
parents:
diff changeset
151 #nbig=as.vector(table(run020[run020$freq>=0.02,1]))
128db16c9399 Uploaded script
nick
parents:
diff changeset
152 #for (i in 1:length(nbig)){text(i,-0.01,nbig[i],cex=0.9, font=2)}
128db16c9399 Uploaded script
nick
parents:
diff changeset
153