opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-10-16 07:19:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fixes to rgManQQ.py to handle chromosome oddities like 'chr1' and X.

Browse Source 
Also cleaned up redundant file recreation - the code passes column numbers to R so the raw data can be read as supplied by Galaxy
Thanks to Alison Harrill for exposing some more input edge cases...
This commit is contained in:
Ross Lazarus
2011-07-20 12:17:43 +10:00
parent 7e4bcb268e
commit 7d2dd59186
2 changed files with 38 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/rgenetics/rgManQQ.py
+37 -66
View File
@@ -1,5 +1,7 @@
#!/usr/local/bin/python
# rgmanqq updated july 19 to deal with x,y and mt
# lots of fixes
# ross lazarus
import sys,math,shutil,subprocess,os,time,tempfile,string
from os.path import abspath
from rgutils import timenow, RRun, galhtmlprefix, galhtmlpostfix, galhtmlattr
@@ -18,7 +20,7 @@ rcode="""
# http://StephenTurner.us/
# http://GettingGeneticsDone.blogspot.com/
# Last updated: Tuesday, December 22, 2009
# Last updated: 19 July 2011 by Ross Lazarus
# R code for making manhattan plots and QQ plots from plink output files.
# With GWAS data this can take a lot of memory. Recommended for use on
# 64bit machines only, for now.
@@ -28,26 +30,30 @@ rcode="""
library(ggplot2)
coloursTouse = c('firebrick','darkblue','goldenrod','darkgreen')
# not too fugly but need a colour expert please...
# not too ugly but need a colour expert please...
manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
suggestiveline=0, genomewide=T, size.x.labels=9, size.y.labels=10, annotate=F, SNPlist=NULL,grey=0) {
DrawManhattan = function(pvals=Null,chrom=Null,offset=Null,title=NULL, max.y="max",suggestiveline=0, genomewide=T, size.x.labels=9,
size.y.labels=10, annotate=F, SNPlist=NULL,grey=0) {
if (annotate & is.null(SNPlist)) stop("You requested annotation but provided no SNPlist!")
genomewideline=NULL # was genomewideline=-log10(5e-8)
if (genomewide) { # use bonferroni since might be only a small region?
genomewideline = -log10(0.05/length(pvals)) }
d=data.frame(CHR=chrom,BP=offset,P=pvals)
#limit to only chrs 1-23?
d=d[d$CHR %in% 1:23, ]
chro = sub('chr','',chrom, ignore.case = T) # just in case
chro = sub(':','',chro, ignore.case = T) # ugh
chro = sub('X',23,chro, ignore.case = T)
chro = sub('Y',24,chro, ignore.case = T)
chro = sub('Mt',25,chro, ignore.case = T)
offset = as.integer(offset)
pvals = as.double(pvals)
chro = as.integer(chro)
d=data.frame(CHR=chro,BP=offset,P=pvals)
#limit to only chrs 1-22, x=23,y=24,Mt=25?
d=d[d$CHR %in% 1:25, ]
if ("CHR" %in% names(d) & "BP" %in% names(d) & "P" %in% names(d) ) {
d=na.omit(d)
#d=na.omit(d)
d=d[d$P>0 & d$P<=1, ]
d$logp = -log10(d$P)
d$logp = as.double(-log10(d$P))
d$pos=NA
ticks=NULL
lastbase=0
@@ -63,7 +69,11 @@ manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
lastchr = chrlist[x-1] # previous whatever the list
lastbase=lastbase+tail(subset(d,CHR==lastchr)$BP, 1)
d[d$CHR==i, ]$pos=d[d$CHR==i, ]$BP+lastbase
if (sum(is.na(lastchr),is.na(lastbase),is.na(d[d$CHR==i, ]$pos))) {
cat(paste('manhattan: For',title,'chrlistx=',i,'lastchr=',lastchr,'lastbase=',lastbase,'pos=',d[d$CHR==i,]$pos))
}
tks=c(tks, d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1])
}
ticklim=c(min(d$pos),max(d$pos))
xlabs = chrlist
@@ -85,8 +95,6 @@ manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
if (max.y=="max") maxy=ceiling(max(d$logp)) else maxy=max.y
maxy = max(maxy,1.1*genomewideline)
# if (maxy<8) maxy=8
# only makes sense if genome wide is assumed - we could have a fine mapping region?
if (annotate) d.annotate=d[as.numeric(substr(d$SNP,3,100)) %in% SNPlist, ]
if (nchr >= 2) {
manplot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
@@ -105,9 +113,6 @@ manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
axis.text.y=theme_text(size=size.y.labels, colour="grey50"),
axis.ticks=theme_segment(colour=NA)
)
#manplot = manplot + opts(panel.grid.y.minor=theme_blank(),panel.grid.y.major=theme_blank())
#manplot = manplot + opts(panel.grid.major=theme_blank())
if (suggestiveline) manplot=manplot+geom_hline(yintercept=suggestiveline,colour="blue", alpha=I(1/3))
if (genomewideline) manplot=manplot+geom_hline(yintercept=genomewideline,colour="red")
manplot
@@ -134,12 +139,13 @@ qq = function(pvector, title=NULL, spartan=F) {
if (spartan) plot=plot+opts(panel.background=theme_rect(col="grey50"), panel.grid.minor=theme_blank())
qq
}
"""
# we need another string to avoid confusion over string substitutions with %in%
# instantiate rcode2 string with infile,chromcol,offsetcol,pvalscols,title before saving and running
rcode2 = """rgqqMan = function(infile="%s",chromcolumn=%s, offsetcolumn=%s, pvalscolumns=%s,
rcode2 = """rgqqMan = function(infile="%s",chromcolumn=%d, offsetcolumn=%d, pvalscolumns=c(%s),
title="%s",grey=%d) {
rawd = read.table(infile,head=T,sep='\\t')
dn = names(rawd)
@@ -167,7 +173,7 @@ if (plen > 0) {
doreorder = 0
}
print(paste('## manhattan on',cname,'starting',chromcolumn,offsetcolumn,pvalscolumn))
mymanplot= manhattan(chrom=rawd[,chromcolumn],offset=rawd[,offsetcolumn],pvals=rawd[,pvalscolumn],title=mytitle,grey=grey)
mymanplot= DrawManhattan(chrom=rawd[,chromcolumn],offset=rawd[,offsetcolumn],pvals=rawd[,pvalscolumn],title=mytitle,grey=grey)
print(paste('## manhattan plot on',cname,'done'))
ggsave(filename=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=6,height=4,dpi=100)
}
@@ -198,50 +204,13 @@ def doManQQ(input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir,
this can be called externally, I guess...for QC eg?
"""
if debug:
print 'doManQQ',input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir
ffd,filtered_fname = tempfile.mkstemp(prefix='rgManQQtemp')
f = open(filtered_fname,'w')
inf = open(input_fname,'r')
ohead = inf.readline().strip().split('\t') # see if we have a header
inf.seek(0) # rewind
newhead = ['pval%d' % (x+1) for x in pval_cols]
newhead.insert(0,'Offset')
newhead.insert(0,'Chrom')
havehead = 0
wewant = [chrom_col,offset_col]
wewant += pval_cols
try:
allnums = ['%d' % x for x in ohead] # this should barf if non numerics == header row?
f.write('\t'.join(newhead)) # for R to read
f.write('\n')
except:
havehead = 1
newhead = [ohead[chrom_col],ohead[offset_col]]
newhead += [ohead[x] for x in pval_cols]
f.write('\t'.join(newhead)) # use the original head
f.write('\n')
for i,row in enumerate(inf):
if i == 0 and havehead:
continue # ignore header
sr = row.strip().split('\t')
if len(sr) > 1:
if sr[chrom_col].lower().find('chr') <> -1:
sr[chrom_col] = sr[chrom_col][3:]
newr = [sr[x] for x in wewant] # grab cols we need
s = '\t'.join(newr)
f.write(s)
f.write('\n')
f.close()
pvc = [x+3 for x in range(len(pval_cols))] # 2 for offset and chrom, 1 for r offset start
pvc = 'c(%s)' % (','.join(map(str,pvc)))
rcmd = '%s%s' % (rcode,rcode2 % (filtered_fname,'1','2',pvc,title,grey))
print 'doManQQ',input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir
rcmd = '%s%s' % (rcode,rcode2 % (input_fname,chrom_col,offset_col,pval_cols,title,grey))
if debug:
print 'running\n%s\n' % rcmd
print 'running\n%s\n' % rcmd
rlog,flist = RRun(rcmd=rcmd,title=ctitle,outdir=outdir)
rlog.append('## R script=')
rlog.append(rcmd)
if beTidy:
os.unlink(filtered_fname)
return rlog,flist
@@ -272,19 +241,20 @@ def main():
offset_col = -1
p = sys.argv[7].strip().split(',')
try:
p = [int(x) for x in p]
q = [int(x) for x in p]
except:
p = [-1]
p = -1
if chrom_col == -1 or offset_col == -1: # was passed as zero - do not do manhattan plots
chrom_col = -1
offset_col = -1
grey = 0
if (sys.argv[8].lower() in ['1','true']):
grey = 1
if p == [-1]:
if p == -1:
print >> sys.stderr,'## Cannot run rgManQQ - missing pval column'
sys.exit(1)
rlog,flist = doManQQ(input_fname,chrom_col,offset_col,p,title,grey,ctitle,outdir)
p = ['%d' % (int(x) + 1) for x in p]
rlog,flist = doManQQ(input_fname,chrom_col+1,offset_col+1,','.join(p),title,grey,ctitle,outdir)
flist.sort()
html = [galhtmlprefix % progname,]
html.append('<h1>%s</h1>' % title)
@@ -294,7 +264,7 @@ def main():
fname,expl = row # RRun returns pairs of filenames fiddled for the log and R script
e = os.path.splitext(fname)[-1]
if e in ['.png','.jpg']:
s= '<tr><td><a href="%s"><img src="%s" alt="%s hspace="10" width="400"><br>(Click to download image %s)</a></td></tr>' \
s= '<tr><td><a href="%s"><img src="%s" title="%s hspace="10" width="400"><br>(Click to download image %s)</a></td></tr>' \
% (fname,fname,expl,expl )
html.append(s)
else:
@@ -317,3 +287,4 @@ def main():
if __name__ == "__main__":
main()
tools/rgenetics/rgManQQ.xml
+1 -1
View File
@@ -1,4 +1,4 @@
<tool id="rgManQQ1" name="Manhattan/QQ:" version="1.0.1">
<tool id="rgManQQ1" name="Manhattan/QQ:" version="1.0.2">
<code file="rgManQQ_code.py"/>
<description>Plots for WGA P values</description>