Mercurial > repos > iuc > limma_voom
changeset 7:e6a4ff41af6b draft
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit cf399638ebca4250bcc15f468238a9964de97b33
| author | iuc |
|---|---|
| date | Wed, 09 May 2018 13:27:14 -0400 |
| parents | 39fa12a6d885 |
| children | 00a42b66e522 |
| files | limma_voom.R limma_voom.xml test-data/out_rscript.txt |
| diffstat | 3 files changed, 106 insertions(+), 24 deletions(-) [+] |
line wrap: on
line diff
--- a/limma_voom.R Tue May 08 18:12:40 2018 -0400 +++ b/limma_voom.R Wed May 09 13:27:14 2018 -0400 @@ -22,13 +22,17 @@ # robOpt", "b", 0, "logical" -String specifying if robust options should be used # trend", "t", 1, "double" -Float for prior.count if limma-trend is used instead of voom # weightOpt", "w", 0, "logical" -String specifying if voomWithQualityWeights should be used -# volhiOpt", "v", 1, "integer" -Integer specifying number of genes to highlight in volcano plot +# topgenes", "G", 1, "integer" -Integer specifying no. of genes to highlight in volcano and heatmap # treatOpt", "T", 0, "logical" -String specifying if TREAT function shuld be used # # OUT: +# Density Plots (if filtering) +# Box Plots (if normalising) # MDS Plot # Voom/SA plot # MD Plot +# Volcano Plot +# Heatmap # Expression Table # HTML file linking to the ouputs # Optional: @@ -37,7 +41,7 @@ # # # Author: Shian Su - registertonysu@gmail.com - Jan 2014 -# Modified by: Maria Doyle - Jun 2017, Jan 2018 +# Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018 # Record starting time timeStart <- as.character(Sys.time()) @@ -50,6 +54,7 @@ library(limma, quietly=TRUE, warn.conflicts=FALSE) library(scales, quietly=TRUE, warn.conflicts=FALSE) library(getopt, quietly=TRUE, warn.conflicts=FALSE) +library(gplots, quietly=TRUE, warn.conflicts=FALSE) ################################################################################ ### Function Declaration @@ -169,7 +174,7 @@ "robOpt", "b", 0, "logical", "trend", "t", 1, "double", "weightOpt", "w", 0, "logical", - "volhiOpt", "v", 1, "integer", + "topgenes", "G", 1, "integer", "treatOpt", "T", 0, "logical"), byrow=TRUE, ncol=4) opt <- getopt(spec) @@ -330,11 +335,13 @@ mdOutPdf <- character() volOutPdf <- character() +heatOutPdf <- character() mdvolOutPng <- character() topOut <- character() for (i in 1:length(contrastData)) { mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf")) volOutPdf[i] <- makeOut(paste0("volplot_", contrastData[i], ".pdf")) + heatOutPdf[i] <- makeOut(paste0("heatmap_", contrastData[i], ".pdf")) mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", contrastData[i], ".png")) topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv")) } @@ -691,12 +698,12 @@ # labels must be in second column currently limma::volcanoplot(fit, coef=i, main=paste("Volcano Plot:", unmake.names(contrastData[i])), - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes[, 2]) } else { limma::volcanoplot(fit, coef=i, - main=paste("Volcano Plot:", unmake.names(contrastData[i])),, - highlight=opt$volhiOpt, + main=paste("Volcano Plot:", unmake.names(contrastData[i])), + highlight=opt$topgenes, names=fit$genes$GeneID) } linkName <- paste0("VolcanoPlot_", contrastData[i], ".pdf") @@ -704,24 +711,53 @@ linkData <- rbind(linkData, c(linkName, linkAddr)) invisible(dev.off()) + # Plot Heatmap + topgenes <- rownames(top[1:opt$topgenes, ]) + if (wantTrend) { + topexp <- plotData[topgenes, ] + } else { + topexp <- plotData$E[topgenes, ] + } + + pdf(heatOutPdf[i]) + par(cex.main=0.8) + if (haveAnno) { + # labels must be in second column currently + heatmap.2(topexp, scale="row", + main=paste("Heatmap:", unmake.names(contrastData[i])), + col="bluered", trace="none", density.info="none", + margin=c(8,6), lhei=c(2,10), dendrogram="column", + cexRow=0.7, cexCol=0.7, srtCol=45, + labRow=top[topgenes, 2]) + } else { + heatmap.2(topexp, scale="row", + main=paste("Heatmap:", unmake.names(contrastData[i])), + col="bluered", trace="none", density.info="none", + margin=c(8,6), lhei=c(2,10), dendrogram="column", + cexRow=0.7, cexCol=0.7, srtCol=45) + } + linkName <- paste0("Heatmap_", contrastData[i], ".pdf") + linkAddr <- paste0("heatmap_", contrastData[i], ".pdf") + linkData <- rbind(linkData, c(linkName, linkAddr)) + invisible(dev.off()) + # PNG of MD and Volcano png(mdvolOutPng[i], width=1200, height=600) par(mfrow=c(1, 2), mar=c(5,4,2,2)+0.1, oma=c(0,0,3,0)) limma::plotMD(fit, status=status[, i], coef=i, main="MD Plot", hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), xlab="Average Expression", ylab="logFC") - abline(h=0, col="grey", lty=2) # Volcano plots if (haveAnno) { # labels must be in second column currently limma::volcanoplot(fit, coef=i, main="Volcano Plot", - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes[, 2]) } else { limma::volcanoplot(fit, coef=i, main="Volcano Plot", - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes$GeneID) }
--- a/limma_voom.xml Tue May 08 18:12:40 2018 -0400 +++ b/limma_voom.xml Wed May 09 13:27:14 2018 -0400 @@ -1,4 +1,4 @@ -<tool id="limma_voom" name="limma" version="3.34.9.2"> +<tool id="limma_voom" name="limma" version="3.34.9.3"> <description> Perform differential expression with limma-voom or limma-trend </description> @@ -10,10 +10,11 @@ <requirement type="package" version="0.5.0">r-scales</requirement> <requirement type="package" version="0.2.15">r-rjson</requirement> <requirement type="package" version="1.20.0">r-getopt</requirement> + <requirement type="package" version="3.0.1">r-gplots</requirement> </requirements> <version_command><![CDATA[ -echo $(R --version | grep version | grep -v GNU)", limma version" $(R --vanilla --slave -e "library(limma); cat(sessionInfo()\$otherPkgs\$limma\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", edgeR version" $(R --vanilla --slave -e "library(edgeR); cat(sessionInfo()\$otherPkgs\$edgeR\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", statmod version" $(R --vanilla --slave -e "library(statmod); cat(sessionInfo()\$otherPkgs\$statmod\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", scales version" $(R --vanilla --slave -e "library(scales); cat(sessionInfo()\$otherPkgs\$scales\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", rjson version" $(R --vanilla --slave -e "library(rjson); cat(sessionInfo()\$otherPkgs\$rjson\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", getopt version" $(R --vanilla --slave -e "library(getopt); cat(sessionInfo()\$otherPkgs\$getopt\$Version)" 2> /dev/null | grep -v -i "WARNING: ") +echo $(R --version | grep version | grep -v GNU)", limma version" $(R --vanilla --slave -e "library(limma); cat(sessionInfo()\$otherPkgs\$limma\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", edgeR version" $(R --vanilla --slave -e "library(edgeR); cat(sessionInfo()\$otherPkgs\$edgeR\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", statmod version" $(R --vanilla --slave -e "library(statmod); cat(sessionInfo()\$otherPkgs\$statmod\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", scales version" $(R --vanilla --slave -e "library(scales); cat(sessionInfo()\$otherPkgs\$scales\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", rjson version" $(R --vanilla --slave -e "library(rjson); cat(sessionInfo()\$otherPkgs\$rjson\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", getopt version" $(R --vanilla --slave -e "library(getopt); cat(sessionInfo()\$otherPkgs\$getopt\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", gplots version" $(R --vanilla --slave -e "library(gplots); cat(sessionInfo()\$otherPkgs\$gplots\$Version)" 2> /dev/null | grep -v -i "WARNING: ") ]]></version_command> <command detect_errors="exit_code"><![CDATA[ @@ -82,7 +83,7 @@ -l '$adv.lfc' -p '$adv.pVal' -d '$adv.pAdjust' --v '$adv.volgenes' +-G '$adv.topgenes' #if $adv.treat: -T #end if @@ -277,9 +278,9 @@ <param name="treat" type="boolean" truevalue="1" falsevalue="0" checked="False" label="Test significance relative to a fold-change threshold (TREAT)" help="If you want to apply a cut-off on a fold change the TREAT function can be used, see Help section below. Default: No"/> - <param name="volgenes" type="integer" value="10" min="0" - label="Number of genes to highlight in Volcano plot" - help="The top DE genes will be highlighted in the Volcano plot for each contrast. Default: 10."/> + <param name="topgenes" type="integer" value="10" min="0" + label="Number of genes to highlight in Volcano plot and Heatmap" + help="The top DE genes will be highlighted in the Volcano plot for each contrast and output in a heatmap PDF. Default: 10."/> <param name="normalisationOption" type="select" label="Normalisation Method" help="Default: TMM"> <option value="TMM" selected="true">TMM</option> <option value="RLE">RLE</option> @@ -317,6 +318,7 @@ <param name="contrast" value="WT-Mut" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output_collection name="outTables" count="2"> <element name="limma-voom_Mut-WT" ftype="tabular" > <assert_contents> @@ -352,6 +354,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output_collection name="outTables" count="1"> <element name="limma-voom_Mut-WT" ftype="tabular" > <assert_contents> @@ -375,6 +378,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output name="outReport" > <assert_contents> <has_text text="RData" /> @@ -398,6 +402,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output_collection name="outTables" count="1" > <element name="limma-voom_Mut-WT" ftype="tabular" > <assert_contents> @@ -417,6 +422,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output_collection name="outTables" count="1"> <element name="limma-voom_Mut-WT" ftype="tabular" > <assert_contents> @@ -439,6 +445,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <output_collection name="outTables" count="2"> <element name="limma-voom_Mut-WT" ftype="tabular" > <assert_contents> @@ -491,6 +498,7 @@ <repeat name="rep_contrast"> <param name="contrast" value="WT-Mut" /> </repeat> + <param name="topgenes" value="6" /> <param name="normCounts" value="true" /> <output_collection name="outTables" count="3"> <element name="limma-voom_Mut-WT" ftype="tabular" > @@ -525,6 +533,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <param name="de_select" value="trend" /> <param name="rdaOption" value="true" /> <output name="outReport" > @@ -555,6 +564,7 @@ <param name="contrast" value="Mut-WT" /> </repeat> <param name="normalisationOption" value="TMM" /> + <param name="topgenes" value="6" /> <param name="de_select" value="trend" /> <param name="rdaOption" value="true" /> <output name="outReport" >
--- a/test-data/out_rscript.txt Tue May 08 18:12:40 2018 -0400 +++ b/test-data/out_rscript.txt Wed May 09 13:27:14 2018 -0400 @@ -22,13 +22,17 @@ # robOpt", "b", 0, "logical" -String specifying if robust options should be used # trend", "t", 1, "double" -Float for prior.count if limma-trend is used instead of voom # weightOpt", "w", 0, "logical" -String specifying if voomWithQualityWeights should be used -# volhiOpt", "v", 1, "integer" -Integer specifying number of genes to highlight in volcano plot +# topgenes", "G", 1, "integer" -Integer specifying no. of genes to highlight in volcano and heatmap # treatOpt", "T", 0, "logical" -String specifying if TREAT function shuld be used # # OUT: +# Density Plots (if filtering) +# Box Plots (if normalising) # MDS Plot # Voom/SA plot # MD Plot +# Volcano Plot +# Heatmap # Expression Table # HTML file linking to the ouputs # Optional: @@ -37,7 +41,7 @@ # # # Author: Shian Su - registertonysu@gmail.com - Jan 2014 -# Modified by: Maria Doyle - Jun 2017, Jan 2018 +# Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018 # Record starting time timeStart <- as.character(Sys.time()) @@ -50,6 +54,7 @@ library(limma, quietly=TRUE, warn.conflicts=FALSE) library(scales, quietly=TRUE, warn.conflicts=FALSE) library(getopt, quietly=TRUE, warn.conflicts=FALSE) +library(gplots, quietly=TRUE, warn.conflicts=FALSE) ################################################################################ ### Function Declaration @@ -169,7 +174,7 @@ "robOpt", "b", 0, "logical", "trend", "t", 1, "double", "weightOpt", "w", 0, "logical", - "volhiOpt", "v", 1, "integer", + "topgenes", "G", 1, "integer", "treatOpt", "T", 0, "logical"), byrow=TRUE, ncol=4) opt <- getopt(spec) @@ -330,11 +335,13 @@ mdOutPdf <- character() volOutPdf <- character() +heatOutPdf <- character() mdvolOutPng <- character() topOut <- character() for (i in 1:length(contrastData)) { mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf")) volOutPdf[i] <- makeOut(paste0("volplot_", contrastData[i], ".pdf")) + heatOutPdf[i] <- makeOut(paste0("heatmap_", contrastData[i], ".pdf")) mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", contrastData[i], ".png")) topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv")) } @@ -691,12 +698,12 @@ # labels must be in second column currently limma::volcanoplot(fit, coef=i, main=paste("Volcano Plot:", unmake.names(contrastData[i])), - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes[, 2]) } else { limma::volcanoplot(fit, coef=i, - main=paste("Volcano Plot:", unmake.names(contrastData[i])),, - highlight=opt$volhiOpt, + main=paste("Volcano Plot:", unmake.names(contrastData[i])), + highlight=opt$topgenes, names=fit$genes$GeneID) } linkName <- paste0("VolcanoPlot_", contrastData[i], ".pdf") @@ -704,24 +711,53 @@ linkData <- rbind(linkData, c(linkName, linkAddr)) invisible(dev.off()) + # Plot Heatmap + topgenes <- rownames(top[1:opt$topgenes, ]) + if (wantTrend) { + topexp <- plotData[topgenes, ] + } else { + topexp <- plotData$E[topgenes, ] + } + + pdf(heatOutPdf[i]) + par(cex.main=0.8) + if (haveAnno) { + # labels must be in second column currently + heatmap.2(topexp, scale="row", + main=paste("Heatmap:", unmake.names(contrastData[i])), + col="bluered", trace="none", density.info="none", + margin=c(8,6), lhei=c(2,10), dendrogram="column", + cexRow=0.7, cexCol=0.7, srtCol=45, + labRow=top[topgenes, 2]) + } else { + heatmap.2(topexp, scale="row", + main=paste("Heatmap:", unmake.names(contrastData[i])), + col="bluered", trace="none", density.info="none", + margin=c(8,6), lhei=c(2,10), dendrogram="column", + cexRow=0.7, cexCol=0.7, srtCol=45) + } + linkName <- paste0("Heatmap_", contrastData[i], ".pdf") + linkAddr <- paste0("heatmap_", contrastData[i], ".pdf") + linkData <- rbind(linkData, c(linkName, linkAddr)) + invisible(dev.off()) + # PNG of MD and Volcano png(mdvolOutPng[i], width=1200, height=600) par(mfrow=c(1, 2), mar=c(5,4,2,2)+0.1, oma=c(0,0,3,0)) limma::plotMD(fit, status=status[, i], coef=i, main="MD Plot", hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), xlab="Average Expression", ylab="logFC") - abline(h=0, col="grey", lty=2) # Volcano plots if (haveAnno) { # labels must be in second column currently limma::volcanoplot(fit, coef=i, main="Volcano Plot", - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes[, 2]) } else { limma::volcanoplot(fit, coef=i, main="Volcano Plot", - highlight=opt$volhiOpt, + highlight=opt$topgenes, names=fit$genes$GeneID) }