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

Fixes to some rgenetics tests and test outputs

Browse Source 
Some odd bugs in the new version of ggplot2 in R
Fixes to the Rrun function in rgutils
This commit is contained in:
Ross Lazarus
2010-08-31 02:27:33 -04:00
parent 2ba1bd61bc
commit 1fe808274d
5 changed files with 41 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/rgenetics/rgGRR.py
+1 -1
View File
@@ -78,7 +78,7 @@ class NullDevice:
pass
tempstderr = sys.stderr # save
sys.stderr = NullDevice()
#sys.stderr = NullDevice()
# need to avoid blather about deprecation and other strange stuff from scipy
# the current galaxy job runner assumes that
# the job is in error if anything appears on sys.stderr
tools/rgenetics/rgManQQ.py
+22 -21
View File
@@ -75,25 +75,25 @@ manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
# 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, ]
plot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
plot=plot+scale_x_continuous(name="Chromosome", breaks=ticks, labels=(unique(d$CHR)))
plot=plot+scale_y_continuous(limits=c(0,maxy), breaks=1:maxy, labels=1:maxy)
plot=plot+scale_colour_manual(value=mycols)
if (annotate) plot=plot + geom_point(data=d.annotate, colour=I("green3"))
plot=plot + opts(legend.position = "none")
plot=plot + opts(title=title)
plot=plot+opts(
manplot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
manplot=manplot+scale_x_continuous(name="Chromosome", breaks=ticks, labels=(unique(d$CHR)))
manplot=manplot+scale_y_continuous(limits=c(0,maxy), breaks=1:maxy, labels=1:maxy)
manplot=manplot+scale_colour_manual(value=mycols)
if (annotate) { manplot=manplot + geom_point(data=d.annotate, colour=I("green3")) }
manplot=manplot + opts(legend.position = "none")
manplot=manplot + opts(title=title)
manplot=manplot+opts(
panel.background=theme_blank(),
axis.text.x=theme_text(size=size.x.labels, colour="grey50"),
axis.text.y=theme_text(size=size.y.labels, colour="grey50"),
axis.ticks=theme_segment(colour=NA)
)
#plot = plot + opts(panel.grid.y.minor=theme_blank(),panel.grid.y.major=theme_blank())
#plot = plot + opts(panel.grid.major=theme_blank())
#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) plot=plot+geom_hline(yintercept=suggestiveline,colour="blue", alpha=I(1/3))
if (genomewideline) plot=plot+geom_hline(yintercept=genomewideline,colour="red")
plot
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
} else {
stop("Make sure your data frame contains columns CHR, BP, and P")
}
@@ -110,12 +110,12 @@ qq = function(pvector, title=NULL, spartan=F) {
# ylab=expression(Observed~~-log[10](italic(p))), xlim=c(0,max(e)), ylim=c(0,max(e)))
# lines(e,e,col="red")
#You'll need ggplot2 installed to do the rest
plot=qplot(e,o, xlim=c(0,max(e)), ylim=c(0,max(o))) + stat_abline(intercept=0,slope=1, col="red")
plot=plot+opts(title=title)
plot=plot+scale_x_continuous(name=expression(Expected~~-log[10](italic(p))))
plot=plot+scale_y_continuous(name=expression(Observed~~-log[10](italic(p))))
qq=qplot(e,o, xlim=c(0,max(e)), ylim=c(0,max(o))) + stat_abline(intercept=0,slope=1, col="red")
qq=qq+opts(title=title)
qq=qq+scale_x_continuous(name=expression(Expected~~-log[10](italic(p))))
qq=qq+scale_y_continuous(name=expression(Observed~~-log[10](italic(p))))
if (spartan) plot=plot+opts(panel.background=theme_rect(col="grey50"), panel.grid.minor=theme_blank())
plot
qq
}
"""
@@ -140,6 +140,7 @@ if (plen > 0) {
mytitle = paste('p=',cname,', ',title,sep='')
myfname = chartr(' ','_',cname)
myqqplot = qq(rawd[,pvalscolumn],title=mytitle)
ggsave(filename=paste(myfname,"qqplot.png",sep='_'),myqqplot,width=8,height=11)
print(paste('## qqplot on',cname,'done'))
if ((chromcolumn > 0) & (offsetcolumn > 0)) {
if (doreorder) {
@@ -151,13 +152,12 @@ if (plen > 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)
print(paste('## manhattan plot on',cname,'done'))
ggsave(file=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=11,height=8,dpi=100)
ggsave(filename=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=11,height=8)
}
else {
print(paste('chrom column =',chromcolumn,'offset column = ',offsetcolumn,
'so no Manhattan plot - supply both chromosome and offset as numerics for Manhattan plots if required'))
}
ggsave(file=paste(myfname,"qqplot.png",sep='_'),myqqplot,width=8,height=11,dpi=100)
}
else {
print(paste('pvalue column =',pvalscolumn,'Cannot parse it so no plots possible'))
@@ -177,9 +177,10 @@ def doManQQ(input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir)
to make them do our bidding - and save the resulting code for posterity
this can be called externally, I guess...for QC eg?
"""
rcmd = '%s%s' % (rcode,rcode2 % (input_fname,chrom_col,offset_col,pval_cols,title,grey))
rlog,flist = RRun(rcmd=rcmd,title=ctitle,outdir=outdir)
rlog.append('## R script=')
rlog.append(rcmd)
return rlog,flist
tools/rgenetics/rgManQQ.xml
+3 -7
View File
@@ -42,14 +42,10 @@
<param name='chrom_col' value='1' />
<param name='offset_col' value='2' />
<param name='grey' value='0' />
<output name='out_html' file='rgtestouts/rgManQQ/rgManQQtest1.html' ftype='html' lines_diff='17'>
<extra_files type="file" name='Allelep_manhattan.png' value='rgtestouts/rgManQQ/Allelep_manhattan.png' compare="sim_size"
delta = "4000" />
<extra_files type="file" name='Allelep_qqplot.png' value='rgtestouts/rgManQQ/Allelep_qqplot.png' compare="sim_size"
delta = "4000" />
<extra_files type="file" name='Armitagep_manhattan.png' value='rgtestouts/rgManQQ/Armitagep_manhattan.png' compare="sim_size"
<output name='out_html' file='rgtestouts/rgManQQ/rgManQQtest1.html' ftype='html' lines_diff='30'>
<extra_files type="file" name='Allelep_manhattan.png' value='rgtestouts/rgManQQ/Allelep_manhattan.png' compare="sim_size"
delta = "4000"/>
<extra_files type="file" name='Armitagep_qqplot.png' value='rgtestouts/rgManQQ/Armitagep_qqplot.png' compare="sim_size"
<extra_files type="file" name='Allelep_qqplot.png' value='rgtestouts/rgManQQ/Allelep_qqplot.png' compare="sim_size"
delta = "4000" />
<extra_files type="file" name='rgManQQtest1.R' value='rgtestouts/rgManQQ/rgManQQtest1.R' compare="diff" lines_diff="160"/>
</output>
tools/rgenetics/rgtest.sh
+1 -1
View File
@@ -41,7 +41,7 @@ CL="python $TOOLPATH/$TOOL.py "$INPATH/smallwgaP.xls" $NPRE ${OUTPATH}/${NPRE}.h
#python /opt/galaxy/tools/rgenetics/rgManQQ.py /opt/galaxy/test-data/smallwgaP.xls rgManQQtest1
#/opt/galaxy/test-data/rgtestouts/rgManQQ/rgManQQtest1.html /opt/galaxy/test-data/rgtestouts/rgManQQ 1 2 5,7
echo "Testing $TOOL using $CL"
python $TOOLPATH/$TOOL.py "$INPATH/smallwgaP.xls" $NPRE ${OUTPATH}/${NPRE}.html $OUTPATH 1 2 5,7 0
python $TOOLPATH/$TOOL.py "$INPATH/smallwgaP.xls" $NPRE ${OUTPATH}/${NPRE}.html $OUTPATH 1 2 7 0
TOOL="rgfakePhe"
NPRE=${TOOL}test1
OUTPATH="$OROOT/$TOOL"
tools/rgenetics/rgutils.py
+14 -18
View File
@@ -122,32 +122,36 @@ def RRun(rcmd=[],outdir=None,title='myR',tidy=True):
tempout=False
rname = '%s.R' % title
stoname = '%s.R.log' % title
cwd = os.getcwd()
if outdir: # want a specific path
try:
os.makedirs(outdir) # might not be there yet...
except:
pass
else:
tempout = True
outdir = tempfile.mkdtemp(prefix=title)
os.chdir(outdir)
if type(rcmd) == type([]):
script = '\n'.join(rcmd)
else: # string
script = rcmd
sto = file(stoname,'w')
rscriptname = os.path.join(outdir,rname)
rscript = file(rscriptname,'w')
rscript = file(rname,'w')
rscript.write(script)
rscript.write('\n#R script autogenerated by rgenetics/rgutils.py on %s\n' % timenow())
rscript.close()
vcl = '%s --slave --vanilla < %s' % (rexe,rscriptname)
x = subprocess.Popen(vcl,shell=True,stderr=sto,stdout=sto,cwd=outdir)
vcl = '%s --slave --vanilla < %s' % (rexe,rname)
if outdir:
x = subprocess.Popen(vcl,shell=True,stderr=sto,stdout=sto,cwd=outdir)
else:
x = subprocess.Popen(vcl,shell=True,stderr=sto,stdout=sto)
retval = x.wait()
sto.close()
rlog = file(stoname,'r').readlines()
if retval <> 0:
rlog.insert(0,'Nonzero exit code = %d' % retval) # indicate failure
flist = os.listdir(outdir)
if outdir:
flist = os.listdir(outdir)
else:
flist = os.listdir(os.getcwd())
flist.sort
flist = [(x,x) for x in flist]
for i,x in enumerate(flist):
@@ -155,16 +159,8 @@ def RRun(rcmd=[],outdir=None,title='myR',tidy=True):
flist[i] = (x,'R script for %s' % title)
elif x == stoname:
flist[i] = (x,'R log for %s' % title)
if tidy and tempout: # we used a temp one - fix that
flist = os.listdir(outdir)
for f in flist:
if os.path.isdir(f):
flist2 = os.listdir(f)
for f2 in flist2:
os.unlink(os.path.abspath(f2))
else:
os.unlink(os.path.abspath(f))
os.rmdir(outdir)
if outdir:
os.chdir(cwd)
return rlog,flist # for html layout
def runPlink(bfn='bar',ofn='foo',logf=None,plinktasks=[],cd='./',vclbase = []):