# HG changeset patch
# User iuc
# Date 1517420742 18000
# Node ID 38aab66ae5cbc72b9c0d9df5ab266a0028e50061
# Parent  a330ddf4386195c2093e20cff7c76e382aecc0bb
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 1640914b9812b0482a3cf684f05465f8d9cfdc65

diff -r a330ddf43861 -r 38aab66ae5cb limma_voom.R
--- a/limma_voom.R	Thu Sep 07 05:27:27 2017 -0400
+++ b/limma_voom.R	Wed Jan 31 12:45:42 2018 -0500
@@ -2,31 +2,40 @@
 # outputs a table of top expressions as well as various plots for differential
 # expression analysis
 #
-# ARGS: 1.countPath       -Path to RData input containing counts
-#       2.annoPath        -Path to input containing gene annotations
-#       3.htmlPath        -Path to html file linking to other outputs
-#       4.outPath         -Path to folder to write all output to
-#       5.rdaOpt          -String specifying if RData should be saved
-#       6.normOpt         -String specifying type of normalisation used
-#       7.weightOpt       -String specifying usage of weights
-#       8.contrastData    -String containing contrasts of interest
-#       9.cpmReq          -Float specifying cpm requirement
-#       10.sampleReq      -Integer specifying cpm requirement
-#       11.pAdjOpt        -String specifying the p-value adjustment method
-#       12.pValReq        -Float specifying the p-value requirement
-#       13.lfcReq         -Float specifying the log-fold-change requirement
-#       14.normCounts     -String specifying if normalised counts should be output
-#       15.factPath       -Path to factor information file
-#       16.factorData     -Strings containing factor names and values if manually input 
+# ARGS: htmlPath", "R", 1, "character"      -Path to html file linking to other outputs
+#       outPath", "o", 1, "character"       -Path to folder to write all output to
+#       filesPath", "j", 2, "character"     -JSON list object if multiple files input
+#       matrixPath", "m", 2, "character"    -Path to count matrix
+#       factFile", "f", 2, "character"      -Path to factor information file
+#       factInput", "i", 2, "character"     -String containing factors if manually input
+#       annoPath", "a", 2, "character"      -Path to input containing gene annotations
+#       contrastData", "C", 1, "character"  -String containing contrasts of interest
+#       cpmReq", "c", 2, "double"           -Float specifying cpm requirement
+#       cntReq", "z", 2, "integer"          -Integer specifying minimum total count requirement
+#       sampleReq", "s", 2, "integer"       -Integer specifying cpm requirement
+#       normCounts", "x", 0, "logical"      -String specifying if normalised counts should be output
+#       rdaOpt", "r", 0, "logical"          -String specifying if RData should be output
+#       lfcReq", "l", 1, "double"           -Float specifying the log-fold-change requirement
+#       pValReq", "p", 1, "double"          -Float specifying the p-value requirement
+#       pAdjOpt", "d", 1, "character"       -String specifying the p-value adjustment method
+#       normOpt", "n", 1, "character"       -String specifying type of normalisation used
+#       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
 #
-# OUT:  Voom Plot
-#       BCV Plot
-#       MA Plot
+# OUT:
+#       MDS Plot
+#       Voom/SA plot
+#       MD Plot
 #       Expression Table
 #       HTML file linking to the ouputs
+# Optional:
+#       Normalised counts Table
+#       RData file
+#
 #
 # Author: Shian Su - registertonysu@gmail.com - Jan 2014
-# Modified by: Maria Doyle - Jun 2017
+# Modified by: Maria Doyle - Jun 2017, Jan 2018
 
 # Record starting time
 timeStart <- as.character(Sys.time())
@@ -38,9 +47,10 @@
 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
 library(limma, quietly=TRUE, warn.conflicts=FALSE)
 library(scales, quietly=TRUE, warn.conflicts=FALSE)
+library(getopt, quietly=TRUE, warn.conflicts=FALSE)
 
 if (packageVersion("limma") < "3.20.1") {
-  stop("Please update 'limma' to version >= 3.20.1 to run this tool")
+    stop("Please update 'limma' to version >= 3.20.1 to run this tool")
 }
 
 ################################################################################
@@ -49,185 +59,276 @@
 # Function to sanitise contrast equations so there are no whitespaces
 # surrounding the arithmetic operators, leading or trailing whitespace
 sanitiseEquation <- function(equation) {
-  equation <- gsub(" *[+] *", "+", equation)
-  equation <- gsub(" *[-] *", "-", equation)
-  equation <- gsub(" *[/] *", "/", equation)
-  equation <- gsub(" *[*] *", "*", equation)
-  equation <- gsub("^\\s+|\\s+$", "", equation)
-  return(equation)
+    equation <- gsub(" *[+] *", "+", equation)
+    equation <- gsub(" *[-] *", "-", equation)
+    equation <- gsub(" *[/] *", "/", equation)
+    equation <- gsub(" *[*] *", "*", equation)
+    equation <- gsub("^\\s+|\\s+$", "", equation)
+    return(equation)
 }
 
 # Function to sanitise group information
 sanitiseGroups <- function(string) {
-  string <- gsub(" *[,] *", ",", string)
-  string <- gsub("^\\s+|\\s+$", "", string)
-  return(string)
+    string <- gsub(" *[,] *", ",", string)
+    string <- gsub("^\\s+|\\s+$", "", string)
+    return(string)
 }
 
 # Function to change periods to whitespace in a string
 unmake.names <- function(string) {
-  string <- gsub(".", " ", string, fixed=TRUE)
-  return(string)
+    string <- gsub(".", " ", string, fixed=TRUE)
+    return(string)
 }
 
 # Generate output folder and paths
 makeOut <- function(filename) {
-  return(paste0(outPath, "/", filename))
+    return(paste0(opt$outPath, "/", filename))
 }
 
 # Generating design information
 pasteListName <- function(string) {
-  return(paste0("factors$", string))
+    return(paste0("factors$", string))
 }
 
 # Create cata function: default path set, default seperator empty and appending
 # true by default (Ripped straight from the cat function with altered argument
 # defaults)
-cata <- function(..., file = htmlPath, sep = "", fill = FALSE, labels = NULL, 
-                 append = TRUE) {
-  if (is.character(file)) 
-    if (file == "") 
-      file <- stdout()
-  else if (substring(file, 1L, 1L) == "|") {
-    file <- pipe(substring(file, 2L), "w")
-    on.exit(close(file))
-  }
-  else {
-    file <- file(file, ifelse(append, "a", "w"))
-    on.exit(close(file))
-  }
-  .Internal(cat(list(...), file, sep, fill, labels, append))
+cata <- function(..., file = opt$htmlPath, sep = "", fill = FALSE, labels = NULL,
+               append = TRUE) {
+    if (is.character(file))
+        if (file == "")
+            file <- stdout()
+        else if (substring(file, 1L, 1L) == "|") {
+            file <- pipe(substring(file, 2L), "w")
+            on.exit(close(file))
+        }
+        else {
+        file <- file(file, ifelse(append, "a", "w"))
+      on.exit(close(file))
+        }
+    .Internal(cat(list(...), file, sep, fill, labels, append))
 }
 
 # Function to write code for html head and title
 HtmlHead <- function(title) {
-  cata("<head>\n")
-  cata("<title>", title, "</title>\n")
-  cata("</head>\n")
+    cata("<head>\n")
+    cata("<title>", title, "</title>\n")
+    cata("</head>\n")
 }
 
 # Function to write code for html links
 HtmlLink <- function(address, label=address) {
-  cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
+    cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
 }
 
 # Function to write code for html images
 HtmlImage <- function(source, label=source, height=600, width=600) {
-  cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
-  cata("\" width=\"", width, "\"/>\n")
+    cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
+    cata("\" width=\"", width, "\"/>\n")
 }
 
 # Function to write code for html list items
 ListItem <- function(...) {
-  cata("<li>", ..., "</li>\n")
+    cata("<li>", ..., "</li>\n")
 }
 
 TableItem <- function(...) {
-  cata("<td>", ..., "</td>\n")
+    cata("<td>", ..., "</td>\n")
 }
 
 TableHeadItem <- function(...) {
-  cata("<th>", ..., "</th>\n")
+    cata("<th>", ..., "</th>\n")
 }
 
 ################################################################################
 ### Input Processing
 ################################################################################
 
-# Collects arguments from command line
-argv <- commandArgs(TRUE)
+# Collect arguments from command line
+args <- commandArgs(trailingOnly=TRUE)
 
-# Grab arguments
-countPath <- as.character(argv[1])
-annoPath <- as.character(argv[2])
-htmlPath <- as.character(argv[3])
-outPath <- as.character(argv[4])
-rdaOpt <- as.character(argv[5])
-normOpt <- as.character(argv[6])
-weightOpt <- as.character(argv[7])
-contrastData <- as.character(argv[8])
-cpmReq <- as.numeric(argv[9])
-sampleReq <- as.numeric(argv[10])
-pAdjOpt <- as.character(argv[11])
-pValReq <- as.numeric(argv[12])
-lfcReq <- as.numeric(argv[13])
-normCounts <- as.character(argv[14])
-factPath <- as.character(argv[15])
-# Process factors
-if (as.character(argv[16])=="None") {
-    factorData <- read.table(factPath, header=TRUE, sep="\t")
-    factors <- factorData[,-1, drop=FALSE]
-}  else { 
-    factorData <- list()
-    for (i in 16:length(argv)) {
-        newFact <- unlist(strsplit(as.character(argv[i]), split="::"))
-        factorData <- rbind(factorData, newFact)
-    } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor.
+# Get options, using the spec as defined by the enclosed list.
+# Read the options from the default: commandArgs(TRUE).
+spec <- matrix(c(
+    "htmlPath", "R", 1, "character",
+    "outPath", "o", 1, "character",
+    "filesPath", "j", 2, "character",
+    "matrixPath", "m", 2, "character",
+    "factFile", "f", 2, "character",
+    "factInput", "i", 2, "character",
+    "annoPath", "a", 2, "character",
+    "contrastData", "C", 1, "character",
+    "cpmReq", "c", 1, "double",
+    "totReq", "y", 0, "logical",
+    "cntReq", "z", 1, "integer",
+    "sampleReq", "s", 1, "integer",
+    "normCounts", "x", 0, "logical",
+    "rdaOpt", "r", 0, "logical",
+    "lfcReq", "l", 1, "double",
+    "pValReq", "p", 1, "double",
+    "pAdjOpt", "d", 1, "character",
+    "normOpt", "n", 1, "character",
+    "robOpt", "b", 0, "logical",
+    "trend", "t", 1, "double",
+    "weightOpt", "w", 0, "logical"),
+    byrow=TRUE, ncol=4)
+opt <- getopt(spec)
 
-    # Set the row names to be the name of the factor and delete first row
-    row.names(factorData) <- factorData[, 1]
-    factorData <- factorData[, -1]
-    factorData <- sapply(factorData, sanitiseGroups)
-    factorData <- sapply(factorData, strsplit, split=",")
-    factorData <- sapply(factorData, make.names)
-    # Transform factor data into data frame of R factor objects
-    factors <- data.frame(factorData)
 
+if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {
+    cat("A counts matrix (or a set of counts files) is required.\n")
+    q(status=1)
+}
+
+if (is.null(opt$cpmReq)) {
+    filtCPM <- FALSE
+} else {
+    filtCPM <- TRUE
 }
 
-# Process other arguments
-if (weightOpt=="yes") {
-  wantWeight <- TRUE
+if (is.null(opt$cntReq) || is.null(opt$sampleReq)) {
+    filtSmpCount <- FALSE
+} else {
+    filtSmpCount <- TRUE
+}
+
+if (is.null(opt$totReq)) {
+    filtTotCount <- FALSE
 } else {
-  wantWeight <- FALSE
+    filtTotCount <- TRUE
+}
+
+if (is.null(opt$rdaOpt)) {
+    wantRda <- FALSE
+} else {
+    wantRda <- TRUE
+}
+
+if (is.null(opt$annoPath)) {
+    haveAnno <- FALSE
+} else {
+    haveAnno <- TRUE
 }
 
-if (rdaOpt=="yes") {
-  wantRda <- TRUE
+if (is.null(opt$normCounts)) {
+    wantNorm <- FALSE
 } else {
-  wantRda <- FALSE
+    wantNorm <- TRUE
+}
+
+if (is.null(opt$robOpt)) {
+    wantRobust <- FALSE
+} else {
+    wantRobust <- TRUE
 }
 
-if (annoPath=="None") {
-  haveAnno <- FALSE
+if (is.null(opt$weightOpt)) {
+    wantWeight <- FALSE
 } else {
-  haveAnno <- TRUE
+    wantWeight <- TRUE
 }
 
-if (normCounts=="yes") {
-  wantNorm <- TRUE
+if (is.null(opt$trend)) {
+    wantTrend <- FALSE
+    deMethod <- "limma-voom"
 } else {
-  wantNorm <- FALSE
+    wantTrend <- TRUE
+    deMethod <- "limma-trend"
+    priorCount <- opt$trend
 }
 
 
+if (!is.null(opt$filesPath)) {
+    # Process the separate count files (adapted from DESeq2 wrapper)
+    library("rjson")
+    parser <- newJSONParser()
+    parser$addData(opt$filesPath)
+    factorList <- parser$getObject()
+    factors <- sapply(factorList, function(x) x[[1]])
+    filenamesIn <- unname(unlist(factorList[[1]][[2]]))
+    sampleTable <- data.frame(sample=basename(filenamesIn),
+                            filename=filenamesIn,
+                            row.names=filenamesIn,
+                            stringsAsFactors=FALSE)
+    for (factor in factorList) {
+        factorName <- factor[[1]]
+        sampleTable[[factorName]] <- character(nrow(sampleTable))
+        lvls <- sapply(factor[[2]], function(x) names(x))
+        for (i in seq_along(factor[[2]])) {
+            files <- factor[[2]][[i]][[1]]
+            sampleTable[files,factorName] <- lvls[i]
+        }
+        sampleTable[[factorName]] <- factor(sampleTable[[factorName]], levels=lvls)
+    }
+    rownames(sampleTable) <- sampleTable$sample
+    rem <- c("sample","filename")
+    factors <- sampleTable[, !(names(sampleTable) %in% rem), drop=FALSE]
+
+    #read in count files and create single table
+    countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)})
+    counts <- do.call("cbind", countfiles)
+
+} else {
+    # Process the single count matrix
+    counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
+    row.names(counts) <- counts[, 1]
+    counts <- counts[ , -1]
+    countsRows <- nrow(counts)
+
+    # Process factors
+    if (is.null(opt$factInput)) {
+            factorData <- read.table(opt$factFile, header=TRUE, sep="\t")
+            factors <- factorData[, -1, drop=FALSE]
+    }  else {
+            factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE))
+            factorData <- list()
+            for (fact in factors) {
+                newFact <- unlist(strsplit(fact, split="::"))
+                factorData <- rbind(factorData, newFact)
+            } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor.
+
+            # Set the row names to be the name of the factor and delete first row
+            row.names(factorData) <- factorData[, 1]
+            factorData <- factorData[, -1]
+            factorData <- sapply(factorData, sanitiseGroups)
+            factorData <- sapply(factorData, strsplit, split=",")
+            factorData <- sapply(factorData, make.names)
+            # Transform factor data into data frame of R factor objects
+            factors <- data.frame(factorData)
+    }
+}
+
+ # if annotation file provided
+if (haveAnno) {
+    geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
+}
+
 #Create output directory
-dir.create(outPath, showWarnings=FALSE)
+dir.create(opt$outPath, showWarnings=FALSE)
 
 # Split up contrasts seperated by comma into a vector then sanitise
-contrastData <- unlist(strsplit(contrastData, split=","))
+contrastData <- unlist(strsplit(opt$contrastData, split=","))
 contrastData <- sanitiseEquation(contrastData)
 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE)
 
