Repository 'limma_voom'
hg clone https://toolshed.g2.bx.psu.edu/repos/iuc/limma_voom

Changeset 9:4255881bebb1 (2018-06-05)
Previous changeset 8:00a42b66e522 (2018-05-15) Next changeset 10:e26047c4562d (2018-06-08)
Commit message:
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 5f8052ab47b789da5ea47feb48e1b8e08b5fc1e0
modified:
limma_voom.R
limma_voom.xml
removed:
test-data/out_rscript.txt
b
diff -r 00a42b66e522 -r 4255881bebb1 limma_voom.R
--- a/limma_voom.R Tue May 15 02:36:36 2018 -0400
+++ b/limma_voom.R Tue Jun 05 04:28:56 2018 -0400
[
@@ -345,6 +345,7 @@
 stripOutPdf <- character()
 mdvolOutPng <- character()
 topOut <- character()
+glimmaOut <- character()
 for (i in 1:length(contrastData)) {
     con <- contrastData[i]
     con <- gsub("\\(|\\)", "", con)
@@ -354,6 +355,7 @@
     stripOutPdf[i] <- makeOut(paste0("stripcharts_", con, ".pdf"))
     mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", con, ".png"))
     topOut[i] <- makeOut(paste0(deMethod, "_", con, ".tsv"))
+    glimmaOut[i] <- makeOut(paste0("glimma_", con, "/MD-Plot.html"))
 }
 
 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv"))
@@ -806,6 +808,17 @@
     linkData <- rbind(linkData, c(linkName, linkAddr))
     invisible(dev.off())
 
+    # Generate Glimma interactive MD plot and table, requires annotation file (assumes gene names in 2nd column)
+    if (haveAnno) {
+        Glimma::glMDPlot(fit, coef=i, counts=y$counts, anno=y$genes, groups=factors[, 1],
+             status=status[, i], sample.cols=col.group,
+             main=paste("MD Plot:", unmake.names(con)), side.main=colnames(y$genes)[2],
+             folder=paste0("glimma_", unmake.names(con)), launch=FALSE)
+        linkName <- paste0("Glimma_MDPlot_", con, ".html")
+        linkAddr <- paste0("glimma_", con, "/MD-Plot.html")
+        linkData <- rbind(linkData, c(linkName, linkAddr))
+    }
+
     # Plot Volcano
     pdf(volOutPdf[i])
     if (haveAnno) {
@@ -1014,6 +1027,13 @@
     }
 }
 
+cata("<h4>Glimma Interactive Results:</h4>\n")
+    for (i in 1:nrow(linkData)) {
+        if (grepl("glimma", linkData$Link[i])) {
+            HtmlLink(linkData$Link[i], linkData$Label[i])
+        }
+    }
+
 cata("<p>Alt-click links to download file.</p>\n")
 cata("<p>Click floppy disc icon associated history item to download ")
 cata("all files.</p>\n")
