# HG changeset patch
# User iuc
# Date 1525886834 14400
# Node ID e6a4ff41af6b9ab029d7ea00dece56e04b26fd33
# Parent  39fa12a6d885d7be6788fe8a835a09db130d6d06
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit cf399638ebca4250bcc15f468238a9964de97b33

diff -r 39fa12a6d885 -r e6a4ff41af6b limma_voom.R
--- 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)
     }
 
diff -r 39fa12a6d885 -r e6a4ff41af6b limma_voom.xml
--- 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" >
diff -r 39fa12a6d885 -r e6a4ff41af6b test-data/out_rscript.txt
--- 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)
     }