-bcvOutPdf <- makeOut("bcvplot.pdf")
-bcvOutPng <- makeOut("bcvplot.png")
-mdsOutPdf <- makeOut("mdsplot.pdf")
-mdsOutPng <- makeOut("mdsplot.png")
-voomOutPdf <- makeOut("voomplot.pdf")
-voomOutPng <- makeOut("voomplot.png") 
+
+mdsOutPdf <- makeOut("mdsplot_nonorm.pdf")
+mdsOutPng <- makeOut("mdsplot_nonorm.png")
+nmdsOutPdf <- makeOut("mdsplot.pdf")
+nmdsOutPng <- makeOut("mdsplot.png")
 maOutPdf <- character()   # Initialise character vector
 maOutPng <- character()
 topOut <- character()
 for (i in 1:length(contrastData)) {
-  maOutPdf[i] <- makeOut(paste0("maplot_", contrastData[i], ".pdf"))
-  maOutPng[i] <- makeOut(paste0("maplot_", contrastData[i], ".png"))
-  topOut[i] <- makeOut(paste0("limma-voom_", contrastData[i], ".tsv"))
-}                         # Save output paths for each contrast as vectors
-normOut <- makeOut("limma-voom_normcounts.tsv")
-rdaOut <- makeOut("RData.rda")
+    maOutPdf[i] <- makeOut(paste0("maplot_", contrastData[i], ".pdf"))
+    maOutPng[i] <- makeOut(paste0("maplot_", contrastData[i], ".png"))
+    topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv"))
+}
+normOut <- makeOut(paste0(deMethod, "_normcounts.tsv"))
+rdaOut <- makeOut(paste0(deMethod, "_analysis.RData"))
 sessionOut <- makeOut("session_info.txt")
 
-# Initialise data for html links and images, data frame with columns Label and 
+# Initialise data for html links and images, data frame with columns Label and
 # Link
 linkData <- data.frame(Label=character(), Link=character(),
                        stringsAsFactors=FALSE)
@@ -238,21 +339,13 @@
 upCount <- numeric()
 downCount <- numeric()
 flatCount <- numeric()
