|
0
|
1 #!/usr/local/bin/python
|
|
|
2 # This is a truly ghastly hack
|
|
|
3 # all of the heavy data cleaning lifting is done in R which is a really dumb place IMHO
|
|
|
4 # Making a new file seems a waste but it would be far easier to set everything up in python
|
|
|
5 # seems to work so I'm leaving it alone
|
|
|
6 # sigh. Should really move this gig to rpy - writing a robust R script is hard.
|
|
|
7 # updated to compress pdf using gs since millions of points = horsechoker pdfs and pdfs are good
|
|
|
8 # updated july 20 to fix sort order - R unique() sorts into strict collating order
|
|
|
9 # so need to sort after unique to revert to lexicographic order for x axis on Manhattan
|
|
|
10 # rgmanqq updated july 19 to deal with x,y and mt
|
|
|
11 # lots of fixes
|
|
|
12 # ross lazarus
|
|
|
13 import sys,math,shutil,subprocess,os,time,tempfile,string
|
|
|
14 from os.path import abspath
|
|
|
15 from rgutils import timenow, RRun, galhtmlprefix, galhtmlpostfix, galhtmlattr
|
|
|
16 progname = os.path.split(sys.argv[0])[1]
|
|
|
17 myversion = 'V000.1 March 2010'
|
|
|
18 verbose = False
|
|
|
19 debug = False
|
|
|
20
|
|
|
21 rcode="""
|
|
|
22 # generalised so 3 core fields passed as parameters ross lazarus March 24 2010 for rgenetics
|
|
|
23 # Originally created as qqman with the following
|
|
|
24 # attribution:
|
|
|
25 #--------------
|
|
|
26 # Stephen Turner
|
|
|
27 # http://StephenTurner.us/
|
|
|
28 # http://GettingGeneticsDone.blogspot.com/
|
|
|
29
|
|
|
30 # Last updated: 19 July 2011 by Ross Lazarus
|
|
|
31 # R code for making manhattan plots and QQ plots from plink output files.
|
|
|
32 # With GWAS data this can take a lot of memory. Recommended for use on
|
|
|
33 # 64bit machines only, for now.
|
|
|
34
|
|
|
35 #
|
|
|
36
|
|
|
37 library(ggplot2)
|
|
|
38
|
|
|
39 coloursTouse = c('firebrick','darkblue','goldenrod','darkgreen')
|
|
|
40 # not too ugly but need a colour expert please...
|
|
|
41
|
|
|
42
|
|
|
43 DrawManhattan = function(pvals=Null,chrom=Null,offset=Null,title=NULL, max.y="max",suggestiveline=0, genomewide=T, size.x.labels=9,
|
|
|
44 size.y.labels=10, annotate=F, SNPlist=NULL,grey=0) {
|
|
|
45 if (annotate & is.null(SNPlist)) stop("You requested annotation but provided no SNPlist!")
|
|
|
46 genomewideline=NULL # was genomewideline=-log10(5e-8)
|
|
|
47 n = length(pvals)
|
|
|
48 if (genomewide) { # use bonferroni since might be only a small region?
|
|
|
49 genomewideline = -log10(0.05/n) }
|
|
|
50 offset = as.integer(offset)
|
|
|
51 if (n > 1000000) { offset = offset/10000 }
|
|
|
52 else if (n > 10000) { offset = offset/1000}
|
|
|
53 chro = as.integer(chrom) # already dealt with X and friends?
|
|
|
54 pvals = as.double(pvals)
|
|
|
55 d=data.frame(CHR=chro,BP=offset,P=pvals)
|
|
|
56 if ("CHR" %in% names(d) & "BP" %in% names(d) & "P" %in% names(d) ) {
|
|
|
57 d=d[!is.na(d$P), ]
|
|
|
58 d=d[!is.na(d$BP), ]
|
|
|
59 d=d[!is.na(d$CHR), ]
|
|
|
60 #limit to only chrs 1-22, x=23,y=24,Mt=25?
|
|
|
61 d=d[d$CHR %in% 1:25, ]
|
|
|
62 d=d[d$P>0 & d$P<=1, ]
|
|
|
63 d$logp = as.double(-log10(d$P))
|
|
|
64 dlen = length(d$P)
|
|
|
65 d$pos=NA
|
|
|
66 ticks=NULL
|
|
|
67 lastbase=0
|
|
|
68 chrlist = unique(d$CHR)
|
|
|
69 chrlist = as.integer(chrlist)
|
|
|
70 chrlist = sort(chrlist) # returns lexical ordering
|
|
|
71 if (max.y=="max") { maxy = ceiling(max(d$logp)) }
|
|
|
72 else { maxy = max.y }
|
|
|
73 nchr = length(chrlist) # may be any number?
|
|
|
74 maxy = max(maxy,1.1*genomewideline)
|
|
|
75 if (nchr >= 2) {
|
|
|
76 for (x in c(1:nchr)) {
|
|
|
77 i = chrlist[x] # need the chrom number - may not == index
|
|
|
78 if (x == 1) { # first time
|
|
|
79 d[d$CHR==i, ]$pos = d[d$CHR==i, ]$BP # initialize to first BP of chr1
|
|
|
80 dsub = subset(d,CHR==i)
|
|
|
81 dlen = length(dsub$P)
|
|
|
82 lastbase = max(dsub$pos) # last one
|
|
|
83 tks = d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1]
|
|
|
84 lastchr = i
|
|
|
85 } else {
|
|
|
86 d[d$CHR==i, ]$pos = d[d$CHR==i, ]$BP+lastbase # one humongous contig
|
|
|
87 if (sum(is.na(lastchr),is.na(lastbase),is.na(d[d$CHR==i, ]$pos))) {
|
|
|
88 cat(paste('manhattan: For',title,'chrlistx=',i,'lastchr=',lastchr,'lastbase=',lastbase,'pos=',d[d$CHR==i,]$pos))
|
|
|
89 }
|
|
|
90 tks=c(tks, d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1])
|
|
|
91 lastchr = i
|
|
|
92 dsub = subset(d,CHR==i)
|
|
|
93 lastbase = max(dsub$pos) # last one
|
|
|
94 }
|
|
|
95 ticklim=c(min(d$pos),max(d$pos))
|
|
|
96 xlabs = chrlist
|
|
|
97 }
|
|
|
98 } else { # nchr is 1
|
|
|
99 nticks = 10
|
|
|
100 last = max(d$BP)
|
|
|
101 first = min(d$BP)
|
|
|
102 tks = c(first)
|
|
|
103 t = (last-first)/nticks # units per tick
|
|
|
104 for (x in c(1:(nticks))) {
|
|
|
105 tks = c(tks,round(x*t)+first) }
|
|
|
106 ticklim = c(first,last)
|
|
|
107 } # else
|
|
|
108 if (grey) {mycols=rep(c("gray10","gray60"),max(d$CHR))
|
|
|
109 } else {
|
|
|
110 mycols=rep(coloursTouse,max(d$CHR))
|
|
|
111 }
|
|
|
112 dlen = length(d$P)
|
|
|
113 d$pranks = rank(d$P)/dlen
|
|
|
114 d$centiles = 100*d$pranks # small are interesting
|
|
|
115 d$sizes = ifelse((d$centile < 1),2,1)
|
|
|
116 if (annotate) d.annotate=d[as.numeric(substr(d$SNP,3,100)) %in% SNPlist, ]
|
|
|
117 if (nchr >= 2) {
|
|
|
118 manplot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR),size=factor(sizes))
|
|
|
119 manplot=manplot+scale_x_continuous(name="Chromosome", breaks=tks, labels=xlabs,limits=ticklim)
|
|
|
120 manplot=manplot+scale_size_manual(values = c(0.5,1.5)) # requires discreet scale - eg factor
|
|
|
121 #manplot=manplot+scale_size(values=c(0.5,2)) # requires continuous
|
|
|
122 }
|
|
|
123 else {
|
|
|
124 manplot=qplot(BP,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
|
|
|
125 manplot=manplot+scale_x_continuous(name=paste("Chromosome",chrlist[1]), breaks=tks, labels=tks,limits=ticklim)
|
|
|
126 }
|
|
|
127 manplot=manplot+scale_y_continuous(limits=c(0,maxy), breaks=1:maxy, labels=1:maxy)
|
|
|
128 manplot=manplot+scale_colour_manual(value=mycols)
|
|
|
129 if (annotate) { manplot=manplot + geom_point(data=d.annotate, colour=I("green3")) }
|
|
|
130 manplot=manplot + opts(legend.position = "none")
|
|
|
131 manplot=manplot + opts(title=title)
|
|
|
132 manplot=manplot+opts(
|
|
|
133 panel.background=theme_blank(),
|
|
|
134 axis.text.x=theme_text(size=size.x.labels, colour="grey50"),
|
|
|
135 axis.text.y=theme_text(size=size.y.labels, colour="grey50"),
|
|
|
136 axis.ticks=theme_segment(colour=NA)
|
|
|
137 )
|
|
|
138 if (suggestiveline) manplot=manplot+geom_hline(yintercept=suggestiveline,colour="blue", alpha=I(1/3))
|
|
|
139 if (genomewideline) manplot=manplot+geom_hline(yintercept=genomewideline,colour="red")
|
|
|
140 manplot
|
|
|
141 } else {
|
|
|
142 stop("Make sure your data frame contains columns CHR, BP, and P")
|
|
|
143 }
|
|
|
144 }
|
|
|
145
|
|
|
146
|
|
|
147
|
|
|
148 qq = function(pvector, title=NULL, spartan=F) {
|
|
|
149 # Thanks to Daniel Shriner at NHGRI for providing this code for creating expected and observed values
|
|
|
150 o = -log10(sort(pvector,decreasing=F))
|
|
|
151 e = -log10( 1:length(o)/length(o) )
|
|
|
152 # you could use base graphics
|
|
|
153 # plot(e,o,pch=19,cex=0.25, xlab=expression(Expected~~-log[10](italic(p))),
|
|
|
154 # ylab=expression(Observed~~-log[10](italic(p))), xlim=c(0,max(e)), ylim=c(0,max(e)))
|
|
|
155 # lines(e,e,col="red")
|
|
|
156 #You'll need ggplot2 installed to do the rest
|
|
|
157 qq=qplot(e,o, xlim=c(0,max(e)), ylim=c(0,max(o))) + stat_abline(intercept=0,slope=1, col="red")
|
|
|
158 qq=qq+opts(title=title)
|
|
|
159 qq=qq+scale_x_continuous(name=expression(Expected~~-log[10](italic(p))))
|
|
|
160 qq=qq+scale_y_continuous(name=expression(Observed~~-log[10](italic(p))))
|
|
|
161 if (spartan) plot=plot+opts(panel.background=theme_rect(col="grey50"), panel.grid.minor=theme_blank())
|
|
|
162 qq
|
|
|
163 }
|
|
|
164
|
|
|
165 """
|
|
|
166
|
|
|
167 # we need another string to avoid confusion over string substitutions with %in%
|
|
|
168 # instantiate rcode2 string with infile,chromcol,offsetcol,pvalscols,title before saving and running
|
|
|
169
|
|
|
170 rcode2 = """rgqqMan = function(infile="%s",chromcolumn=%d, offsetcolumn=%d, pvalscolumns=c(%s),
|
|
|
171 title="%s",grey=%d) {
|
|
|
172 rawd = read.table(infile,head=T,sep='\\t')
|
|
|
173 dn = names(rawd)
|
|
|
174 cc = dn[chromcolumn]
|
|
|
175 oc = dn[offsetcolumn]
|
|
|
176 rawd[,cc] = sub('chr','',rawd[,cc],ignore.case = T) # just in case
|
|
|
177 rawd[,cc] = sub(':','',rawd[,cc],ignore.case = T) # ugh
|
|
|
178 rawd[,cc] = sub('X',23,rawd[,cc],ignore.case = T)
|
|
|
179 rawd[,cc] = sub('Y',24,rawd[,cc],ignore.case = T)
|
|
|
180 rawd[,cc] = sub('Mt',25,rawd[,cc], ignore.case = T)
|
|
|
181 nams = c(cc,oc) # for sorting
|
|
|
182 plen = length(rawd[,1])
|
|
|
183 print(paste('###',plen,'values read from',infile,'read - now running plots',sep=' '))
|
|
|
184 rawd = rawd[do.call(order,rawd[nams]),]
|
|
|
185 # mmmf - suggested by http://onertipaday.blogspot.com/2007/08/sortingordering-dataframe-according.html
|
|
|
186 # in case not yet ordered
|
|
|
187 if (plen > 0) {
|
|
|
188 for (pvalscolumn in pvalscolumns) {
|
|
|
189 if (pvalscolumn > 0)
|
|
|
190 {
|
|
|
191 cname = names(rawd)[pvalscolumn]
|
|
|
192 mytitle = paste('p=',cname,', ',title,sep='')
|
|
|
193 myfname = chartr(' ','_',cname)
|
|
|
194 myqqplot = qq(rawd[,pvalscolumn],title=mytitle)
|
|
|
195 ggsave(filename=paste(myfname,"qqplot.png",sep='_'),myqqplot,width=8,height=6,dpi=96)
|
|
|
196 ggsave(filename=paste(myfname,"qqplot.pdf",sep='_'),myqqplot,width=8,height=6,dpi=96)
|
|
|
197 print(paste('## qqplot on',cname,'done'))
|
|
|
198 if ((chromcolumn > 0) & (offsetcolumn > 0)) {
|
|
|
199 print(paste('## manhattan on',cname,'starting',chromcolumn,offsetcolumn,pvalscolumn))
|
|
|
200 mymanplot= DrawManhattan(chrom=rawd[,chromcolumn],offset=rawd[,offsetcolumn],pvals=rawd[,pvalscolumn],title=mytitle,grey=grey)
|
|
|
201 ggsave(filename=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=8,height=6,dpi=96)
|
|
|
202 ggsave(filename=paste(myfname,"manhattan.pdf",sep='_'),mymanplot,width=8,height=6,dpi=96)
|
|
|
203 print(paste('## manhattan plot on',cname,'done'))
|
|
|
204 }
|
|
|
205 else {
|
|
|
206 print(paste('chrom column =',chromcolumn,'offset column = ',offsetcolumn,
|
|
|
207 'so no Manhattan plot - supply both chromosome and offset as numerics for Manhattan plots if required'))
|
|
|
208 }
|
|
|
209 }
|
|
|
210 else {
|
|
|
211 print(paste('pvalue column =',pvalscolumn,'Cannot parse it so no plots possible'))
|
|
|
212 }
|
|
|
213 } # for pvalscolumn
|
|
|
214 } else { print('## Problem - no values available to plot - was there really a chromosome and offset column?') }
|
|
|
215 }
|
|
|
216
|
|
|
217 rgqqMan()
|
|
|
218 # execute with defaults as substituted
|
|
|
219 """
|
|
|
220
|
|
|
221
|
|
|
222 def doManQQ(input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir,beTidy=False):
|
|
|
223 """
|
|
|
224 we may have an interval file or a tabular file - if interval, will have chr1... so need to adjust
|
|
|
225 to chrom numbers
|
|
|
226 draw a qq for pvals and a manhattan plot if chrom/offset <> 0
|
|
|
227 contains some R scripts as text strings - we substitute defaults into the calls
|
|
|
228 to make them do our bidding - and save the resulting code for posterity
|
|
|
229 this can be called externally, I guess...for QC eg?
|
|
|
230 """
|
|
|
231 if debug:
|
|
|
232 print 'doManQQ',input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir
|
|
|
233 rcmd = '%s%s' % (rcode,rcode2 % (input_fname,chrom_col,offset_col,pval_cols,title,grey))
|
|
|
234 if debug:
|
|
|
235 print 'running\n%s\n' % rcmd
|
|
|
236 rlog,flist = RRun(rcmd=rcmd,title=ctitle,outdir=outdir)
|
|
|
237 rlog.append('## R script=')
|
|
|
238 rlog.append(rcmd)
|
|
|
239 return rlog,flist
|
|
|
240
|
|
|
241 def compressPDF(inpdf=None):
|
|
|
242 """need absolute path to pdf
|
|
|
243 """
|
|
|
244 assert os.path.isfile(inpdf), "## Input %s supplied to compressPDF not found" % inpdf
|
|
|
245 outpdf = '%s_compressed' % inpdf
|
|
|
246 cl = ["gs", "-sDEVICE=pdfwrite", "-dNOPAUSE", "-dBATCH", "-sOutputFile=%s" % outpdf,inpdf]
|
|
|
247 retval = subprocess.call(cl)
|
|
|
248 if retval == 0:
|
|
|
249 os.unlink(inpdf)
|
|
|
250 shutil.move(outpdf,inpdf)
|
|
|
251 return retval
|
|
|
252
|
|
|
253 def main():
|
|
|
254 u = """<command interpreter="python">
|
|
|
255 rgManQQ.py '$input_file' "$name" '$out_html' '$out_html.files_path' '$chrom_col' '$offset_col' '$pval_col'
|
|
|
256 </command>
|
|
|
257 """
|
|
|
258 npar = 8
|
|
|
259 if len(sys.argv) < npar:
|
|
|
260 print >> sys.stdout, '## error - too few command line parameters - wanting %d' % npar
|
|
|
261 print >> sys.stdout, u
|
|
|
262 sys.exit(1)
|
|
|
263 input_fname = sys.argv[1]
|
|
|
264 title = sys.argv[2]
|
|
|
265 killme = string.punctuation + string.whitespace
|
|
|
266 trantab = string.maketrans(killme,'_'*len(killme))
|
|
|
267 ctitle = title.translate(trantab)
|
|
|
268 outhtml = sys.argv[3]
|
|
|
269 outdir = sys.argv[4]
|
|
|
270 try:
|
|
|
271 chrom_col = int(sys.argv[5])
|
|
|
272 except:
|
|
|
273 chrom_col = -1
|
|
|
274 try:
|
|
|
275 offset_col = int(sys.argv[6])
|
|
|
276 except:
|
|
|
277 offset_col = -1
|
|
|
278 p = sys.argv[7].strip().split(',')
|
|
|
279 try:
|
|
|
280 q = [int(x) for x in p]
|
|
|
281 except:
|
|
|
282 p = -1
|
|
|
283 if chrom_col == -1 or offset_col == -1: # was passed as zero - do not do manhattan plots
|
|
|
284 chrom_col = -1
|
|
|
285 offset_col = -1
|
|
|
286 grey = 0
|
|
|
287 if (sys.argv[8].lower() in ['1','true']):
|
|
|
288 grey = 1
|
|
|
289 if p == -1:
|
|
|
290 print >> sys.stderr,'## Cannot run rgManQQ - missing pval column'
|
|
|
291 sys.exit(1)
|
|
|
292 p = ['%d' % (int(x) + 1) for x in p]
|
|
|
293 rlog,flist = doManQQ(input_fname,chrom_col+1,offset_col+1,','.join(p),title,grey,ctitle,outdir)
|
|
|
294 flist.sort()
|
|
|
295 html = [galhtmlprefix % progname,]
|
|
|
296 html.append('<h1>%s</h1>' % title)
|
|
|
297 if len(flist) > 0:
|
|
|
298 html.append('<table>\n')
|
|
|
299 for row in flist:
|
|
|
300 fname,expl = row # RRun returns pairs of filenames fiddled for the log and R script
|
|
|
301 n,e = os.path.splitext(fname)
|
|
|
302 if e in ['.png','.jpg']:
|
|
|
303 pdf = '%s.pdf' % n
|
|
|
304 pdff = os.path.join(outdir,pdf)
|
|
|
305 if os.path.exists(pdff):
|
|
|
306 rval = compressPDF(inpdf=pdff)
|
|
|
307 if rval <> 0:
|
|
|
308 pdf = '%s(not_compressed)' % pdf
|
|
|
309 else:
|
|
|
310 pdf = '%s(not_found)' % pdf
|
|
|
311 s= '<tr><td><a href="%s"><img src="%s" title="%s" hspace="10" width="800"></a></td></tr>' \
|
|
|
312 % (pdf,fname,expl)
|
|
|
313 html.append(s)
|
|
|
314 else:
|
|
|
315 html.append('<tr><td><a href="%s">%s</a></td></tr>' % (fname,expl))
|
|
|
316 html.append('</table>\n')
|
|
|
317 else:
|
|
|
318 html.append('<h2>### Error - R returned no files - please confirm that parameters are sane</h1>')
|
|
|
319 html.append('<h3>R log follows below</h3><hr><pre>\n')
|
|
|
320 html += rlog
|
|
|
321 html.append('</pre>\n')
|
|
|
322 html.append(galhtmlattr % (progname,timenow()))
|
|
|
323 html.append(galhtmlpostfix)
|
|
|
324 htmlf = file(outhtml,'w')
|
|
|
325 htmlf.write('\n'.join(html))
|
|
|
326 htmlf.write('\n')
|
|
|
327 htmlf.close()
|
|
|
328
|
|
|
329
|
|
|
330
|
|
|
331 if __name__ == "__main__":
|
|
|
332 main()
|
|
|
333
|
|
|
334
|