b
diff -r 00a42b66e522 -r 4255881bebb1 limma_voom.xml
--- a/limma_voom.xml Tue May 15 02:36:36 2018 -0400
+++ b/limma_voom.xml Tue Jun 05 04:28:56 2018 -0400
[
@@ -1,4 +1,4 @@
-<tool id="limma_voom" name="limma" version="3.34.9.4">
+<tool id="limma_voom" name="limma" version="3.34.9.5">
     <description>
         Perform differential expression with limma-voom or limma-trend
     </description>
@@ -11,10 +11,11 @@
         <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>
+        <requirement type="package" version="1.6.0">bioconductor-glimma</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: ")", gplots version" $(R --vanilla --slave -e "library(gplots); cat(sessionInfo()\$otherPkgs\$gplots\$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: ")", Glimma version" $(R --vanilla --slave -e "library(Glimma); cat(sessionInfo()\$otherPkgs\$Glimma\$Version)" 2> /dev/null | grep -v -i "WARNING: ")
     ]]></version_command>
 
     <command detect_errors="exit_code"><![CDATA[
@@ -110,7 +111,10 @@
 mkdir ./output_dir
 
 &&
-cp '$outReport.files_path'/*.tsv output_dir/
+#if $anno.annoOpt=='yes':
+    cp -r ./glimma* '$outReport.files_path' &&
+#end if
+cp '$outReport.files_path'/*tsv output_dir/
 
 #if $out.rscript:
     && cp '$__tool_directory__/limma_voom.R' '$rscript'
@@ -194,7 +198,7 @@
         <!-- Gene Annotations -->
         <conditional name="anno">
             <param name="annoOpt" type="select" label="Use Gene Annotations?"
-                    help="If you provide an annotation file, annotations will be added to the table(s) of differential expression results to provide descriptions for each gene, and used to label the top genes in the Volcano plot. See Help section below.">
+                    help="If you provide an annotation file, annotations will be added to the table(s) of differential expression results to provide descriptions for each gene, and used to label the top genes in the Volcano plot. An interactive Glimma MD plot and table will also be generated. See Help section below.">
                 <option value="no">No</option>
                 <option value="yes">Yes</option>
             </param>
@@ -397,7 +401,11 @@
                     <has_text text="RData" />
                 </assert_contents>
             </output>
-            <output name="rscript" value="out_rscript.txt"/>
+            <output name="rscript">
+                <assert_contents>
+                    <has_text_matching expression="Task run time" />
+                </assert_contents>
+            </output>
         </test>
         <!-- Ensure secondary factors work -->
         <test>
@@ -601,7 +609,7 @@
 
 **What it does**
 
-Given a matrix of counts (e.g. from featureCounts) and optional information about the genes, performs differential expression (DE) using the limma_ Bioconductor package and produces plots and tables useful in DE analysis.
+Given a matrix of counts (e.g. from featureCounts) and optional information about the genes, performs differential expression (DE) using the limma_ Bioconductor package and produces plots and tables useful in DE analysis. If an annotation file is provided, interactive Glimma_ plots and a table of differentially expressed genes will also be generated. See an example workflow here_.
 
 In the `limma approach`_ to RNA-seq, read counts are converted to log2-counts-per-million (logCPM) and the mean-variance relationship is modelled either with precision weights or with an empirical Bayes prior trend. The precision weights approach is called “voom” and the prior trend approach is called “limma-trend”. For more information, see the Help section below.
 
@@ -647,7 +655,7 @@
 **Gene Annotations:**
 Optional input for gene annotations, this can contain more
 information about the genes than just an ID number. The annotations will
-be available in the differential expression results table and the optional normalised counts table. The file must contain a header row and have the gene IDs in the first column. The second column will be used to label the genes in the Volcano plot instead of the default Gene IDs. The number of rows should match that of the counts files, add NA for any gene IDs with no annotation. The Galaxy tool **annotateMyIDs** can be used to obtain annotations for human, mouse, fly and zebrafish.
+be available in the differential expression results table and the optional normalised counts table. They will also be used to generate interactive Glimma_ MD plots and table of differential expression, a link to the Glimma plots will be provided in the report. The input annotation file must contain a header row and have the gene IDs in the first column. The second column will be used to label the genes in the Volcano plot and interactive Glimma plots, additional columns will be available in the Glimma interactive table. The number of rows should match that of the counts files, add NA for any gene IDs with no annotation. The Galaxy tool **annotateMyIDs** can be used to obtain annotations for human, mouse, fly and zebrafish.
 
 Example:
 
@@ -783,6 +791,7 @@
 
     * a table of differentially expressed genes for each contrast of interest
     * a HTML report with plots and additional information
+    * an interactive Glimma MD plot and table (if annotation file provided)
 
 Optionally, under **Output Options** you can choose to output
 
@@ -844,6 +853,8 @@
 Please report problems or suggestions to: su.s@wehi.edu.au
 
 .. _limma: http://www.bioconductor.org/packages/release/bioc/html/limma.html
+.. _Glimma: https://bioconductor.org/packages/release/bioc/html/Glimma.html
+.. _here: https://f1000research.com/articles/5-1408/v2
 .. _limma approach: https://www.ncbi.nlm.nih.gov/pubmed/25605792
 .. _limma User's Guide: http://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf
 .. _edgeR: http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
b
diff -r 00a42b66e522 -r 4255881bebb1 test-data/out_rscript.txt
--- a/test-data/out_rscript.txt Tue May 15 02:36:36 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
[
b'@@ -1,1211 +0,0 @@\n-# This tool takes in a matrix of feature counts as well as gene annotations and\n-# outputs a table of top expressions as well as various plots for differential\n-# expression analysis\n-#\n-# ARGS: htmlPath", "R", 1, "character"      -Path to html file linking to other outputs\n-#       outPath", "o", 1, "character"       -Path to folder to write all output to\n-#       filesPath", "j", 2, "character"     -JSON list object if multiple files input\n-#       matrixPath", "m", 2, "character"    -Path to count matrix\n-#       factFile", "f", 2, "character"      -Path to factor information file\n-#       factInput", "i", 2, "character"     -String containing factors if manually input\n-#       annoPath", "a", 2, "character"      -Path to input containing gene annotations\n-#       contrastData", "C", 1, "character"  -String containing contrasts of interest\n-#       cpmReq", "c", 2, "double"           -Float specifying cpm requirement\n-#       cntReq", "z", 2, "integer"          -Integer specifying minimum total count requirement\n-#       sampleReq", "s", 2, "integer"       -Integer specifying cpm requirement\n-#       normCounts", "x", 0, "logical"      -String specifying if normalised counts should be output\n-#       rdaOpt", "r", 0, "logical"          -String specifying if RData should be output\n-#       lfcReq", "l", 1, "double"           -Float specifying the log-fold-change requirement\n-#       pValReq", "p", 1, "double"          -Float specifying the p-value requirement\n-#       pAdjOpt", "d", 1, "character"       -String specifying the p-value adjustment method\n-#       normOpt", "n", 1, "character"       -String specifying type of normalisation used\n-#       robOpt", "b", 0, "logical"          -String specifying if robust options should be used\n-#       trend", "t", 1, "double"            -Float for prior.count if limma-trend is used instead of voom\n-#       weightOpt", "w", 0, "logical"       -String specifying if voomWithQualityWeights should be used\n-#       topgenes", "G", 1, "integer"        -Integer specifying no. of genes to highlight in volcano and heatmap\n-#       treatOpt", "T", 0, "logical"        -String specifying if TREAT function should be used\n-#       plots, "P", 1, "character"          -String specifying additional plots to be created\n-#\n-# OUT:\n-#       Density Plots (if filtering)\n-#       Box Plots (if normalising)\n-#       MDS Plot\n-#       Voom/SA plot\n-#       MD Plot\n-#       Volcano Plot\n-#       Heatmap\n-#       Expression Table\n-#       HTML file linking to the ouputs\n-# Optional:\n-#       Normalised counts Table\n-#       RData file\n-#\n-#\n-# Author: Shian Su - registertonysu@gmail.com - Jan 2014\n-# Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018\n-\n-# Record starting time\n-timeStart <- as.character(Sys.time())\n-\n-# Load all required libraries\n-library(methods, quietly=TRUE, warn.conflicts=FALSE)\n-library(statmod, quietly=TRUE, warn.conflicts=FALSE)\n-library(splines, quietly=TRUE, warn.conflicts=FALSE)\n-library(edgeR, quietly=TRUE, warn.conflicts=FALSE)\n-library(limma, quietly=TRUE, warn.conflicts=FALSE)\n-library(scales, quietly=TRUE, warn.conflicts=FALSE)\n-library(getopt, quietly=TRUE, warn.conflicts=FALSE)\n-library(gplots, quietly=TRUE, warn.conflicts=FALSE)\n-\n-################################################################################\n-### Function Declaration\n-################################################################################\n-# Function to sanitise contrast equations so there are no whitespaces\n-# surrounding the arithmetic operators, leading or trailing whitespace\n-sanitiseEquation <- function(equation) {\n-    equation <- gsub(" *[+] *", "+", equation)\n-    equation <- gsub(" *[-] *", "-", equation)\n-    equation <- gsub(" *[/] *", "/", equation)\n-    equation <- gsub(" *[*] *", "*", equation)\n-    equation <- gsub("^\\\\s+|\\\\s+$", "", equation)\n-    return(equation)\n-}\n-\n-# Function to sanitise group information\n-sanitiseGroups <- function(string) {\n-    s'..b'("<a href=\\"",\n-                  "http://www.bioconductor.org/packages/release/bioc/",\n-                  "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf",\n-                  "\\">", "edgeR User\'s Guide", "</a>")\n-\n-cit[1] <- paste("Please cite the following paper for this tool:")\n-\n-cit[2] <- paste("Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME,",\n-                "Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? ",\n-                "Modelling sample and observational level variability improves power ",\n-                "in RNA-seq analyses. Nucleic Acids Research, 43(15), e97.")\n-\n-cit[3] <- paste("Please cite the paper below for the limma software itself.",\n-                "Please also try to cite the appropriate methodology articles",\n-                "that describe the statistical methods implemented in limma,",\n-                "depending on which limma functions you are using. The",\n-                "methodology articles are listed in Section 2.1 of the",\n-                link[1],\n-                "Cite no. 3 only if sample weights were used.")\n-cit[4] <- paste("Smyth GK (2005). Limma: linear models for microarray data.",\n-                "In: \'Bioinformatics and Computational Biology Solutions using",\n-                "R and Bioconductor\'. R. Gentleman, V. Carey, S. doit,.",\n-                "Irizarry, W. Huber (eds), Springer, New York, pages 397-420.")\n-cit[5] <- paste("Please cite the first paper for the software itself and the",\n-                "other papers for the various original statistical methods",\n-                "implemented in edgeR.  See Section 1.2 in the", link[2],\n-                "for more detail.")\n-cit[6] <- paste("Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a",\n-                "Bioconductor package for differential expression analysis",\n-                "of digital gene expression data. Bioinformatics 26, 139-140")\n-cit[7] <- paste("Robinson MD and Smyth GK (2007). Moderated statistical tests",\n-                "for assessing differences in tag abundance. Bioinformatics",\n-                "23, 2881-2887")\n-cit[8] <- paste("Robinson MD and Smyth GK (2008). Small-sample estimation of",\n-                "negative binomial dispersion, with applications to SAGE data.",\n-                "Biostatistics, 9, 321-332")\n-cit[9] <- paste("McCarthy DJ, Chen Y and Smyth GK (2012). Differential",\n-                "expression analysis of multifactor RNA-Seq experiments with",\n-                "respect to biological variation. Nucleic Acids Research 40,",\n-                "4288-4297")\n-cit[10] <- paste("Law CW, Chen Y, Shi W, and Smyth GK (2014). Voom:",\n-                "precision weights unlock linear model analysis tools for",\n-                "RNA-seq read counts. Genome Biology 15, R29.")\n-cit[11] <- paste("Ritchie ME, Diyagama D, Neilson J, van Laar R,",\n-                "Dobrovic A, Holloway A and Smyth GK (2006).",\n-                "Empirical array quality weights for microarray data.",\n-                "BMC Bioinformatics 7, Article 261.")\n-cata("<h3>Citations</h3>\\n")\n-cata(cit[1], "\\n")\n-cata("<br>\\n")\n-cata(cit[2], "\\n")\n-\n-cata("<h4>limma</h4>\\n")\n-cata(cit[3], "\\n")\n-cata("<ol>\\n")\n-ListItem(cit[4])\n-ListItem(cit[10])\n-ListItem(cit[11])\n-cata("</ol>\\n")\n-\n-cata("<h4>edgeR</h4>\\n")\n-cata(cit[5], "\\n")\n-cata("<ol>\\n")\n-ListItem(cit[6])\n-ListItem(cit[7])\n-ListItem(cit[8])\n-ListItem(cit[9])\n-cata("</ol>\\n")\n-\n-cata("<p>Please report problems or suggestions to: su.s@wehi.edu.au</p>\\n")\n-\n-for (i in 1:nrow(linkData)) {\n-    if (grepl("session_info", linkData$Link[i])) {\n-        HtmlLink(linkData$Link[i], linkData$Label[i])\n-    }\n-}\n-\n-cata("<table border=\\"0\\">\\n")\n-cata("<tr>\\n")\n-TableItem("Task started at:"); TableItem(timeStart)\n-cata("</tr>\\n")\n-cata("<tr>\\n")\n-TableItem("Task ended at:"); TableItem(timeEnd)\n-cata("</tr>\\n")\n-cata("<tr>\\n")\n-TableItem("Task run time:"); TableItem(timeTaken)\n-cata("<tr>\\n")\n-cata("</table>\\n")\n-\n-cata("</body>\\n")\n-cata("</html>")\n'