-                        
-# Read in counts and geneanno data
-counts <- read.table(countPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
-row.names(counts) <- counts[, 1]
-counts <- counts[ , -1]
-countsRows <- nrow(counts)
-if (haveAnno) {
-  geneanno <- read.table(annoPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
-}
 
 ################################################################################
 ### Data Processing
 ################################################################################
 
 # Extract counts and annotation data
+print("Extracting counts")
 data <- list()
 data$counts <- counts
 if (haveAnno) {
@@ -261,12 +354,24 @@
   data$genes <- data.frame(GeneID=row.names(counts))
 }
 
-# Filter out genes that do not have a required cpm in a required number of
-# samples
+# If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of
+# samples. Default is no filtering
 preFilterCount <- nrow(data$counts)
-sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
-data$counts <- data$counts[sel, ]
-data$genes <- data$genes[sel, ,drop = FALSE]
+
+if (filtCPM || filtSmpCount || filtTotCount) {
+
+    if (filtTotCount) {
+        keep <- rowSums(data$counts) >= opt$cntReq
+    } else if (filtSmpCount) {
+        keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
+    } else if (filtCPM) {
+        keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq
+    }
+
+    data$counts <- data$counts[keep, ]
+    data$genes <- data$genes[keep, , drop=FALSE]
+}
+
 postFilterCount <- nrow(data$counts)
 filteredCount <- preFilterCount-postFilterCount
 
@@ -276,173 +381,231 @@
 
 
 # Generating the DGEList object "data"
+print("Generating DGEList object")
 data$samples <- sampleanno
 data$samples$lib.size <- colSums(data$counts)
 data$samples$norm.factors <- 1
 row.names(data$samples) <- colnames(data$counts)
 data <- new("DGEList", data)
 
+print("Generating Design")
+# Name rows of factors according to their sample
+row.names(factors) <- names(data$counts)
 factorList <- sapply(names(factors), pasteListName)
-
-formula <- "~0" 
+formula <- "~0"
 for (i in 1:length(factorList)) {
     formula <- paste(formula,factorList[i], sep="+")
 }
-
 formula <- formula(formula)
 design <- model.matrix(formula)
-
 for (i in 1:length(factorList)) {
     colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
 }
 
-# Calculating normalising factor, estimating dispersion
-data <- calcNormFactors(data, method=normOpt)
-#data <- estimateDisp(data, design=design, robust=TRUE)
+# Calculating normalising factors
+print("Calculating Normalisation Factors")
+data <- calcNormFactors(data, method=opt$normOpt)
 
 # Generate contrasts information
+print("Generating Contrasts")
 contrasts <- makeContrasts(contrasts=contrastData, levels=design)
 
-# Name rows of factors according to their sample
-row.names(factors) <- names(data$counts)
-
 ################################################################################
 ### Data Output
 ################################################################################
-
-# BCV Plot
-#png(bcvOutPng, width=600, height=600)
-#plotBCV(data, main="BCV Plot")
-#imageData[1, ] <- c("BCV Plot", "bcvplot.png")
-#invisible(dev.off())
-
-#pdf(bcvOutPdf)
-#plotBCV(data, main="BCV Plot")
-#invisible(dev.off())
-
-if (wantWeight) {
-  # Creating voom data object and plot
-  png(voomOutPng, width=1000, height=600)
-  vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
-  imageData[1, ] <- c("Voom Plot", "voomplot.png")
-  invisible(dev.off())
-  
-  pdf(voomOutPdf, width=14)
-  vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
-  linkData[1, ] <- c("Voom Plot (.pdf)", "voomplot.pdf")
-  invisible(dev.off())
-  
-  # Generating fit data and top table with weights
-  wts <- vData$weights
-  voomFit <- lmFit(vData, design, weights=wts)
-  
-} else {
-  # Creating voom data object and plot
-  png(voomOutPng, width=600, height=600)
-  vData <- voom(data, design=design, plot=TRUE)
-  imageData[1, ] <- c("Voom Plot", "voomplot.png")
-  invisible(dev.off())
-  
-  pdf(voomOutPdf)
-  vData <- voom(data, design=design, plot=TRUE)
-  linkData[1, ] <- c("Voom Plot (.pdf)", "voomplot.pdf")
-  invisible(dev.off())
-  
-  # Generate voom fit
-  voomFit <- lmFit(vData, design)
-  
-}
-
- # Save normalised counts (log2cpm)
-if (wantNorm) {   
-    norm_counts <- data.frame(vData$genes, vData$E)
-    write.table (norm_counts, file=normOut, row.names=FALSE, sep="\t")
-    linkData <- rbind(linkData, c("limma-voom_normcounts.tsv", "limma-voom_normcounts.tsv"))
-}
-
-# Fit linear model and estimate dispersion with eBayes
-voomFit <- contrasts.fit(voomFit, contrasts)
-voomFit <- eBayes(voomFit)
-
 # Plot MDS
+print("Generating MDS plot")
 labels <- names(counts)
 png(mdsOutPng, width=600, height=600)
 # Currently only using a single factor
-plotMDS(vData, labels=labels, col=as.numeric(factors[, 1]), cex=0.8)
-imgName <- "Voom Plot"
+plotMDS(data, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main="MDS Plot (unnormalised)")
+imageData[1, ] <- c("MDS Plot (unnormalised)", "mdsplot_nonorm.png")
+invisible(dev.off())
+
+pdf(mdsOutPdf)
+plotMDS(data, labels=labels, cex=0.5)
+linkData[1, ] <- c("MDS Plot (unnormalised).pdf", "mdsplot_nonorm.pdf")
+invisible(dev.off())
+
+if (wantTrend) {
+    # limma-trend approach
+    logCPM <- cpm(data, log=TRUE, prior.count=opt$trend)
+    fit <- lmFit(logCPM, design)
+    fit <- contrasts.fit(fit, contrasts)
+    if (wantRobust) {
+        fit <- eBayes(fit, trend=TRUE, robust=TRUE)
+    } else {
+        fit <- eBayes(fit, trend=TRUE, robust=FALSE)
+    }
+    # plot fit with plotSA
+    saOutPng <- makeOut("saplot.png")
+    saOutPdf <- makeOut("saplot.pdf")
+
+    png(saOutPng, width=600, height=600)
+    plotSA(fit, main="SA Plot")
+    imgName <- "SA Plot.png"
+    imgAddr <- "saplot.png"
+    imageData <- rbind(imageData, c(imgName, imgAddr))
+    invisible(dev.off())
+
+    pdf(saOutPdf, width=14)
+    plotSA(fit, main="SA Plot")
+    linkName <- paste0("SA Plot.pdf")
+    linkAddr <- paste0("saplot.pdf")
+    linkData <- rbind(linkData, c(linkName, linkAddr))
+    invisible(dev.off())
+
+    plotData <- logCPM
+
+    # Save normalised counts (log2cpm)
+    if (wantNorm) {
+        write.table(logCPM, file=normOut, row.names=TRUE, sep="\t")
+        linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv"))))
+    }
+} else {
+    # limma-voom approach
+    voomOutPdf <- makeOut("voomplot.pdf")
+    voomOutPng <- makeOut("voomplot.png")
+
+    if (wantWeight) {
+        # Creating voom data object and plot
+        png(voomOutPng, width=1000, height=600)
+        vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
+        imgName <- "Voom Plot.png"
+        imgAddr <- "voomplot.png"
+        imageData <- rbind(imageData, c(imgName, imgAddr))
+        invisible(dev.off())
+
+        pdf(voomOutPdf, width=14)
+        vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
+        linkName <- paste0("Voom Plot.pdf")
+        linkAddr <- paste0("voomplot.pdf")
+        linkData <- rbind(linkData, c(linkName, linkAddr))
+        invisible(dev.off())
+
+        # Generating fit data and top table with weights
+        wts <- vData$weights
+        voomFit <- lmFit(vData, design, weights=wts)
+
+    } else {
+        # Creating voom data object and plot
+        png(voomOutPng, width=600, height=600)
+        vData <- voom(data, design=design, plot=TRUE)
+        imgName <- "Voom Plot"
+        imgAddr <- "voomplot.png"
+        imageData <- rbind(imageData, c(imgName, imgAddr))
+        invisible(dev.off())
+
+        pdf(voomOutPdf)
+        vData <- voom(data, design=design, plot=TRUE)
+        linkName <- paste0("Voom Plot.pdf")
+        linkAddr <- paste0("voomplot.pdf")
+        linkData <- rbind(linkData, c(linkName, linkAddr))
+        invisible(dev.off())
+
+        # Generate voom fit
+        voomFit <- lmFit(vData, design)
+    }
+
+     # Save normalised counts (log2cpm)
+    if (wantNorm) {
+        norm_counts <- data.frame(vData$genes, vData$E)
+        write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t")
+        linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv"))))
+    }
+
+    # Fit linear model and estimate dispersion with eBayes
+    voomFit <- contrasts.fit(voomFit, contrasts)
+    if (wantRobust) {
+        fit <- eBayes(voomFit, robust=TRUE)
+    } else {
+        fit <- eBayes(voomFit, robust=FALSE)
+    }
+    plotData <- vData
+}
+
+print("Generating normalised MDS plot")
+png(nmdsOutPng, width=600, height=600)
+# Currently only using a single factor
+plotMDS(plotData, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main="MDS Plot (normalised)")
+imgName <- "MDS Plot (normalised)"
 imgAddr <- "mdsplot.png"
 imageData <- rbind(imageData, c(imgName, imgAddr))
 invisible(dev.off())
 
-pdf(mdsOutPdf)
-plotMDS(vData, labels=labels, cex=0.5)
-linkName <- paste0("MDS Plot (.pdf)")
+pdf(nmdsOutPdf)
+plotMDS(plotData, labels=labels, cex=0.5)
+linkName <- paste0("MDS Plot (normalised).pdf")
 linkAddr <- paste0("mdsplot.pdf")
 linkData <- rbind(linkData, c(linkName, linkAddr))
 invisible(dev.off())
 
 
+print("Generating DE results")
+status = decideTests(fit, adjust.method=opt$pAdjOpt, p.value=opt$pValReq,
+                       lfc=opt$lfcReq)
+sumStatus <- summary(status)
+
 for (i in 1:length(contrastData)) {
+    # Collect counts for differential expression
+    upCount[i] <- sumStatus["Up", i]
+    downCount[i] <- sumStatus["Down", i]
+    flatCount[i] <- sumStatus["NotSig", i]
+
+    # Write top expressions table
+    top <- topTable(fit, coef=i, number=Inf, sort.by="P")
+    if (wantTrend) {
+        write.table(top, file=topOut[i], row.names=TRUE, sep="\t")
+    } else {
+        write.table(top, file=topOut[i], row.names=FALSE, sep="\t")
+    }
+
+    linkName <- paste0(deMethod, "_", contrastData[i], ".tsv")
+    linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv")
+    linkData <- rbind(linkData, c(linkName, linkAddr))
 
-  status = decideTests(voomFit[, i], adjust.method=pAdjOpt, p.value=pValReq,
-                       lfc=lfcReq)
-                       
-  sumStatus <- summary(status)
-  
-  # Collect counts for differential expression
-  upCount[i] <- sumStatus["1",]
-  downCount[i] <- sumStatus["-1",]
-  flatCount[i] <- sumStatus["0",]
-                       
-  # Write top expressions table
-  top <- topTable(voomFit, coef=i, number=Inf, sort.by="P")
-  write.table(top, file=topOut[i], row.names=FALSE, sep="\t")
-  
-  linkName <- paste0("limma-voom_", contrastData[i], ".tsv")
-  linkAddr <- paste0("limma-voom_", contrastData[i], ".tsv")
-  linkData <- rbind(linkData, c(linkName, linkAddr))
-  
-  # Plot MA (log ratios vs mean average) using limma package on weighted 
-  pdf(maOutPdf[i])
-  limma::plotMA(voomFit, status=status, coef=i,
-                main=paste("MA Plot:", unmake.names(contrastData[i])), 
-                col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
-                xlab="Average Expression", ylab="logFC")
-  
-  abline(h=0, col="grey", lty=2)
-  
-  linkName <- paste0("MA Plot_", contrastData[i], " (.pdf)")
-  linkAddr <- paste0("maplot_", contrastData[i], ".pdf")
-  linkData <- rbind(linkData, c(linkName, linkAddr))
-  invisible(dev.off())
-  
-  png(maOutPng[i], height=600, width=600)
-  limma::plotMA(voomFit, status=status, coef=i,
-                main=paste("MA Plot:", unmake.names(contrastData[i])), 
-                col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
-                xlab="Average Expression", ylab="logFC")
-  
-  abline(h=0, col="grey", lty=2)
-  
-  imgName <- paste0("MA Plot_", contrastData[i])
-  imgAddr <- paste0("maplot_", contrastData[i], ".png")
-  imageData <- rbind(imageData, c(imgName, imgAddr))
-  invisible(dev.off())
+    # Plot MA (log ratios vs mean average) using limma package on weighted
+    pdf(maOutPdf[i])
+    limma::plotMD(fit, status=status, coef=i,
+                  main=paste("MA Plot:", unmake.names(contrastData[i])),
+                  col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
+                  xlab="Average Expression", ylab="logFC")
+
+    abline(h=0, col="grey", lty=2)
+
+    linkName <- paste0("MA Plot_", contrastData[i], " (.pdf)")
+    linkAddr <- paste0("maplot_", contrastData[i], ".pdf")
+    linkData <- rbind(linkData, c(linkName, linkAddr))
+    invisible(dev.off())
+
+    png(maOutPng[i], height=600, width=600)
+    limma::plotMD(fit, status=status, coef=i,
+                  main=paste("MA Plot:", unmake.names(contrastData[i])),
+                  col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
+                  xlab="Average Expression", ylab="logFC")
+
+    abline(h=0, col="grey", lty=2)
+
+    imgName <- paste0("MA Plot_", contrastData[i])
+    imgAddr <- paste0("maplot_", contrastData[i], ".png")
+    imageData <- rbind(imageData, c(imgName, imgAddr))
+    invisible(dev.off())
 }
 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
 row.names(sigDiff) <- contrastData
 
 # Save relevant items as rda object
 if (wantRda) {
-  if (wantWeight) {
-    save(data, status, vData, labels, factors, wts, voomFit, top, contrasts, 
-         design,
-         file=rdaOut, ascii=TRUE)
-  } else {
-    save(data, status, vData, labels, factors, voomFit, top, contrasts, design,
-         file=rdaOut, ascii=TRUE)
-  }
-  linkData <- rbind(linkData, c("RData (.rda)", "RData.rda"))
+    print("Saving RData")
+    if (wantWeight) {
+      save(data, status, plotData, labels, factors, wts, fit, top, contrasts,
+           design,
+           file=rdaOut, ascii=TRUE)
+    } else {
+      save(data, status, plotData, labels, factors, fit, top, contrasts, design,
+           file=rdaOut, ascii=TRUE)
+    }
+    linkData <- rbind(linkData, c((paste0(deMethod, "_analysis.RData")), (paste0(deMethod, "_analysis.RData"))))
 }
 
 # Record session info
@@ -458,21 +621,21 @@
 ################################################################################
 
 # Clear file
-cat("", file=htmlPath)
+cat("", file=opt$htmlPath)
 
 cata("<html>\n")
 
 cata("<body>\n")
-cata("<h3>Limma-voom Analysis Output:</h3>\n")
+cata("<h3>Limma Analysis Output:</h3>\n")
 cata("PDF copies of JPEGS available in 'Plots' section.<br />\n")
 if (wantWeight) {
-  HtmlImage(imageData$Link[1], imageData$Label[1], width=1000)
+    HtmlImage(imageData$Link[1], imageData$Label[1], width=1000)
 } else {
-  HtmlImage(imageData$Link[1], imageData$Label[1])
+    HtmlImage(imageData$Link[1], imageData$Label[1])
 }
 
 for (i in 2:nrow(imageData)) {
-  HtmlImage(imageData$Link[i], imageData$Label[i])
+    HtmlImage(imageData$Link[i], imageData$Label[i])
 }
 
 cata("<h4>Differential Expression Counts:</h4>\n")
@@ -481,40 +644,40 @@
 cata("<tr>\n")
 TableItem()
 for (i in colnames(sigDiff)) {
-  TableHeadItem(i)
+    TableHeadItem(i)
 }
 cata("</tr>\n")
 for (i in 1:nrow(sigDiff)) {
-  cata("<tr>\n")
-  TableHeadItem(unmake.names(row.names(sigDiff)[i]))
-  for (j in 1:ncol(sigDiff)) {
-    TableItem(as.character(sigDiff[i, j]))
-  }
-  cata("</tr>\n")
+    cata("<tr>\n")
+    TableHeadItem(unmake.names(row.names(sigDiff)[i]))
+    for (j in 1:ncol(sigDiff)) {
+        TableItem(as.character(sigDiff[i, j]))
+    }
+    cata("</tr>\n")
 }
 cata("</table>")
 
 cata("<h4>Plots:</h4>\n")
 for (i in 1:nrow(linkData)) {
-  if (grepl(".pdf", linkData$Link[i])) {
-    HtmlLink(linkData$Link[i], linkData$Label[i])
+    if (grepl(".pdf", linkData$Link[i])) {
+        HtmlLink(linkData$Link[i], linkData$Label[i])
   }
 }
 
 cata("<h4>Tables:</h4>\n")
 for (i in 1:nrow(linkData)) {
-  if (grepl(".tsv", linkData$Link[i])) {
-    HtmlLink(linkData$Link[i], linkData$Label[i])
-  }
+    if (grepl(".tsv", linkData$Link[i])) {
+        HtmlLink(linkData$Link[i], linkData$Label[i])
+    }
 }
 
 if (wantRda) {
-  cata("<h4>R Data Object:</h4>\n")
-  for (i in 1:nrow(linkData)) {
-    if (grepl(".rda", linkData$Link[i])) {
-      HtmlLink(linkData$Link[i], linkData$Label[i])
+    cata("<h4>R Data Object:</h4>\n")
+    for (i in 1:nrow(linkData)) {
+        if (grepl(".RData", linkData$Link[i])) {
+            HtmlLink(linkData$Link[i], linkData$Label[i])
+        }
     }
-  }
 }
 
 cata("<p>Alt-click links to download file.</p>\n")
@@ -524,40 +687,60 @@
 
 cata("<h4>Additional Information</h4>\n")
 cata("<ul>\n")
-if (cpmReq!=0 && sampleReq!=0) {
-  tempStr <- paste("Genes without more than", cpmReq,
-                   "CPM in at least", sampleReq, "samples are insignificant",
-                   "and filtered out.")
-  ListItem(tempStr)
-  filterProp <- round(filteredCount/preFilterCount*100, digits=2)
-  tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
-                   "%) genes were filtered out for low expression.")
-  ListItem(tempStr)
-}
-ListItem(normOpt, " was the method used to normalise library sizes.")
-if (wantWeight) {
-  ListItem("Weights were applied to samples.")
-} else {
-  ListItem("Weights were not applied to samples.")
+
+if (filtCPM || filtSmpCount || filtTotCount) {
+    if (filtCPM) {
+    tempStr <- paste("Genes without more than", opt$cmpReq,
+                                     "CPM in at least", opt$sampleReq, "samples are insignificant",
+                                     "and filtered out.")
+    } else if (filtSmpCount) {
+        tempStr <- paste("Genes without more than", opt$cntReq,
+                                     "counts in at least", opt$sampleReq, "samples are insignificant",
+                                     "and filtered out.")
+    } else if (filtTotCount) {
+            tempStr <- paste("Genes without more than", opt$cntReq,
+                                     "counts, after summing counts for all samples, are insignificant",
+                                     "and filtered out.")
+    }
+
+    ListItem(tempStr)
+    filterProp <- round(filteredCount/preFilterCount*100, digits=2)
+    tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
+                                     "%) genes were filtered out for low expression.")
+    ListItem(tempStr)
 }
-if (pAdjOpt!="none") {
-  if (pAdjOpt=="BH" || pAdjOpt=="BY") {
-    tempStr <- paste0("MA-Plot highlighted genes are significant at FDR ",
-                      "of ", pValReq," and exhibit log2-fold-change of at ", 
-                      "least ", lfcReq, ".")
-    ListItem(tempStr)
-  } else if (pAdjOpt=="holm") {
-    tempStr <- paste0("MA-Plot highlighted genes are significant at adjusted ",
-                      "p-value of ", pValReq,"  by the Holm(1979) ",
-                      "method, and exhibit log2-fold-change of at least ", 
-                      lfcReq, ".")
-    ListItem(tempStr)
-  }
+ListItem(opt$normOpt, " was the method used to normalise library sizes.")
+if (wantTrend) {
+    ListItem("The limma-trend method was used.")
+} else {
+    ListItem("The limma-voom method was used.")
+}
+if (wantWeight) {
+    ListItem("Weights were applied to samples.")
 } else {
-  tempStr <- paste0("MA-Plot highlighted genes are significant at p-value ",
-                    "of ", pValReq," and exhibit log2-fold-change of at ", 
-                    "least ", lfcReq, ".")
-  ListItem(tempStr)
+    ListItem("Weights were not applied to samples.")
+}
+if (wantRobust) {
+    ListItem("eBayes was used with robust settings (robust=TRUE).")
+}
+if (opt$pAdjOpt!="none") {
+    if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") {
+        tempStr <- paste0("MA-Plot highlighted genes are significant at FDR ",
+                        "of ", opt$pValReq," and exhibit log2-fold-change of at ",
+                        "least ", opt$lfcReq, ".")
+        ListItem(tempStr)
+    } else if (opt$pAdjOpt=="holm") {
+        tempStr <- paste0("MA-Plot highlighted genes are significant at adjusted ",
+                        "p-value of ", opt$pValReq,"  by the Holm(1979) ",
+                        "method, and exhibit log2-fold-change of at least ",
+                        opt$lfcReq, ".")
+        ListItem(tempStr)
+    }
+  } else {
+        tempStr <- paste0("MA-Plot highlighted genes are significant at p-value ",
+                      "of ", opt$pValReq," and exhibit log2-fold-change of at ",
+                      "least ", opt$lfcReq, ".")
+        ListItem(tempStr)
 }
 cata("</ul>\n")
 
@@ -570,19 +753,18 @@
 TableHeadItem("SampleID")
 TableHeadItem(names(factors)[1]," (Primary Factor)")
 
-  if (ncol(factors) > 1) {
-
+if (ncol(factors) > 1) {
     for (i in names(factors)[2:length(names(factors))]) {
-      TableHeadItem(i)
+        TableHeadItem(i)
     }
     cata("</tr>\n")
-  }
+}
 
 for (i in 1:nrow(factors)) {
-  cata("<tr>\n")
-  TableHeadItem(row.names(factors)[i])
-  for (j in 1:ncol(factors)) {
-    TableItem(as.character(unmake.names(factors[i, j])))
+    cata("<tr>\n")
+    TableHeadItem(row.names(factors)[i])
+    for (j in 1:ncol(factors)) {
+        TableItem(as.character(unmake.names(factors[i, j])))
   }
   cata("</tr>\n")
 }
@@ -638,7 +820,7 @@
 cit[10] <- paste("Law CW, Chen Y, Shi W, and Smyth GK (2014). Voom:",
                 "precision weights unlock linear model analysis tools for",
                 "RNA-seq read counts. Genome Biology 15, R29.")
-cit[11] <- paste("Ritchie ME, Diyagama D, Neilson J, van Laar R,", 
+cit[11] <- paste("Ritchie ME, Diyagama D, Neilson J, van Laar R,",
                 "Dobrovic A, Holloway A and Smyth GK (2006).",
                 "Empirical array quality weights for microarray data.",
                 "BMC Bioinformatics 7, Article 261.")
@@ -667,9 +849,9 @@
 cata("<p>Please report problems or suggestions to: su.s@wehi.edu.au</p>\n")
 
 for (i in 1:nrow(linkData)) {
-  if (grepl("session_info", linkData$Link[i])) {
-    HtmlLink(linkData$Link[i], linkData$Label[i])
-  }
+    if (grepl("session_info", linkData$Link[i])) {
+        HtmlLink(linkData$Link[i], linkData$Label[i])
+    }
 }
 
 cata("<table border=\"0\">\n")
diff -r a330ddf43861 -r 38aab66ae5cb limma_voom.xml
--- a/limma_voom.xml	Thu Sep 07 05:27:27 2017 -0400
+++ b/limma_voom.xml	Wed Jan 31 12:45:42 2018 -0500
@@ -1,65 +1,100 @@
-<tool id="limma_voom" name="limma-voom" version="1.2.0">
+<tool id="limma_voom" name="limma" version="3.34.6.0">
     <description>
-        Differential expression with optional sample weights
+        Perform differential expression with limma-voom or limma-trend
     </description>
 
     <requirements>
-        <requirement type="package" version="3.16.5">bioconductor-edger</requirement>
-        <requirement type="package" version="3.30.13">bioconductor-limma</requirement>
-        <requirement type="package" version="1.4.29">r-statmod</requirement>
-        <requirement type="package" version="0.4.1">r-scales</requirement>
+        <requirement type="package" version="3.34.6">bioconductor-limma</requirement>
+        <requirement type="package" version="3.20.7">bioconductor-edger</requirement>
+        <requirement type="package" version="1.4.30">r-statmod</requirement>
+        <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>
     </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: ")
+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: ")
     ]]></version_command>
 
     <command detect_errors="exit_code"><![CDATA[
+#import json
 Rscript '$__tool_directory__/limma_voom.R'
-'$counts'
+
+-R '$outReport'
+-o '$outReport.files_path'
+
+#if $input.format=="files":
+
+    ## Adapted from DESeq2 wrapper
+    #set $temp_factor_names = list()
+    #for $fact in $input.rep_factor:
+        #set $temp_factor = list()
+        #for $g in $fact.rep_group:
+            #set $count_files = list()
+            #for $file in $g.countsFile:
+                $count_files.append(str($file))
+            #end for
+            $temp_factor.append( {str($g.groupName): $count_files} )
+        #end for
+
+        $temp_factor.reverse()
+        $temp_factor_names.append([str($fact.factorName), $temp_factor])
+    #end for
+    -j '#echo json.dumps(temp_factor_names)#'
+
+#elif $input.format=="matrix":
+    -m '$input.counts'
+    #if $input.fact.ffile=='yes':
+        -f '$input.fact.finfo'
+    #else:
+        -i '${ '|'.join( ['%s::%s' % ($x.factorName, $x.groupNames) for x in $input.fact.rep_factor] ) }'
+    #end if
+#end if
 
 #if $anno.annoOpt=='yes':
-  '$geneanno'
-#else:
-  None
+    -a '$anno.geneanno'
 #end if
 
-'$outReport'
-'$outReport.files_path'
-$rdaOption
-$normalisationOption
-$weightOption
-'$contrast'
+-C '${ ','.join( ['%s' % $x.contrast for x in $rep_contrast] ) }'
 
-#if $filterCPM.filterLowCPM=='yes':
-  '$filterCPM.cpmReq'
-  '$filterCPM.sampleReq'
-#else:
-  0
-  0
+#if $f.filt.filt_select == 'yes':
+    #if $f.filt.cformat.format_select == 'cpm':
+        -c '$f.filt.cformat.cpmReq'
+        -s '$f.filt.cformat.cpmSampleReq'
+    #elif $f.filt.cformat.format_select == 'counts':
+            -z '$f.filt.cformat.cntReq'
+        #if $f.filt.cformat.samples.count_select == 'total':
+            -y
+        #elif $f.filt.cformat.samples.count_select == 'sample':
+            -s '$f.filt.cformat.samples.cntSampleReq'
+        #end if
+    #end if
+#end if
+
+#if $out.normCounts:
+    -x
 #end if
 
-#if $testOpt.wantOpt=='yes':
-  '$testOpt.pAdjust'
-  '$testOpt.pVal'
-  '$testOpt.lfc'
-#else:
-  "BH"
-  0.05
-  0
+#if $out.rdaOption:
+    -r
 #end if
 
-$normCounts
+-l '$adv.lfc'
+-p '$adv.pVal'
+-d '$adv.pAdjust'
 
-#if $fact.ffile=='yes':
-    '$finfo'
-    'None'
-#else:
-    'None'
-    '$fact.pfactName::$fact.pfactLevel'
-    #for $sfact in $fact.sfactors:
-        '$sfact.sfactName::$sfact.sfactLevel'
-    #end for
+#if $deMethod.de_select == 'voom':
+    #if $deMethod.weightOption:
+        -w
+    #end if
+#elif $deMethod.de_select == 'trend':
+    -t $deMethod.prior_count
+#end if
+
+-n '$adv.normalisationOption'
+
+#if $adv.robOption:
+    -b
 #end if
 
 &&
@@ -70,11 +105,83 @@
     ]]></command>
 
     <inputs>
-        <param name="counts" type="data" format="tabular" label="Counts Data"/>
+
+        <!-- DE Method Option -->
+        <conditional name="deMethod">
+            <param name="de_select" type="select" label="Differential Expression Method" help="Select the limma-voom or limma-trend method. See Help section below for more information. Default: limma-voom">
+                <option value="voom" selected="True">limma-voom</option>
+                <option value="trend">limma-trend</option>
+            </param>
+            <when value="voom">
+                <param name="weightOption" type="boolean" truevalue="1" falsevalue="0" checked="false" label="Apply voom with sample quality weights?"
+                help="Apply weights if outliers are present (voomWithQualityWeights). Default: False.">
+                </param>
+            </when>
+            <when value="trend">
+                <param name="prior_count" type="float" min="0" value="3" label="Prior count" help="Average count to be added to each observation to avoid taking log of zero. Default: 3." />
+            </when>
+        </conditional>
+        <!-- Counts and Factors -->
+        <conditional name="input">
+            <param name="format" type="select" label="Count Files or Matrix?"
+                help="You can choose to input either separate count files (one per sample) or a single count matrix">
+                <option value="files">Separate Count Files</option>
+                <option value="matrix">Single Count Matrix</option>
+            </param>
 
+            <when value="files">
+                <repeat name="rep_factor" title="Factor" min="1">
+                    <param name="factorName" type="text" label="Name" help="Name of experiment factor of interest (e.g. Genotype). One factor must be entered and there must be two or more groups per factor. Optional additional factors (e.g. Batch) can be entered using the Insert Factor button below, see Help section for more information. NOTE: Please only use letters, numbers or underscores.">
+                    <sanitizer>
+                        <valid initial="string.letters,string.digits"><add value="_" /></valid>
+                    </sanitizer>
+                    </param>
+                    <repeat name="rep_group" title="Group" min="2" default="2">
+                        <param name="groupName" type="text" label="Name"
+                        help="Name of group that the counts files(s) belong to (e.g. WT or Mut). NOTE: Please only use letters, numbers or underscores (case sensitive).">
+                        <sanitizer>
+                            <valid initial="string.letters,string.digits"><add value="_" /></valid>
+                        </sanitizer>
+                        </param>
+                        <param name="countsFile" type="data" format="tabular" multiple="true" label="Counts file(s)"/>
+                    </repeat>
+                </repeat>
+            </when>
+
+            <when value="matrix">
+                <param name="counts" type="data" format="tabular" label="Count Matrix"/>
+
+                <conditional name="fact">
+                    <param name="ffile" type="select" label="Input factor information from file?"
+                        help="You can choose to input the factor and group information for the samples from a file or manually enter below.">
+                        <option value="no">No</option>
+                        <option value="yes">Yes</option>
+                    </param>
+                    <when value="yes">
+                        <param name="finfo" type="data" format="tabular" label="Factor File"/>
+                    </when>
+                    <when value="no" >
+                        <repeat name="rep_factor" title="Factor" min="1">
+                            <param name="factorName" type="text" label="Factor Name"
+                                help="Name of experiment factor of interest (e.g. Genotype). One factor must be entered and there must be two or more groups per factor. Additional factors (e.g. Batch) can be entered using the Insert Factor button below, see Help section below. NOTE: Please only use letters, numbers or underscores.">
+                                <validator type="empty_field" />
+                                <validator type="regex" message="Please only use letters, numbers or underscores">^[\w]+$</validator>
+                            </param>
+                            <param name="groupNames" type="text" label="Groups"
+                                help="Enter the group names for the samples separated with commas e.g. WT,WT,WT,Mut,Mut,Mut. The order of the names must match the order of the samples in the columns of the count matrix. NOTE: Please only use letters, numbers or underscores (case sensitive).">
+                                <validator type="empty_field" />
+                                <validator type="regex" message="Please only use letters, numbers or underscores, and separate levels by commas">^[\w,]+$</validator>
+                            </param>
+                        </repeat>
+                    </when>
+                </conditional>
+            </when>
+        </conditional>
+
+        <!-- Gene Annotations -->
         <conditional name="anno">
             <param name="annoOpt" type="select" label="Use Gene Annotations?"
-                help="If an annotation file is provided, annotations will be added to the table of differential expression results to provide descriptions for each gene.">
+                    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. See Help section below.">
                 <option value="no">No</option>
                 <option value="yes">Yes</option>
             </param>
@@ -84,102 +191,88 @@
             <when value="no" />
         </conditional>
 
-        <conditional name="fact">
-            <param name="ffile" type="select" label="Input Factor Information from file?"
-                help="You can choose to input the factor information from a file or manually enter below.">
-                <option value="no">No</option>
-                <option value="yes">Yes</option>
-            </param>
-            <when value="yes">
-                <param name="finfo" type="data" format="tabular" label="Factor Information"/>
-            </when>
-            <when value="no" >
-                <param name="pfactName" type="text" label="Primary Factor Name"
-                    help="Eg. Genotype NOTE: Please only use letters, numbers or underscores.">
-                    <validator type="empty_field" />
-                    <validator type="regex" message="Please only use letters, numbers or underscores">^[\w]+$</validator>
-                </param>
-                <param name="pfactLevel" type="text" label="Primary Factor Levels"
-                    help="Eg. WT,WT,WT,Mut,Mut,Mut NOTE: Please only use letters, numbers or underscores and ensure that the same levels are typed identically with cases matching.">
-                    <validator type="empty_field" />
-                    <validator type="regex" message="Please only use letters, numbers or underscores, and separate levels by commas">^[\w,]+$</validator>
-                </param>
-                <repeat name="sfactors" title="Secondary Factor" >
-                    <param name="sfactName" type="text" label="Secondary Factor Name" help="Eg. Batch">
-                        <validator type="empty_field" />
-                        <validator type="regex" message="Please only use letters, numbers or underscores">^[\w]+$</validator>
-                    </param>
-                    <param name="sfactLevel" type="text" label="Secondary Factor Levels"
-                        help="Eg. b1,b2,b3,b1,b2,b3 NOTE: Please only use letters, numbers or underscores and ensure that the same levels are typed identically with cases matching.">
-                        <validator type="empty_field" />
-                        <validator type="regex" message="Please only use letters, numbers or underscores">^[\w,]+$</validator>
-                    </param>
-                </repeat>
-            </when>
-        </conditional>
-
-        <param name="contrast" type="text" label="Contrasts of interest" help="Eg. Mut-WT,KD-Control">
-            <validator type="empty_field" />
-            <validator type="regex" message="Please only use letters, numbers or underscores">^[\w,-]+$</validator>
-        </param>
-
-        <conditional name="filterCPM">
-            <param name="filterLowCPM" type="select" label="Filter Low CPM?"
-                help="Treat genes with very low expression as unexpressed and filter out to speed up computation.">
-                <option value="yes" selected="True">Yes</option>
-                <option value="no">No</option>
+        <!-- Contrasts -->
+        <repeat name="rep_contrast" title="Contrast" min="1" default="1">
+            <param name="contrast" type="text" label="Contrast of Interest" help="Names of two groups to compare separated by a hyphen e.g. Mut-WT. If the order is Mut-WT the fold changes in the results will be up/down in Mut relative to WT. If you have more than one contrast enter each separately using the Insert Contrast button below. For more info, see Chapter 8 in the limma User's guide: https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf">
+                <validator type="empty_field" />
+                <validator type="regex" message="Please only use letters, numbers or underscores">^[\w-]+$</validator>
             </param>
-            <when value="yes">
-                <param name="cpmReq" type="float" value="0.5" min="0" label="Minimum CPM"/>
-                <param name="sampleReq" type="integer" value="1" min="0" label="Minimum Samples"
-                    help="Filter out all the genes that do not meet the minimum CPM in at least this many samples."/>
-            </when>
-            <when value="no"/>
-        </conditional>
-
-        <param name="weightOption" type="boolean" truevalue="yes" falsevalue="no" checked="false" label="Apply sample weights?"
-            help="Apply weights if outliers are present.">
-        </param>
-
-        <param name="normalisationOption" type="select" label="Normalisation Method">
-            <option value="TMM">TMM</option>
-            <option value="RLE">RLE</option>
-            <option value="upperquartile">Upperquartile</option>
-            <option value="none">None (Don't normalise)</option>
-        </param>
-
-        <param name="normCounts" type="boolean" truevalue="yes" falsevalue="no" checked="false"
-            label="Output normalised counts table?"
-            help="Output a file containing the normalised counts, these are in log2 counts per million (logCPM).">
-        </param>
+        </repeat>
 
-        <param name="rdaOption" type="boolean" truevalue="yes" falsevalue="no" checked="false"
-            label="Output RData?"
-            help="Output all the data used by R to construct the plots and tables, can be loaded into R. A link to the RData file will be provided in the HTML report.">
-        </param>
+        <!-- Filter Options -->
+        <section name="f" expanded="false" title="Filter Low Counts">
+            <conditional name="filt">
+                <param name="filt_select" type="select" label="Filter lowly expressed genes?" help="Treat genes with very low expression as unexpressed and filter out. See the Filter Low Counts section below for more information. Default: No">
+                    <option value="no" selected="true">No</option>
+                    <option value="yes">Yes</option>
+                </param>
+                <when value="yes">
+                    <conditional name="cformat">
+                        <param name="format_select" type="select" label="Filter on CPM or Count values?" help="It is slightly better to base the filtering on count-per-million (CPM) rather than the raw count values so as to avoid favoring genes expressed in samples sequenced to a higher depth. ">
+                            <option value="cpm">CPM</option>
+                            <option value="counts">Counts</option>
+                        </param>
+                        <when value="cpm">
+                            <param name="cpmReq" type="float" value="1" min="0" label="Minimum CPM" help="Treat genes with CPM below this value as unexpressed and filter out. See the Filter Low Counts section below for more information."/>
+                             <param name="cpmSampleReq" type="integer" value="0" min="0" label="Minimum Samples"
+                             help="Filter out all genes that do not meet the Minimum CPM in at least this many samples. See the Filter Low Counts section below for more information."/>
+                        </when>
+                        <when value="counts">
+                            <param name="cntReq" type="integer" value="0" min="0" label="Minimum Count" help="Filter out all genes that do not meet this minimum count. You can choose below to apply this filter to the total count for all samples or specify the number of samples under Minimum Samples. See the Filter Low Counts section below for more information." />
+                            <conditional name="samples">
+                                <param name="count_select" type="select" label="Filter on Total Count or per Sample Count values?" >
+                                    <option value="total">Total</option>
+                                    <option value="sample">Sample</option>
+                                </param>
+                                <when value="total">
+                                    <param name="totReq" type="boolean" truevalue="1" falsevalue="0" checked="false" label="Filter on Total Count" help="Apply the Minimum Count filter to genes after summing counts for all samples. See the Filter Low Counts section below for more information." />
+                                </when>
+                                <when value="sample">
+                                    <param name="cntSampleReq" type="integer" value="0" min="0" label="Minimum Samples"
+                                    help="Filter out all genes that do not meet the Minimum Count in at least this many samples. See the Filter Low Counts section below for more information."/>
+                                </when>
+                            </conditional>
+                        </when>
+                    </conditional>
+                </when>
+                <when value="no" />
+            </conditional>
+        </section>
 
-        <conditional name="testOpt">
-            <param name="wantOpt" type="select" label="Use Advanced Testing Options?"
-                help="Enable choices for p-value adjustment method, p-value threshold and log2-fold-change threshold.">
-                <option value="no" selected="True">No</option>
-                <option value="yes">Yes</option>
+        <!-- Output Options -->
+        <section name="out" expanded="false" title="Output Options">
+            <param name="normCounts" type="boolean" truevalue="1" falsevalue="0" checked="false"
+                label="Output Normalised Counts Table?"
+                help="Output a file containing the normalised counts, these are in log2 counts per million (logCPM). Default: No">
             </param>
-            <when value="yes">
-                <param name="pAdjust" type="select" label="P-Value Adjustment Method.">
-                    <option value="BH">Benjamini and Hochberg (1995)</option>
-                    <option value="BY">Benjamini and Yekutieli (2001)</option>
-                    <option value="holm">Holm (1979)</option>
-                    <option value="none">None</option>
-                </param>
-                <param name="pVal" type="float" value="0.05" min="0" max="1"
-                    label="Adjusted Threshold"
-                    help="Genes below this threshold are considered significant and highlighted in the MA plot. If either BH(1995) or BY(2001) were selected then this value is a false-discovery-rate control. If Holm(1979) was selected then this is an adjusted p-value for family-wise error rate."/>
-                <param name="lfc" type="float" value="0" min="0"
-                    label="Minimum log2-fold-change Required"
-                    help="Genes above this threshold and below the p-value threshold are considered significant and highlighted in the MA plot."/>
-            </when>
-            <when value="no"/>
-        </conditional>
+            <param name="rdaOption" type="boolean" truevalue="1" falsevalue="0" checked="false"
+                label="Output RData file?"
+                help="Output all the data used by R to construct the plots and tables, can be loaded into R. A link to the RData file will be provided in the HTML report. Default: No">
+            </param>
+        </section>
+
+        <!-- Advanced Options -->
+        <section name="adv" expanded="false" title="Advanced Options">
+            <param name="lfc" type="float" value="0" min="0"
+                label="Minimum Log2 Fold Change"
+                help="Genes above this threshold and below the p-value threshold are considered significant and highlighted in the MD plot. Default: 0."/>
+            <param name="pVal" type="float" value="0.05" min="0" max="1"
+                label="P-Value Adjusted Threshold"
+                help="Genes below this threshold are considered significant and highlighted in the MD plot. If either BH(1995) or BY(2001) are selected then this value is a false-discovery-rate control. If Holm(1979) is selected then this is an adjusted p-value for family-wise error rate. Default: 0.05."/>
+            <param name="pAdjust" type="select" label="P-Value Adjustment Method" help="Default: BH">
+                <option value="BH" selected="true">Benjamini and Hochberg (1995)</option>
+                <option value="BY">Benjamini and Yekutieli (2001)</option>
+                <option value="holm">Holm (1979)</option>
+                <option value="none">None</option>
+            </param>
+            <param name="normalisationOption" type="select" label="Normalisation Method" help="Default: TMM">
+                <option value="TMM" selected="true">TMM</option>
+                <option value="RLE">RLE</option>
+                <option value="upperquartile">Upperquartile</option>
+                <option value="none">None (Don't normalise)</option>
+            </param>
+            <param name="robOption" type="boolean" truevalue="1" falsevalue="0" checked="true" label="Use Robust Settings?" help="Using robust settings is usually recommended to protect against outlier genes. Default: Yes" />
+        </section>
 
     </inputs>
 
@@ -191,11 +284,20 @@
     </outputs>
 
     <tests>
+        <!-- Ensure report is output -->
         <test>
+            <param name="format" value="matrix" />
             <param name="counts" value="matrix.txt" />
-            <param name="pfactName" value="Genotype" />
-            <param name="pfactLevel" value="WT,WT,WT,Mut,Mut,Mut" />
-            <param name="contrast" value="Mut-WT,WT-Mut" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="WT-Mut" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
             <output_collection name="outTables" count="2">
                 <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT.tsv" />
@@ -203,29 +305,41 @@
             </output_collection>
             <output name="outReport" >
                 <assert_contents>
-                    <has_text text="Limma-voom Analysis Output" />
+                    <has_text text="Limma Analysis Output" />
                     <not_has_text text="RData" />
                 </assert_contents>
             </output>
         </test>
+        <!-- Ensure annotation file input works -->
         <test>
+            <param name="format" value="matrix" />
             <param name="annoOpt" value="yes" />
             <param name="geneanno" value="anno.txt" />
             <param name="counts" value="matrix.txt" />
-            <param name="pfactName" value="Genotype" />
-            <param name="pfactLevel" value="WT,WT,WT,Mut,Mut,Mut" />
-            <param name="contrast" value="Mut-WT" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
-            <output_collection name="outTables" >
-                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WTanno.tsv" />
+            <output_collection name="outTables" count="1">
+                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT_anno.tsv" />
             </output_collection>
         </test>
+        <!-- Ensure RData file can be output -->
         <test>
-            <param name="rdaOption" value="yes" />
+            <param name="format" value="matrix" />
+            <param name="rdaOption" value="true" />
             <param name="counts" value="matrix.txt" />
-            <param name="pfactName" value="Genotype" />
-            <param name="pfactLevel" value="WT,WT,WT,Mut,Mut,Mut" />
-            <param name="contrast" value="Mut-WT" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
             <output name="outReport" >
                 <assert_contents>
@@ -233,42 +347,125 @@
                 </assert_contents>
             </output>
         </test>
+        <!-- Ensure secondary factors work -->
         <test>
+            <param name="format" value="matrix" />
             <param name="counts" value="matrix.txt" />
-            <param name="pfactName" value="Genotype"/>
-            <param name="pfactLevel" value="WT,WT,WT,Mut,Mut,Mut"/>
-            <repeat name="sfactors">
-                <param name="sfactName" value="Batch"/>
-                <param name="sfactLevel" value="b1,b2,b3,b1,b2,b3"/>
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
             </repeat>
-            <param name="contrast" value="Mut-WT" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Batch"/>
+                <param name="groupNames" value="b1,b2,b3,b1,b2,b3"/>
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
-            <output_collection name="outTables" >
-                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WTmultifact.tsv" />
+            <output_collection name="outTables" count="1" >
+                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT_2fact.tsv" />
             </output_collection>
         </test>
+        <!-- Ensure factors file input works -->
         <test>
+            <param name="format" value="matrix" />
             <param name="ffile" value="yes" />
             <param name="finfo" value="factorinfo.txt" />
             <param name="counts" value="matrix.txt" />
-            <param name="contrast" value="Mut-WT" />
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
-            <output_collection name="outTables" >
-                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WTmultifact.tsv" />
+            <output_collection name="outTables" count="1">
+                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT_2fact.tsv" />
             </output_collection>
         </test>
+        <!-- Ensure normalised counts file output works-->
         <test>
-            <param name="normCounts" value="yes" />
+            <param name="format" value="matrix" />
+            <param name="normCounts" value="true" />
             <param name="counts" value="matrix.txt" />
-            <param name="pfactName" value="Genotype" />
-            <param name="pfactLevel" value="WT,WT,WT,Mut,Mut,Mut" />
-            <param name="contrast" value="Mut-WT" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
             <param name="normalisationOption" value="TMM" />
             <output_collection name="outTables" count="2">
                 <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT.tsv" />
                 <element name="limma-voom_normcounts" ftype="tabular" file="limma-voom_normcounts.tsv" />
             </output_collection>
         </test>
+        <!-- Ensure multiple counts files input works -->
+        <test>
+            <param name="format" value="files" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <repeat name="rep_group">
+                    <param name="groupName" value="WT"/>
+                    <param name="countsFile" value="WT1.counts,WT2.counts,WT3.counts"/>
+                </repeat>
+                <repeat name="rep_group">
+                    <param name="groupName" value="Mut"/>
+                    <param name="countsFile" value="Mut1.counts,Mut2.counts,Mut3.counts"/>
+                </repeat>
+            </repeat>
+            <repeat name="rep_factor">
+                <param name="factorName" value="Batch"/>
+                <repeat name="rep_group">
+                    <param name="groupName" value="b1"/>
+                    <param name="countsFile" value="WT1.counts,Mut1.counts"/>
+                </repeat>
+                <repeat name="rep_group">
+                    <param name="groupName" value="b2"/>
+                    <param name="countsFile" value="WT2.counts,Mut2.counts"/>
+                </repeat>
+                <repeat name="rep_group">
+                    <param name="groupName" value="b3"/>
+                    <param name="countsFile" value="WT3.counts,Mut3.counts"/>
+                </repeat>
+            </repeat>
+            <param name="annoOpt" value="yes" />
+            <param name="geneanno" value="anno.txt" />
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="WT-Mut" />
+            </repeat>
+            <param name="normCounts" value="true" />
+            <output_collection name="outTables" count="3">
+                <element name="limma-voom_Mut-WT" ftype="tabular" file="limma-voom_Mut-WT_2fact_anno.tsv" />
+                <element name="limma-voom_WT-Mut" ftype="tabular" file="limma-voom_WT-Mut_2fact_anno.tsv" />
+                <element name="limma-voom_normcounts" ftype="tabular" file="limma-voom_normcounts_anno.tsv" />
+            </output_collection>
+        </test>
+        <!-- Ensure limma-trend option works -->
+        <test>
+            <param name="format" value="matrix" />
+            <param name="counts" value="matrix.txt" />
+            <repeat name="rep_factor">
+                <param name="factorName" value="Genotype"/>
+                <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT" />
+            </repeat>
+            <repeat name="rep_contrast">
+                <param name="contrast" value="Mut-WT" />
+            </repeat>
+            <param name="normalisationOption" value="TMM" />
+            <param name="de_select" value="trend" />
+            <param name="rdaOption" value="true" />
+            <output name="outReport" >
+                <assert_contents>
+                    <has_text text="The limma-trend method was used" />
+                </assert_contents>
+            </output>
+            <output_collection name="outTables" count="1">
+                <element name="limma-trend_Mut-WT" ftype="tabular" file="limma-trend_Mut-WT.tsv" />
+            </output_collection>
+        </test>
     </tests>
 
     <help><![CDATA[
@@ -276,21 +473,38 @@
 
 **What it does**
 
-Given a matrix of counts (e.g. from featureCounts) and optional information about the genes, this tool
-produces plots and tables useful in the analysis of differential gene
-expression.
+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.
+
+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.
 
 -----
 
 **Inputs**
 
+**Differential Expression Method:**
+Option to use the limma-voom or limma-trend approach for differential expression. The default is limma-voom.
+If the sequencing depth is reasonably consistent across the RNA samples, then the simplest and most
+robust approach to differential expression is to use limma-trend. This approach will usually work well if the
+ratio of the largest library size to the smallest is not more than about 3-fold. When the library sizes are quite variable between samples, then the voom approach is theoretically more powerful than limma-trend. For more information see the excellent `limma User's Guide`_.
+
 **Counts Data:**
-A matrix of counts, with rows corresponding to genes
-and columns corresponding to counts for the samples.
-Values must be tab separated, with the first row containing the sample/column
-labels and the first column containing the row/gene labels.
+
+The counts data can either be input as separate counts files (one sample per file) or a single count matrix (one sample per column). The rows correspond to genes, and columns correspond to the counts for the samples. Values must be tab separated, with the first row containing the sample/column labels and the first column containing the row/gene labels. Gene identifiers can be of any type but must be unique and not repeated within a counts file.
+
+Example - **Separate Count Files**:
 
-Example:
+    ========== =======
+    **GeneID** **WT1**
+    ---------- -------
+    11287      1699
+    11298      1905
+    11302      6
+    11303      2099
+    11304      356
+    11305      2528
+    ========== =======
+
+Example - **Single Count Matrix**:
 
     ========== ======= ======= ======= ======== ======== ========
     **GeneID** **WT1** **WT2** **WT3** **Mut1** **Mut2** **Mut3**
@@ -322,7 +536,7 @@
     ==========  ==========  ===================================================
 
 **Factor Information:**
-Enter Factor Names and Levels in the tool form or provide a tab-separated file that has the samples in the same order as listed in the columns of the counts matrix. The second column should contain the Primary Factor levels (e.g. Genotype) with optional additional columns for any Secondary Factors (e.g. Batch).
+Enter factor names and groups in the tool form, or provide a tab-separated file that has the samples in the same order as listed in the columns of the counts matrix. The second column should contain the primary factor levels (e.g. WT, Mut) with optional additional columns for any secondary factors.
 
 Example:
 
@@ -337,88 +551,145 @@
     Mut3       Mut          b3
     ========== ============ =========
 
-**Primary Factor Name:** The name of the primary factor being investigated e.g. Genotype. One primary factor must be entered and spaces must not be used.
+*Factor Name:* The name of the experimental factor being investigated e.g. Genotype, Treatment. One factor must be entered and spaces must not be used. Optionally, additional factors can be included, these are variables that might influence your experiment e.g. Batch, Gender, Subject. If additional factors are entered, edgeR will fit an additive linear model.
+
+*Groups:* The names of the groups for the factor. These must be entered in the same order as the samples (to which the groups correspond) are listed in the columns of the counts matrix. Spaces must not be used and if entered into the tool form above, the values should be separated by commas.
+
 
-**Primary Factor Levels:** The levels of the primary factor of interest, these must be entered in the same order as the samples to which the levels correspond, as listed in the columns of the counts matrix. Spaces must not be used and if entered in the tool form the values should be separated by commas.
+**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.
 
-**Secondary Factor Name:** Optionally, one or more secondary factors can be included. These are variables that might influence your experiment e.g. Batch, Gender. Spaces must not be used.
+Example:
 
-**Secondary Factor Levels:** The levels of the secondary factor of interest, these must be entered in the same order as the samples to which the levels correspond, as listed in the columns of the counts matrix. Spaces must not be used and if entered in the tool form the values should be separated by commas.
-
+    ==========  ==========  ===================================================
+    **GeneID**  **Symbol**  **GeneName**
+    ----------  ----------  ---------------------------------------------------
+    1287        Pzp         pregnancy zone protein
+    1298        Aanat       arylalkylamine N-acetyltransferase
+    1302        Aatk        apoptosis-associated tyrosine kinase
+    1303        Abca1       ATP-binding cassette, sub-family A (ABC1), member 1
+    1304        Abca4       ATP-binding cassette, sub-family A (ABC1), member 4
+    1305        Abca2       ATP-binding cassette, sub-family A (ABC1), member 2
+    ==========  ==========  ===================================================
 
 **Contrasts of Interest:**
 The contrasts you wish to make between levels.
 A common contrast would be a simple difference between two levels: "Mut-WT"
 represents the difference between the mutant and wild type genotypes.
-Multiple contrasts should be separated by commas and spaces must not be used.
+Multiple contrasts must be entered separately using the Insert Contrast button, spaces must not be used.
 
-**Filter Low CPM:**
-Option to ignore the genes that do not show significant levels of
-expression, this filtering is dependent on two criteria:
+**Filter Low Counts:**
+Genes with very low counts across all libraries provide little evidence for differential expression.
+In the biological point of view, a gene must be expressed at some minimal level before
+it is likely to be translated into a protein or to be biologically important. In addition, the
+pronounced discreteness of these counts interferes with some of the statistical approximations
+that are used later in the pipeline. These genes should be filtered out prior to further
+analysis.
+As a rule of thumb, genes are dropped if they can’t possibly be expressed in all the samples
+for any of the conditions. Users can set their own definition of genes being expressed. Usually
+a gene is required to have a count of 5-10 in a library to be considered expressed in that
+library. Users should also filter with count-per-million (CPM) rather than filtering on the
+counts directly, as the latter does not account for differences in library sizes between samples.
 
-    * **Minimum CPM:** This is the counts per million that a gene must have in at
-      least some specified number of samples.
+Option to ignore the genes that do not show significant levels of
+expression, this filtering is dependent on two criteria: CPM/count and number of samples. You can specify to filter on CPM (Minimum CPM) or count (Minimum Count) values:
+
+    * **Minimum CPM:** This is the minimum count per million that a gene must have in at
+      least the number of samples specified under Minimum Samples.
 
-    * **Minumum Samples:** This is the number of samples in which the CPM
-      requirement must be met in order for that gene to be acknowledged.
+    * **Minimum Count:** This is the minimum count that a gene must have. It can be combined with either Filter
+      on Total Count or Minimum Samples.
 
-Only genes that exhibit a CPM greater than the required amount in at least the
-number of samples specified will be used for analysis. Care should be taken to
+    * **Filter on Total Count:** This can be used with the Minimum Count filter to keep genes
+      with a minimum total read count.
+
+    * **Minimum Samples:** This is the number of samples in which the Minimum CPM/Count
+      requirement must be met in order for that gene to be kept.
+
+If the Minimum Samples filter is applied, only genes that exhibit a CPM/count greater than the required amount in at least the number of samples specified will be used for analysis. Care should be taken to
 ensure that the sample requirement is appropriate. In the case of an experiment
 with two experimental groups each with two members, if there is a change from
-insignificant cpm to significant cpm but the sample requirement is set to 3,
+insignificant CPM/count to significant CPM/count but the sample requirement is set to 3,
 then this will cause that gene to fail the criteria. When in doubt simply do not
-filter.
+filter or consult the `limma User's Guide`_ for filtering recommendations.
 
-**Normalisation Method:**
-Option for using different methods to rescale the raw library
-size. For more information, see calcNormFactor section in the edgeR_ user's
-manual.
+**Advanced Options:**
 
-**Apply Sample Weights:**
-Option to downweight outlier samples such that their information is still
-used in the statistical analysis but their impact is reduced. Use this
-whenever significant outliers are present. The MDS plotting tool in this package
-is useful for identifying outliers. For more information on this option see Liu et al. (2015).
-
-**Use Advanced Testing Options?:**
 By default error rate for multiple testing is controlled using Benjamini and
 Hochberg's false discovery rate control at a threshold value of 0.05. However
 there are options to change this to custom values.
 
+    * **Minimum log2-fold-change Required:**
+      In addition to meeting the requirement for the adjusted statistic for
+      multiple testing, the observation must have an absolute log2-fold-change
+      greater than this threshold to be considered significant, thus highlighted
+      in the MD plot.
+
+    * **Adjusted Threshold:**
+      Set the threshold for the resulting value of the multiple testing control
+      method. Only observations whose statistic falls below this value is
+      considered significant, thus highlighted in the MD plot.
+
     * **P-Value Adjustment Method:**
       Change the multiple testing control method, the options are BH(1995) and
       BY(2001) which are both false discovery rate controls. There is also
       Holm(1979) which is a method for family-wise error rate control.
 
-    * **Adjusted Threshold:**
-      Set the threshold for the resulting value of the multiple testing control
-      method. Only observations whose statistic falls below this value is
-      considered significant, thus highlighted in the MA plot.
+**Normalisation Method:**
+The most obvious technical factor that affects the read counts, other than gene expression
+levels, is the sequencing depth of each RNA sample. edgeR adjusts any differential expression
+analysis for varying sequencing depths as represented by differing library sizes. This is
+part of the basic modeling procedure and flows automatically into fold-change or p-value
+calculations. It is always present, and doesn’t require any user intervention.
+The second most important technical influence on differential expression is one that is less
+obvious. RNA-seq provides a measure of the relative abundance of each gene in each RNA
+sample, but does not provide any measure of the total RNA output on a per-cell basis.
+This commonly becomes important when a small number of genes are very highly expressed
+in one sample, but not in another. The highly expressed genes can consume a substantial
+proportion of the total library size, causing the remaining genes to be under-sampled in that
+sample. Unless this RNA composition effect is adjusted for, the remaining genes may falsely
+appear to be down-regulated in that sample . The edgeR `calcNormFactors` function normalizes for RNA composition by finding a set of scaling factors for the library sizes that minimize the log-fold changes between the samples for most genes. The default method for computing these scale factors uses a trimmed mean of M values (TMM) between each pair of samples. We call the product of the original library size and the scaling factor the *effective library size*. The effective library size replaces the original library size in all downsteam analyses. TMM is the recommended method for most RNA-Seq data where the majority (more than half) of the genes are believed not differentially expressed between any pair of the samples. You can change the normalisation method under **Advanced Options** above. For more information, see the `calcNormFactors` section in the `edgeR User's Guide`_.
 
-    * **Minimum log2-fold-change Required:**
-      In addition to meeting the requirement for the adjusted statistic for
-      multiple testing, the observation must have an absolute log2-fold-change
-      greater than this threshold to be considered significant, thus highlighted
-      in the MA plot.
+**Robust Settings**
+Option to use robust settings with eBayes, used by both liamm-voom and limma-trend. Using robust settings is usually recommended to protect against outlier genes, for more information see the `limma User's Guide`_. This is turned on by default.
+
+**Prior Count:**
+If the limma-trend method is used, a count (`prior.count`) is added to all counts to avoid taking a log of zero, and damp down the variances of logarithms of low counts. A default of 3 is used, as recommended in the `limma User's Guide`_.
+
+**Apply Sample Weights:**
+If the limma-voom method is used, an option is available to downweight outlier samples, such that their information is still
+used in the statistical analysis but their impact is reduced. Use this
+whenever significant outliers are present. The MDS plotting tool in this package
+is useful for identifying outliers. For more information on this option see Liu et al. (2015).
+
 
 **Outputs**
 
-This tool outputs a table of differentially expressed genes for each contrast of interest and a HTML report with plots and additional information. Optionally you can choose to output the normalised counts table and the RData file.
+This tool outputs
+
+    * a table of differentially expressed genes for each contrast of interest
+    * a HTML report with plots and additional information
+
+Optionally, under **Output Options** you can choose to output
+
+    * a normalised counts table
+    * an RData file
 
 -----
 
 **Citations:**
 
-.. class:: infomark
+Please try to cite the appropriate articles when you publish results obtained using software, as such citation is the main means by which the authors receive credit for their work.
 
 limma
 
 Please cite the paper below for the limma software itself. Please also try
 to cite the appropriate methodology articles that describe the statistical
 methods implemented in limma, depending on which limma functions you are
-using.  The methodology articles are listed in Section 2.1 of the limma
-User's Guide.
+using.  The methodology articles are listed in Section 2.1 of the `limma
+User's Guide`_.
 
     * Smyth GK (2005). Limma: linear models for microarray data. In:
       'Bioinformatics and Computational Biology Solutions using R and
@@ -435,13 +706,12 @@
       A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights
       for microarray data. BMC Bioinformatics 7, Article 261.
 
-.. class:: infomark
 
 edgeR
 
 Please cite the first paper for the software itself and the other papers for
 the various original statistical methods implemented in edgeR.  See
-Section 1.2 in the User's Guide for more detail.
+Section 1.2 in the `edgeR User's Guide`_ for more detail.
 
     * Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor
       package for differential expression analysis of digital gene expression
@@ -460,10 +730,14 @@
 
 Please report problems or suggestions to: su.s@wehi.edu.au
 
+.. _limma: http://www.bioconductor.org/packages/release/bioc/html/limma.html
+.. _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
-.. _limma: http://www.bioconductor.org/packages/release/bioc/html/limma.html
+.. _edgeR User's Guide: https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
     ]]></help>
     <citations>
+        <citation type="doi">10.1186/gb-2014-15-2-r29</citation>
         <citation type="doi">10.1093/nar/gkv412</citation>
     </citations>
-</tool>
+</tool>
\ No newline at end of file
diff -r a330ddf43861 -r 38aab66ae5cb test-data/Mut1.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/Mut1.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	Mut1
+11287	1463
+11298	1345
+11302	5
+11303	1574
+11304	361
+11305	1762
diff -r a330ddf43861 -r 38aab66ae5cb test-data/Mut2.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/Mut2.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	Mut2
+11287	1441
+11298	1291
+11302	6
+11303	1519
+11304	397
+11305	1942
diff -r a330ddf43861 -r 38aab66ae5cb test-data/Mut3.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/Mut3.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	Mut3
+11287	1495
+11298	1346
+11302	5
+11303	1654
+11304	346
+11305	2027
diff -r a330ddf43861 -r 38aab66ae5cb test-data/WT1.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/WT1.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	WT1
+11287	1699
+11298	1905
+11302	6
+11303	2099
+11304	356
+11305	2528
diff -r a330ddf43861 -r 38aab66ae5cb test-data/WT2.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/WT2.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	WT2
+11287	1528
+11298	1744
+11302	8
+11303	1974
+11304	312
+11305	2438
diff -r a330ddf43861 -r 38aab66ae5cb test-data/WT3.counts
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/WT3.counts	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+GeneID	WT3
+11287	1601
+11298	1834
+11302	7
+11303	2100
+11304	337
+11305	2493
diff -r a330ddf43861 -r 38aab66ae5cb test-data/factorinfo.txt
--- a/test-data/factorinfo.txt	Thu Sep 07 05:27:27 2017 -0400
+++ b/test-data/factorinfo.txt	Wed Jan 31 12:45:42 2018 -0500
@@ -1,7 +1,7 @@
 Sample	Genotype	Batch
+Mut1	Mut	b1
+Mut2	Mut	b2
+Mut3	Mut	b3
 WT1	WT	b1
 WT2	WT	b2
 WT3	WT	b3
-Mut1	Mut	b1
-Mut2	Mut	b2
-Mut3	Mut	b3
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-trend_Mut-WT.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-trend_Mut-WT.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
+"11304"	0.454006866294872	15.5347733798565	5.73783569533779	6.51380077190898e-06	3.90828046314539e-05	9.50615253604686
+"11287"	0.189211879829905	17.6560194022488	4.52571052740792	0.000138730279836171	0.000416190839508512	3.37443017677096
+"11298"	-0.137906115120717	17.6760627799956	-3.28370943058883	0.00313397342215802	0.00626794684431604	-1.40612492217495
+"11303"	-0.0565048762006306	17.887811868748	-1.22181580605621	0.233642901495507	0.35046435224326	-5.9848993971045
+"11305"	-0.0597744562444475	18.1592280697576	-1.04836609029276	0.304913314924737	0.365895977909685	-6.17897539114193
+"11302"	-0.0455235218343741	10.2501491061291	-0.373618198821297	0.711968519770644	0.711968519770644	-6.65188039014484
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WT.tsv
--- a/test-data/limma-voom_Mut-WT.tsv	Thu Sep 07 05:27:27 2017 -0400
+++ b/test-data/limma-voom_Mut-WT.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -1,7 +1,7 @@
 "GeneID"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
-"11304"	0.457332061341026	15.5254133001226	6.50459574633681	9.98720685006039e-07	5.99232411003624e-06	14.0741948485896
-"11287"	0.190749727701785	17.6546448244617	5.09535410066402	3.26518807654125e-05	9.79556422962375e-05	5.46773893802392
-"11298"	-0.138014418336201	17.6747285193431	-3.33168485842331	0.00278753263633162	0.00557506527266324	-1.84301342041449
-"11303"	-0.0558958943606989	17.886791401216	-1.30108531275576	0.205582481502297	0.254491025872973	-6.4924124057801
-"11305"	-0.0606991650996633	18.1585474109909	-1.28203791127299	0.212075854894144	0.254491025872973	-6.42090197700503
-"11302"	-0.0350239682204432	9.78883119065989	-0.236945963165269	0.814709535394087	0.814709535394087	-6.09497670655944
+"11304"	0.457332061341022	15.5254133001226	6.50459574633635	9.98720685007144e-07	5.99232411004286e-06	14.0741948485867
+"11287"	0.190749727701785	17.6546448244617	5.09535410066427	3.26518807653921e-05	9.79556422961764e-05	5.46773893802514
+"11298"	-0.138014418336201	17.6747285193431	-3.3316848584234	0.00278753263633103	0.00557506527266206	-1.84301342041422
+"11303"	-0.0558958943607024	17.886791401216	-1.30108531275582	0.205582481502277	0.254491025872966	-6.49241240578
+"11305"	-0.0606991650996669	18.1585474109909	-1.28203791127301	0.212075854894138	0.254491025872966	-6.42090197700496
+"11302"	-0.035023968220445	9.78883119065989	-0.236945963165271	0.814709535394086	0.814709535394086	-6.09497670655939
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WT_2fact.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-voom_Mut-WT_2fact.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"GeneID"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
+"11287"	0.190521683937302	17.6546448244617	14.4538942795641	5.93483727833391e-09	3.56090236700035e-08	96.5864491179047
+"11298"	-0.13987741320075	17.6747285193431	-7.7190169464239	5.41139783114199e-06	1.6234193493426e-05	22.3495333961476
+"11304"	0.459011254726059	15.5254133001226	5.6408319840902	0.00010884906092378	0.000217698121847559	8.76657226635715
+"11303"	-0.0641599901594248	17.886791401216	-2.9900818153928	0.0112725655419901	0.0169088483129851	-2.70089035288878
+"11305"	-0.0651479753016169	18.1585474109909	-2.28935282063866	0.0409794711446066	0.0491753653735279	-4.2713352660451
+"11302"	-0.0358817134644287	9.78883119065989	-0.439436789626843	0.668154169055813	0.668154169055813	-5.75838315483931
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WT_2fact_anno.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-voom_Mut-WT_2fact_anno.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"EntrezID"	"Symbol"	"GeneName"	"Chr"	"Length"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
+11287	"Pzp"	"pregnancy zone protein"	6	4681	0.190521683937302	17.6546448244617	15.1152647604854	3.5605743515492e-09	2.13634461092952e-08	106.319792304784
+11298	"Aanat"	"arylalkylamine N-acetyltransferase"	11	1455	-0.139877413200757	17.6747285193431	-8.49164953574146	2.03124167559177e-06	6.0937250267753e-06	28.5167891813121
+11304	"Abca4"	"ATP-binding cassette, sub-family A (ABC1), member 4"	3	7248	0.459011254726057	15.5254133001226	5.60316118539212	0.000115556485215252	0.000231112970430504	8.56249981126828
+11303	"Abca1"	"ATP-binding cassette, sub-family A (ABC1), member 1"	4	10260	-0.0641599901594141	17.886791401216	-2.92798536696768	0.0126512435710263	0.0189768653565394	-2.86386681629442
+11305	"Abca2"	"ATP-binding cassette, sub-family A (ABC1), member 2"	2	8061	-0.0651479753016169	18.1585474109909	-2.03911882030581	0.0640903140266962	0.0769083768320354	-4.69855556350315
+11302	"Aatk"	"apoptosis-associated tyrosine kinase"	11	5743	-0.0358817134644287	9.78883119065989	-0.43940476739557	0.668176726641997	0.668176726641997	-5.75841999159995
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WT_anno.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-voom_Mut-WT_anno.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"EntrezID"	"Symbol"	"GeneName"	"Chr"	"Length"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
+11304	"Abca4"	"ATP-binding cassette, sub-family A (ABC1), member 4"	3	7248	0.457332061341022	15.5254133001226	6.50459574633635	9.98720685007144e-07	5.99232411004286e-06	14.0741948485867
+11287	"Pzp"	"pregnancy zone protein"	6	4681	0.190749727701785	17.6546448244617	5.09535410066427	3.26518807653921e-05	9.79556422961764e-05	5.46773893802514
+11298	"Aanat"	"arylalkylamine N-acetyltransferase"	11	1455	-0.138014418336201	17.6747285193431	-3.3316848584234	0.00278753263633103	0.00557506527266206	-1.84301342041422
+11303	"Abca1"	"ATP-binding cassette, sub-family A (ABC1), member 1"	4	10260	-0.0558958943607024	17.886791401216	-1.30108531275582	0.205582481502277	0.254491025872966	-6.49241240578
+11305	"Abca2"	"ATP-binding cassette, sub-family A (ABC1), member 2"	2	8061	-0.0606991650996669	18.1585474109909	-1.28203791127301	0.212075854894138	0.254491025872966	-6.42090197700496
+11302	"Aatk"	"apoptosis-associated tyrosine kinase"	11	5743	-0.035023968220445	9.78883119065989	-0.236945963165271	0.814709535394086	0.814709535394086	-6.09497670655939
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WTanno.tsv
--- a/test-data/limma-voom_Mut-WTanno.tsv	Thu Sep 07 05:27:27 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-"EntrezID"	"Symbol"	"GeneName"	"Chr"	"Length"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
-11304	"Abca4"	"ATP-binding cassette, sub-family A (ABC1), member 4"	3	7248	0.457332061341026	15.5254133001226	6.50459574633681	9.98720685006039e-07	5.99232411003624e-06	14.0741948485896
-11287	"Pzp"	"pregnancy zone protein"	6	4681	0.190749727701785	17.6546448244617	5.09535410066402	3.26518807654125e-05	9.79556422962375e-05	5.46773893802392
-11298	"Aanat"	"arylalkylamine N-acetyltransferase"	11	1455	-0.138014418336201	17.6747285193431	-3.33168485842331	0.00278753263633162	0.00557506527266324	-1.84301342041449
-11303	"Abca1"	"ATP-binding cassette, sub-family A (ABC1), member 1"	4	10260	-0.0558958943606989	17.886791401216	-1.30108531275576	0.205582481502297	0.254491025872973	-6.4924124057801
-11305	"Abca2"	"ATP-binding cassette, sub-family A (ABC1), member 2"	2	8061	-0.0606991650996633	18.1585474109909	-1.28203791127299	0.212075854894144	0.254491025872973	-6.42090197700503
-11302	"Aatk"	"apoptosis-associated tyrosine kinase"	11	5743	-0.0350239682204432	9.78883119065989	-0.236945963165269	0.814709535394087	0.814709535394087	-6.09497670655944
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_Mut-WTmultifact.tsv
--- a/test-data/limma-voom_Mut-WTmultifact.tsv	Thu Sep 07 05:27:27 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-"GeneID"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
-"11287"	0.190521683937291	17.6546448244617	14.4538942795625	5.93483727834149e-09	3.56090236700489e-08	96.5864491178815
-"11298"	-0.139877413200757	17.6747285193431	-7.71901694642401	5.4113978311412e-06	1.62341934934236e-05	22.3495333961486
-"11304"	0.459011254726061	15.5254133001226	5.64083198409048	0.000108849060923731	0.000217698121847462	8.76657226635859
-"11303"	-0.0641599901594212	17.886791401216	-2.99008181539252	0.0112725655419959	0.0169088483129939	-2.70089035288952
-"11305"	-0.0651479753016204	18.1585474109909	-2.28935282063873	0.0409794711446019	0.0491753653735222	-4.27133526604488
-"11302"	-0.035881713464434	9.78883119065989	-0.439436789626921	0.668154169055758	0.668154169055758	-5.75838315483926
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_WT-Mut.tsv
--- a/test-data/limma-voom_WT-Mut.tsv	Thu Sep 07 05:27:27 2017 -0400
+++ b/test-data/limma-voom_WT-Mut.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -1,7 +1,7 @@
 "GeneID"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
-"11304"	-0.457332061341026	15.5254133001226	-6.50459574633681	9.98720685006039e-07	5.99232411003624e-06	14.0741948485896
-"11287"	-0.190749727701785	17.6546448244617	-5.09535410066402	3.26518807654125e-05	9.79556422962375e-05	5.46773893802392
-"11298"	0.138014418336201	17.6747285193431	3.33168485842331	0.00278753263633162	0.00557506527266324	-1.84301342041449
-"11303"	0.0558958943606989	17.886791401216	1.30108531275576	0.205582481502297	0.254491025872973	-6.4924124057801
-"11305"	0.0606991650996633	18.1585474109909	1.28203791127299	0.212075854894144	0.254491025872973	-6.42090197700503
-"11302"	0.0350239682204432	9.78883119065989	0.236945963165269	0.814709535394087	0.814709535394087	-6.09497670655944
+"11304"	-0.457332061341022	15.5254133001226	-6.50459574633635	9.98720685007144e-07	5.99232411004286e-06	14.0741948485867
+"11287"	-0.190749727701785	17.6546448244617	-5.09535410066427	3.26518807653921e-05	9.79556422961764e-05	5.46773893802514
+"11298"	0.138014418336201	17.6747285193431	3.3316848584234	0.00278753263633103	0.00557506527266206	-1.84301342041422
+"11303"	0.0558958943607024	17.886791401216	1.30108531275582	0.205582481502277	0.254491025872966	-6.49241240578
+"11305"	0.0606991650996669	18.1585474109909	1.28203791127301	0.212075854894138	0.254491025872966	-6.42090197700496
+"11302"	0.035023968220445	9.78883119065989	0.236945963165271	0.814709535394086	0.814709535394086	-6.09497670655939
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_WT-Mut_2fact_anno.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-voom_WT-Mut_2fact_anno.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"EntrezID"	"Symbol"	"GeneName"	"Chr"	"Length"	"logFC"	"AveExpr"	"t"	"P.Value"	"adj.P.Val"	"B"
+11287	"Pzp"	"pregnancy zone protein"	6	4681	-0.190521683937302	17.6546448244617	-15.1152647604854	3.5605743515492e-09	2.13634461092952e-08	106.319792304784
+11298	"Aanat"	"arylalkylamine N-acetyltransferase"	11	1455	0.139877413200757	17.6747285193431	8.49164953574146	2.03124167559177e-06	6.0937250267753e-06	28.5167891813121
+11304	"Abca4"	"ATP-binding cassette, sub-family A (ABC1), member 4"	3	7248	-0.459011254726057	15.5254133001226	-5.60316118539212	0.000115556485215252	0.000231112970430504	8.56249981126828
+11303	"Abca1"	"ATP-binding cassette, sub-family A (ABC1), member 1"	4	10260	0.0641599901594141	17.886791401216	2.92798536696768	0.0126512435710263	0.0189768653565394	-2.86386681629442
+11305	"Abca2"	"ATP-binding cassette, sub-family A (ABC1), member 2"	2	8061	0.0651479753016169	18.1585474109909	2.03911882030581	0.0640903140266962	0.0769083768320354	-4.69855556350315
+11302	"Aatk"	"apoptosis-associated tyrosine kinase"	11	5743	0.0358817134644287	9.78883119065989	0.43940476739557	0.668176726641997	0.668176726641997	-5.75841999159995
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_normcounts.tsv
--- a/test-data/limma-voom_normcounts.tsv	Thu Sep 07 05:27:27 2017 -0400
+++ b/test-data/limma-voom_normcounts.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -1,7 +1,7 @@
-"GeneID"	"WT1"	"WT2"	"WT3"	"Mut1"	"Mut2"	"Mut3"
-"11287"	17.6076545522668	17.5079905058358	17.5639197719462	17.7719335903598	17.7105244453357	17.7658460810258
-"11298"	17.7727137985373	17.6986875511432	17.7598828784201	17.6506532355915	17.5520012519588	17.6144324004081
-"11302"	9.57719962386619	10.0175525101992	9.82560228478005	9.71615817858834	9.91760904198188	9.67886550454371
-"11303"	17.9125899785601	17.8773613214092	17.955228759658	17.8774046020864	17.7865502833631	17.911613462219
-"11304"	15.354518172169	15.2178020484983	15.3174553811097	15.7545783969022	15.8519803740125	15.6561454280436
-"11305"	18.1808259688524	18.1818679259847	18.2026681768366	18.0401341021246	18.1408682071359	18.204920085011
+"GeneID"	"Mut1"	"Mut2"	"Mut3"	"WT1"	"WT2"	"WT3"
+"11287"	17.7719335903598	17.7105244453357	17.7658460810258	17.6076545522668	17.5079905058358	17.5639197719462
+"11298"	17.6506532355915	17.5520012519588	17.6144324004081	17.7727137985373	17.6986875511432	17.7598828784201
+"11302"	9.71615817858834	9.91760904198188	9.67886550454371	9.57719962386619	10.0175525101992	9.82560228478005
+"11303"	17.8774046020864	17.7865502833631	17.911613462219	17.9125899785601	17.8773613214092	17.955228759658
+"11304"	15.7545783969022	15.8519803740125	15.6561454280436	15.354518172169	15.2178020484983	15.3174553811097
+"11305"	18.0401341021246	18.1408682071359	18.204920085011	18.1808259688524	18.1818679259847	18.2026681768366
diff -r a330ddf43861 -r 38aab66ae5cb test-data/limma-voom_normcounts_anno.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/limma-voom_normcounts_anno.tsv	Wed Jan 31 12:45:42 2018 -0500
@@ -0,0 +1,7 @@
+"EntrezID"	"Symbol"	"GeneName"	"Chr"	"Length"	"Mut1"	"Mut2"	"Mut3"	"WT1"	"WT2"	"WT3"
+11287	"Pzp"	"pregnancy zone protein"	6	4681	17.7719335903598	17.7105244453357	17.7658460810258	17.6076545522668	17.5079905058358	17.5639197719462
+11298	"Aanat"	"arylalkylamine N-acetyltransferase"	11	1455	17.6506532355915	17.5520012519588	17.6144324004081	17.7727137985373	17.6986875511432	17.7598828784201
+11302	"Aatk"	"apoptosis-associated tyrosine kinase"	11	5743	9.71615817858834	9.91760904198188	9.67886550454371	9.57719962386619	10.0175525101992	9.82560228478005
+11303	"Abca1"	"ATP-binding cassette, sub-family A (ABC1), member 1"	4	10260	17.8774046020864	17.7865502833631	17.911613462219	17.9125899785601	17.8773613214092	17.955228759658
+11304	"Abca4"	"ATP-binding cassette, sub-family A (ABC1), member 4"	3	7248	15.7545783969022	15.8519803740125	15.6561454280436	15.354518172169	15.2178020484983	15.3174553811097
+11305	"Abca2"	"ATP-binding cassette, sub-family A (ABC1), member 2"	2	8061	18.0401341021246	18.1408682071359	18.204920085011	18.1808259688524	18.1818679259847	18.2026681768366
diff -r a330ddf43861 -r 38aab66ae5cb test-data/matrix.txt
--- a/test-data/matrix.txt	Thu Sep 07 05:27:27 2017 -0400
+++ b/test-data/matrix.txt	Wed Jan 31 12:45:42 2018 -0500
@@ -1,7 +1,7 @@
-GeneID	WT1	WT2	WT3	Mut1	Mut2	Mut3
-11287	1699	1528	1601	1463	1441	1495
-11298	1905	1744	1834	1345	1291	1346
-11302	6	8	7	5	6	5
-11303	2099	1974	2100	1574	1519	1654
-11304	356	312	337	361	397	346
-11305	2528	2438	2493	1762	1942	2027
+GeneID	Mut1	Mut2	Mut3	WT1	WT2	WT3
+11287	1463	1441	1495	1699	1528	1601
+11298	1345	1291	1346	1905	1744	1834
+11302	5	6	5	6	8	7
+11303	1574	1519	1654	2099	1974	2100
+11304	361	397	346	356	312	337
+11305	1762	1942	2027	2528	2438	2493