Mercurial > repos > iuc > edger
changeset 14:c5fa04118f83 draft default tip
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/edger commit 4b17306763bc8eb92c8967c7b4b57abcc514e670
| author | iuc |
|---|---|
| date | Wed, 04 Sep 2024 15:49:53 +0000 |
| parents | 838b481dc6f9 |
| children | |
| files | edger.R edger.xml test-data/contrasts_file.txt |
| diffstat | 3 files changed, 474 insertions(+), 357 deletions(-) [+] |
line wrap: on
line diff
--- a/edger.R Wed Sep 27 19:52:53 2023 +0000 +++ b/edger.R Wed Sep 04 15:49:53 2024 +0000 @@ -8,6 +8,7 @@ # 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 +# formula", "F", 2, "character". -String containing a formula to override default use of factInput # 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 @@ -41,8 +42,8 @@ # setup R error handling to go to stderr options(show.error.messages = FALSE, error = function() { - cat(geterrmessage(), file = stderr()) - q("no", 1, FALSE) + cat(geterrmessage(), file = stderr()) + q("no", 1, FALSE) }) # we need that to not crash galaxy with an UTF8 error on German LC settings. @@ -63,35 +64,41 @@ # Function to sanitise contrast equations so there are no whitespaces # surrounding the arithmetic operators, leading or trailing whitespace sanitise_equation <- 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 sanitise_groups <- 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) +} + +# Sanitise file base names coming from factors or contrasts +sanitise_basename <- function(string) { + string <- gsub("[/^]", "_", string) + return(string) } # Generate output folder and paths make_out <- function(filename) { - return(paste0(out_path, "/", filename)) + return(paste0(out_path, "/", filename)) } # Generating design information paste_listname <- function(string) { - return(paste0("factors$", string)) + return(paste0("factors$", string)) } # Create cata function: default path set, default seperator empty and appending @@ -99,49 +106,49 @@ # defaults) 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)) + 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)) + .Internal(cat(list(...), file, sep, fill, labels, append)) } # Function to write code for html head and title html_head <- 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 html_link <- 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 html_image <- 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 list_item <- function(...) { - cata("<li>", ..., "</li>\n") + cata("<li>", ..., "</li>\n") } table_item <- function(...) { - cata("<td>", ..., "</td>\n") + cata("<td>", ..., "</td>\n") } table_head_item <- function(...) { - cata("<th>", ..., "</th>\n") + cata("<th>", ..., "</th>\n") } ################################################################################ @@ -153,167 +160,174 @@ # 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", - "lrtOpt", "t", 0, "logical" -), -byrow = TRUE, ncol = 4 +spec <- matrix( + c( + "htmlPath", "R", 1, "character", + "outPath", "o", 1, "character", + "filesPath", "j", 2, "character", + "matrixPath", "m", 2, "character", + "factFile", "f", 2, "character", + "formula", "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", + "lrtOpt", "t", 0, "logical" + ), + byrow = TRUE, ncol = 4 ) opt <- getopt(spec) 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) + cat("A counts matrix (or a set of counts files) is required.\n") + q(status = 1) } if (is.null(opt$cpmReq)) { - filt_cpm <- FALSE + filt_cpm <- FALSE } else { - filt_cpm <- TRUE + filt_cpm <- TRUE } if (is.null(opt$cntReq) || is.null(opt$sampleReq)) { - filt_smpcount <- FALSE + filt_smpcount <- FALSE } else { - filt_smpcount <- TRUE + filt_smpcount <- TRUE } if (is.null(opt$totReq)) { - filt_totcount <- FALSE + filt_totcount <- FALSE } else { - filt_totcount <- TRUE + filt_totcount <- TRUE } if (is.null(opt$lrtOpt)) { - want_lrt <- FALSE + want_lrt <- FALSE } else { - want_lrt <- TRUE + want_lrt <- TRUE } if (is.null(opt$rdaOpt)) { - want_rda <- FALSE + want_rda <- FALSE } else { - want_rda <- TRUE + want_rda <- TRUE } if (is.null(opt$annoPath)) { - have_anno <- FALSE + have_anno <- FALSE } else { - have_anno <- TRUE + have_anno <- TRUE } if (is.null(opt$normCounts)) { - want_norm <- FALSE + want_norm <- FALSE } else { - want_norm <- TRUE + want_norm <- TRUE } if (is.null(opt$robOpt)) { - want_robust <- FALSE + want_robust <- FALSE } else { - want_robust <- TRUE + want_robust <- TRUE } if (!is.null(opt$filesPath)) { - # Process the separate count files (adapted from DESeq2 wrapper) - library("rjson") - parser <- newJSONParser() - parser$addData(opt$filesPath) - factor_list <- parser$getObject() - factors <- sapply(factor_list, function(x) x[[1]]) - filenames_in <- unname(unlist(factor_list[[1]][[2]])) - sampletable <- data.frame( - sample = basename(filenames_in), - filename = filenames_in, - row.names = filenames_in, - stringsAsFactors = FALSE - ) - for (factor in factor_list) { - 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] + # Process the separate count files (adapted from DESeq2 wrapper) + library("rjson") + parser <- newJSONParser() + parser$addData(opt$filesPath) + factor_list <- parser$getObject() + factors <- sapply(factor_list, function(x) x[[1]]) + filenames_in <- unname(unlist(factor_list[[1]][[2]])) + sampletable <- data.frame( + sample = basename(filenames_in), + filename = filenames_in, + row.names = filenames_in, + stringsAsFactors = FALSE + ) + for (factor in factor_list) { + 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) } - sampletable[[factorname]] <- factor(sampletable[[factorname]], levels = lvls) - } - rownames(sampletable) <- sampletable$sample - rem <- c("sample", "filename") - factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE] + 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) + # 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", strip.white = TRUE, stringsAsFactors = FALSE) - row.names(counts) <- counts[, 1] - counts <- counts[, -1] - countsrows <- nrow(counts) + # Process the single count matrix + counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, 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", strip.white = TRUE, stringsAsFactors = TRUE) - # check samples names match - if (!any(factordata[, 1] %in% colnames(counts))) { - stop("Sample IDs in factors file and count matrix don't match") + # Process factors + if (is.null(opt$factInput)) { + factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE) + # check samples names match + if (!any(factordata[, 1] %in% colnames(counts))) { + stop("Sample IDs in factors file and count matrix don't match") + } + # order samples as in counts matrix + factordata <- factordata[match(colnames(counts), factordata[, 1]), ] + factors <- data.frame(sapply(factordata[, -1, drop = FALSE], make.names)) + } 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, sanitise_groups) + factordata <- sapply(factordata, strsplit, split = ",") + factordata <- sapply(factordata, make.names) + # Transform factor data into data frame of R factor objects + factors <- data.frame(factordata) } - # order samples as in counts matrix - factordata <- factordata[match(colnames(counts), factordata[, 1]), ] - factors <- data.frame(sapply(factordata[, -1, drop = FALSE], make.names)) - } 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, sanitise_groups) - 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 (have_anno) { - geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE) + geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE) } # Create output directory out_path <- opt$outPath dir.create(out_path, showWarnings = FALSE) -# Split up contrasts separated by comma into a vector then sanitise -contrast_data <- unlist(strsplit(opt$contrastData, split = ",")) +# Check if contrastData is a file or not +if (file.exists(opt$contrastData)) { + contrast_data <- unlist(read.table(opt$contrastData, sep = "\t", header = TRUE)[[1]]) +} else { + # Split up contrasts separated by comma into a vector then sanitise + contrast_data <- unlist(strsplit(opt$contrastData, split = ",")) +} contrast_data <- sanitise_equation(contrast_data) contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE) @@ -324,16 +338,16 @@ mds_pdf <- character() # Initialise character vector mds_png <- character() for (i in seq_len(ncol(factors))) { - mds_pdf[i] <- make_out(paste0("mdsplot_", names(factors)[i], ".pdf")) - mds_png[i] <- make_out(paste0("mdsplot_", names(factors)[i], ".png")) + mds_pdf[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf")) + mds_png[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png")) } md_pdf <- character() md_png <- character() top_out <- character() for (i in seq_along(contrast_data)) { - md_pdf[i] <- make_out(paste0("mdplot_", contrast_data[i], ".pdf")) - md_png[i] <- make_out(paste0("mdplot_", contrast_data[i], ".png")) - top_out[i] <- make_out(paste0("edgeR_", contrast_data[i], ".tsv")) + md_pdf[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf")) + md_png[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png")) + top_out[i] <- make_out(paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")) } # Save output paths for each contrast as vectors norm_out <- make_out("edgeR_normcounts.tsv") rda_out <- make_out("edgeR_analysis.RData") @@ -357,11 +371,11 @@ data <- list() data$counts <- counts if (have_anno) { - # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) - annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ] - data$genes <- annoord + # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) + annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ] + data$genes <- annoord } else { - data$genes <- data.frame(GeneID = row.names(counts)) + data$genes <- data.frame(GeneID = row.names(counts)) } # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of @@ -369,16 +383,16 @@ prefilter_count <- nrow(data$counts) if (filt_cpm || filt_smpcount || filt_totcount) { - if (filt_totcount) { - keep <- rowSums(data$counts) >= opt$cntReq - } else if (filt_smpcount) { - keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq - } else if (filt_cpm) { - keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq - } + if (filt_totcount) { + keep <- rowSums(data$counts) >= opt$cntReq + } else if (filt_smpcount) { + keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq + } else if (filt_cpm) { + keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq + } - data$counts <- data$counts[keep, ] - data$genes <- data$genes[keep, , drop = FALSE] + data$counts <- data$counts[keep, ] + data$genes <- data$genes[keep, , drop = FALSE] } postfilter_count <- nrow(data$counts) @@ -397,26 +411,33 @@ data$genes <- genes - -formula <- "~0" -for (i in seq_along(factor_list)) { - formula <- paste(formula, factor_list[i], sep = "+") +if (!is.null(opt$formula)) { + formula <- opt$formula + # sanitisation can be getting rid of the "~" + if (!startsWith(formula, "~")) { + formula <- paste0("~", formula) + } +} else { + formula <- "~0" + for (i in seq_along(factor_list)) { + formula <- paste(formula, factor_list[i], sep = "+") + } } formula <- formula(formula) design <- model.matrix(formula, factors) for (i in seq_along(factor_list)) { - colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE) + colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE) } # Calculating normalising factor, estimating dispersion data <- calcNormFactors(data, method = opt$normOpt) if (want_robust) { - data <- estimateDisp(data, design = design, robust = TRUE) + data <- estimateDisp(data, design = design, robust = TRUE) } else { - data <- estimateDisp(data, design = design) + data <- estimateDisp(data, design = design) } # Generate contrasts information @@ -432,35 +453,35 @@ # MDS plot png(mds_png, width = 600, height = 600) plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1])) -img_name <- paste0("MDS Plot_", names(factors)[1], ".png") -img_addr <- paste0("mdsplot_", names(factors)[1], ".png") +img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[1]), ".png") +img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[1]), ".png") image_data[1, ] <- c(img_name, img_addr) invisible(dev.off()) pdf(mds_pdf) plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1])) -link_name <- paste0("MDS Plot_", names(factors)[1], ".pdf") -link_addr <- paste0("mdsplot_", names(factors)[1], ".pdf") +link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[1]), ".pdf") +link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[1]), ".pdf") link_data[1, ] <- c(link_name, link_addr) invisible(dev.off()) # If additional factors create additional MDS plots coloured by factor if (ncol(factors) > 1) { - for (i in 2:ncol(factors)) { - png(mds_png[i], width = 600, height = 600) - plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) - img_name <- paste0("MDS Plot_", names(factors)[i], ".png") - img_addr <- paste0("mdsplot_", names(factors)[i], ".png") - image_data <- rbind(image_data, c(img_name, img_addr)) - invisible(dev.off()) + for (i in 2:ncol(factors)) { + png(mds_png[i], width = 600, height = 600) + plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) + img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".png") + img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png") + image_data <- rbind(image_data, c(img_name, img_addr)) + invisible(dev.off()) - pdf(mds_pdf[i]) - plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) - link_name <- paste0("MDS Plot_", names(factors)[i], ".pdf") - link_addr <- paste0("mdsplot_", names(factors)[i], ".pdf") - link_data <- rbind(link_data, c(link_name, link_addr)) - invisible(dev.off()) - } + pdf(mds_pdf[i]) + plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) + link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".pdf") + link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf") + link_data <- rbind(link_data, c(link_name, link_addr)) + invisible(dev.off()) + } } # BCV Plot @@ -480,111 +501,111 @@ # Generate fit if (want_lrt) { - fit <- glmFit(data, design) + fit <- glmFit(data, design) } else { - if (want_robust) { - fit <- glmQLFit(data, design, robust = TRUE) - } else { - fit <- glmQLFit(data, design) - } + if (want_robust) { + fit <- glmQLFit(data, design, robust = TRUE) + } else { + fit <- glmQLFit(data, design) + } - # Plot QL dispersions - png(ql_png, width = 600, height = 600) - plotQLDisp(fit, main = "QL Plot") - img_name <- "QL Plot" - img_addr <- "qlplot.png" - image_data <- rbind(image_data, c(img_name, img_addr)) - invisible(dev.off()) + # Plot QL dispersions + png(ql_png, width = 600, height = 600) + plotQLDisp(fit, main = "QL Plot") + img_name <- "QL Plot" + img_addr <- "qlplot.png" + image_data <- rbind(image_data, c(img_name, img_addr)) + invisible(dev.off()) - pdf(ql_pdf) - plotQLDisp(fit, main = "QL Plot") - link_name <- "QL Plot.pdf" - link_addr <- "qlplot.pdf" - link_data <- rbind(link_data, c(link_name, link_addr)) - invisible(dev.off()) + pdf(ql_pdf) + plotQLDisp(fit, main = "QL Plot") + link_name <- "QL Plot.pdf" + link_addr <- "qlplot.pdf" + link_data <- rbind(link_data, c(link_name, link_addr)) + invisible(dev.off()) } # Save normalised counts (log2cpm) if (want_norm) { - normalised_counts <- cpm(data, normalized.lib.sizes = TRUE, log = TRUE) - normalised_counts <- data.frame(data$genes, normalised_counts) - write.table(normalised_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE) - link_data <- rbind(link_data, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv")) + normalised_counts <- cpm(data, normalized.lib.sizes = TRUE, log = TRUE) + normalised_counts <- data.frame(data$genes, normalised_counts) + write.table(normalised_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE) + link_data <- rbind(link_data, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv")) } for (i in seq_along(contrast_data)) { - if (want_lrt) { - res <- glmLRT(fit, contrast = contrasts[, i]) - } else { - res <- glmQLFTest(fit, contrast = contrasts[, i]) - } + if (want_lrt) { + res <- glmLRT(fit, contrast = contrasts[, i]) + } else { + res <- glmQLFTest(fit, contrast = contrasts[, i]) + } - status <- decideTestsDGE(res, - adjust.method = opt$pAdjOpt, p.value = opt$pValReq, - lfc = opt$lfcReq - ) - sum_status <- summary(status) + status <- decideTestsDGE(res, + adjust.method = opt$pAdjOpt, p.value = opt$pValReq, + lfc = opt$lfcReq + ) + sum_status <- summary(status) - # Collect counts for differential expression - up_count[i] <- sum_status["Up", ] - down_count[i] <- sum_status["Down", ] - flat_count[i] <- sum_status["NotSig", ] + # Collect counts for differential expression + up_count[i] <- sum_status["Up", ] + down_count[i] <- sum_status["Down", ] + flat_count[i] <- sum_status["NotSig", ] - # Write top expressions table - top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue") - write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE) + # Write top expressions table + top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue") + write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE) - link_name <- paste0("edgeR_", contrast_data[i], ".tsv") - link_addr <- paste0("edgeR_", contrast_data[i], ".tsv") - link_data <- rbind(link_data, c(link_name, link_addr)) + link_name <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv") + link_addr <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv") + link_data <- rbind(link_data, c(link_name, link_addr)) - # Plot MD (log ratios vs mean difference) using limma package - pdf(md_pdf[i]) - limma::plotMD(res, - status = status, - main = paste("MD Plot:", unmake_names(contrast_data[i])), - hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), - xlab = "Average Expression", ylab = "logFC" - ) + # Plot MD (log ratios vs mean difference) using limma package + pdf(md_pdf[i]) + limma::plotMD(res, + status = status, + main = paste("MD Plot:", unmake_names(contrast_data[i])), + hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), + xlab = "Average Expression", ylab = "logFC" + ) - abline(h = 0, col = "grey", lty = 2) + abline(h = 0, col = "grey", lty = 2) - link_name <- paste0("MD Plot_", contrast_data[i], ".pdf") - link_addr <- paste0("mdplot_", contrast_data[i], ".pdf") - link_data <- rbind(link_data, c(link_name, link_addr)) - invisible(dev.off()) + link_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".pdf") + link_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf") + link_data <- rbind(link_data, c(link_name, link_addr)) + invisible(dev.off()) - png(md_png[i], height = 600, width = 600) - limma::plotMD(res, - status = status, - main = paste("MD Plot:", unmake_names(contrast_data[i])), - hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), - xlab = "Average Expression", ylab = "logFC" - ) + png(md_png[i], height = 600, width = 600) + limma::plotMD(res, + status = status, + main = paste("MD Plot:", unmake_names(contrast_data[i])), + hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), + xlab = "Average Expression", ylab = "logFC" + ) - abline(h = 0, col = "grey", lty = 2) + abline(h = 0, col = "grey", lty = 2) - img_name <- paste0("MD Plot_", contrast_data[i], ".png") - img_addr <- paste0("mdplot_", contrast_data[i], ".png") - image_data <- rbind(image_data, c(img_name, img_addr)) - invisible(dev.off()) + img_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".png") + img_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png") + image_data <- rbind(image_data, c(img_name, img_addr)) + invisible(dev.off()) } sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count) row.names(sig_diff) <- contrast_data # Save relevant items as rda object if (want_rda) { - if (want_norm) { - save(counts, data, status, normalised_counts, labels, factors, fit, res, top, contrasts, design, - file = rda_out, ascii = TRUE - ) - } else { - save(counts, data, status, labels, factors, fit, res, top, contrasts, design, - file = rda_out, ascii = TRUE - ) - } - link_data <- rbind(link_data, c("edgeR_analysis.RData", "edgeR_analysis.RData")) + if (want_norm) { + save(counts, data, status, normalised_counts, labels, factors, fit, res, top, contrasts, design, + file = rda_out, ascii = TRUE + ) + } else { + save(counts, data, status, labels, factors, fit, res, top, contrasts, design, + file = rda_out, ascii = TRUE + ) + } + link_data <- rbind(link_data, c("edgeR_analysis.RData", "edgeR_analysis.RData")) } # Record session info @@ -612,7 +633,7 @@ html_image(image_data$Link[1], image_data$Label[1]) for (i in 2:nrow(image_data)) { - html_image(image_data$Link[i], image_data$Label[i]) + html_image(image_data$Link[i], image_data$Label[i]) } cata("<h4>Differential Expression Counts:</h4>\n") @@ -621,40 +642,40 @@ cata("<tr>\n") table_item() for (i in colnames(sig_diff)) { - table_head_item(i) + table_head_item(i) } cata("</tr>\n") for (i in seq_len(nrow(sig_diff))) { - cata("<tr>\n") - table_head_item(unmake_names(row.names(sig_diff)[i])) - for (j in seq_len(ncol(sig_diff))) { - table_item(as.character(sig_diff[i, j])) - } - cata("</tr>\n") + cata("<tr>\n") + table_head_item(unmake_names(row.names(sig_diff)[i])) + for (j in seq_len(ncol(sig_diff))) { + table_item(as.character(sig_diff[i, j])) + } + cata("</tr>\n") } cata("</table>") cata("<h4>Plots:</h4>\n") for (i in seq_len(nrow(link_data))) { - if (grepl(".pdf", link_data$Link[i])) { - html_link(link_data$Link[i], link_data$Label[i]) - } + if (grepl(".pdf", link_data$Link[i])) { + html_link(link_data$Link[i], link_data$Label[i]) + } } cata("<h4>Tables:</h4>\n") for (i in seq_len(nrow(link_data))) { - if (grepl(".tsv", link_data$Link[i])) { - html_link(link_data$Link[i], link_data$Label[i]) - } + if (grepl(".tsv", link_data$Link[i])) { + html_link(link_data$Link[i], link_data$Label[i]) + } } if (want_rda) { - cata("<h4>R Data Objects:</h4>\n") - for (i in seq_len(nrow(link_data))) { - if (grepl(".RData", link_data$Link[i])) { - html_link(link_data$Link[i], link_data$Label[i]) + cata("<h4>R Data Objects:</h4>\n") + for (i in seq_len(nrow(link_data))) { + if (grepl(".RData", link_data$Link[i])) { + html_link(link_data$Link[i], link_data$Label[i]) + } } - } } cata("<p>Alt-click links to download file.</p>\n") @@ -666,68 +687,68 @@ cata("<ul>\n") if (filt_cpm || filt_smpcount || filt_totcount) { - if (filt_cpm) { - temp_str <- paste( - "Genes without more than", opt$cpmReq, - "CPM in at least", opt$sampleReq, "samples are insignificant", - "and filtered out." - ) - } else if (filt_smpcount) { - temp_str <- paste( - "Genes without more than", opt$cntReq, - "counts in at least", opt$sampleReq, "samples are insignificant", - "and filtered out." + if (filt_cpm) { + temp_str <- paste( + "Genes without more than", opt$cpmReq, + "CPM in at least", opt$sampleReq, "samples are insignificant", + "and filtered out." + ) + } else if (filt_smpcount) { + temp_str <- paste( + "Genes without more than", opt$cntReq, + "counts in at least", opt$sampleReq, "samples are insignificant", + "and filtered out." + ) + } else if (filt_totcount) { + temp_str <- paste( + "Genes without more than", opt$cntReq, + "counts, after summing counts for all samples, are insignificant", + "and filtered out." + ) + } + + list_item(temp_str) + filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2) + temp_str <- paste0( + filtered_count, " of ", prefilter_count, " (", filter_prop, + "%) genes were filtered out for low expression." ) - } else if (filt_totcount) { - temp_str <- paste( - "Genes without more than", opt$cntReq, - "counts, after summing counts for all samples, are insignificant", - "and filtered out." - ) - } - - list_item(temp_str) - filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2) - temp_str <- paste0( - filtered_count, " of ", prefilter_count, " (", filter_prop, - "%) genes were filtered out for low expression." - ) - list_item(temp_str) + list_item(temp_str) } list_item(opt$normOpt, " was the method used to normalise library sizes.") if (want_lrt) { - list_item("The edgeR likelihood ratio test was used.") + list_item("The edgeR likelihood ratio test was used.") } else { - if (want_robust) { - list_item("The edgeR quasi-likelihood test was used with robust settings (robust=TRUE with estimateDisp and glmQLFit).") - } else { - list_item("The edgeR quasi-likelihood test was used.") - } + if (want_robust) { + list_item("The edgeR quasi-likelihood test was used with robust settings (robust=TRUE with estimateDisp and glmQLFit).") + } else { + list_item("The edgeR quasi-likelihood test was used.") + } } if (opt$pAdjOpt != "none") { - if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") { + if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") { + temp_str <- paste0( + "MD-Plot highlighted genes are significant at FDR ", + "of ", opt$pValReq, " and exhibit log2-fold-change of at ", + "least ", opt$lfcReq, "." + ) + list_item(temp_str) + } else if (opt$pAdjOpt == "holm") { + temp_str <- paste0( + "MD-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, "." + ) + list_item(temp_str) + } +} else { temp_str <- paste0( - "MD-Plot highlighted genes are significant at FDR ", - "of ", opt$pValReq, " and exhibit log2-fold-change of at ", - "least ", opt$lfcReq, "." + "MD-Plot highlighted genes are significant at p-value ", + "of ", opt$pValReq, " and exhibit log2-fold-change of at ", + "least ", opt$lfcReq, "." ) list_item(temp_str) - } else if (opt$pAdjOpt == "holm") { - temp_str <- paste0( - "MD-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, "." - ) - list_item(temp_str) - } -} else { - temp_str <- paste0( - "MD-Plot highlighted genes are significant at p-value ", - "of ", opt$pValReq, " and exhibit log2-fold-change of at ", - "least ", opt$lfcReq, "." - ) - list_item(temp_str) } cata("</ul>\n") @@ -741,26 +762,26 @@ table_head_item(names(factors)[1], " (Primary Factor)") if (ncol(factors) > 1) { - for (i in names(factors)[2:length(names(factors))]) { - table_head_item(i) - } - cata("</tr>\n") + for (i in names(factors)[2:length(names(factors))]) { + table_head_item(i) + } + cata("</tr>\n") } for (i in seq_len(nrow((factors)))) { - cata("<tr>\n") - table_head_item(row.names(factors)[i]) - for (j in seq_len(ncol(factors))) { - table_item(as.character(unmake_names(factors[i, j]))) - } - cata("</tr>\n") + cata("<tr>\n") + table_head_item(row.names(factors)[i]) + for (j in seq_len(ncol(factors))) { + table_item(as.character(unmake_names(factors[i, j]))) + } + cata("</tr>\n") } cata("</table>") for (i in seq_len(nrow(link_data))) { - if (grepl("session_info", link_data$Link[i])) { - html_link(link_data$Link[i], link_data$Label[i]) - } + if (grepl("session_info", link_data$Link[i])) { + html_link(link_data$Link[i], link_data$Label[i]) + } } cata("<table border=\"0\">\n")
--- a/edger.xml Wed Sep 27 19:52:53 2023 +0000 +++ b/edger.xml Wed Sep 04 15:49:53 2024 +0000 @@ -4,7 +4,7 @@ </description> <macros> <token name="@TOOL_VERSION@">3.36.0</token> - <token name="@VERSION_SUFFIX@">2</token> + <token name="@VERSION_SUFFIX@">5</token> </macros> <edam_topics> <edam_topic>topic_3308</edam_topic> @@ -68,7 +68,15 @@ -a '$anno.geneanno' #end if --C '${ ','.join( ['%s' % $x.contrast for x in $rep_contrast] ) }' +#if $formula: + -F '$formula' +#end if + +#if $contrasts.contrastOpt == 'file': + -C '$contrasts.cinfo' +#else: + -C '${ ','.join( ['%s' % $x.contrast for x in $contrasts.rep_contrast] ) }' +#end if #if $f.filt.filt_select == 'yes': #if $f.filt.cformat.format_select == 'cpm': @@ -122,11 +130,8 @@ <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, and the first character of each factor must be a letter"> - <sanitizer> - <valid initial="string.letters,string.digits"> - <add value="_"/> - </valid> - </sanitizer> + <validator type="empty_field"/> + <validator type="regex" message="Please only use letters, numbers or underscores">^[\w]+$</validator> </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 belong to (e.g. WT or Mut). NOTE: Please only use letters, numbers or underscores (case sensitive), and the first character of each group must be a letter"> @@ -176,13 +181,45 @@ </when> <when value="no"/> </conditional> + <!-- Optional formula --> + <param name="formula" type="text" optional="true" label="Formula for linear model" help="An optional formula for the EdgeR linear model, this will override the use of the fields in factors as a simple sum. The formula can only use elements available in the factors file. This needs to be exactly as EdgeR expect the formula, ie. `~ 0 + factor_A + factor_B:factor_C`. See EdgeR documentation for more details."> + <sanitizer invalid_char=""> + <valid initial="string.letters,string.digits"> + <add value="_"/> + <add value="-"/> + <add value="+"/> + <add value="*"/> + <add value="/"/> + <add value="^"/> + <add value=":"/> + <add value="."/> + <add value="~"/> + <add value=" "/> + <add value="("/> + <add value=")"/> + <add value="@"/> + <add value="$"/> + </valid> + </sanitizer> + </param> <!-- 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 differences between contrasts use e.g. (MT.t1-MT.t0)-(WT.t1-WT.t0). 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 or https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf page 36 for nested comparisons."> - <validator type="empty_field"/> - <validator type="regex" message="Please only use letters, numbers, parentheses or underscores">^[\w\-()]+$</validator> + <conditional name="contrasts"> + <param name="contrastOpt" type="select" label="Input contrasts manually or through a file"> + <option value="manual">manually</option> + <option value="file">file</option> </param> - </repeat> + <when value="manual"> + <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 differences between contrasts use e.g. (MT.t1-MT.t0)-(WT.t1-WT.t0). 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 or https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf page 36 for nested comparisons."> + <validator type="empty_field"/> + <validator type="regex" message="Please only use letters, numbers, parentheses or underscores">^[\w\-()]+$</validator> + </param> + </repeat> + </when> + <when value="file"> + <param name="cinfo" optional="true" type="data" format="tabular" label="Contrasts File" help="Setting this file will ignore any manually added contrasts above, make sure to remove any contrast fields above pressing the trash bin icon, or the tool will fail. First line of the file must be a header, below that each separate contrast should be on a line. Contrast formulas need to be based on ther factors data and potentially the formula provided. See EdgeR documentation on contrasts for more details."/> + </when> + </conditional> <!-- Filter Options --> <section name="f" expanded="false" title="Filter Low Counts"> <conditional name="filt"> @@ -264,6 +301,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -302,6 +340,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="MutA,MutA,MutA,MutB,MutB,MutB,WTA,WTA,WTA,WTB,WTB,WTB"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="(MutA-MutB)-(WTA-WTB)"/> </repeat> @@ -333,6 +372,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -356,6 +396,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -383,6 +424,7 @@ <param name="factorName" value="Batch"/> <param name="groupNames" value="b1,b2,b3,b1,b2,b3"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -402,6 +444,7 @@ <param name="ffile" value="yes"/> <param name="finfo" value="factorinfo.txt"/> <param name="counts" value="matrix.txt"/> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -424,6 +467,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -451,6 +495,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -494,6 +539,7 @@ </repeat> <param name="annoOpt" value="yes"/> <param name="geneanno" value="anno.txt"/> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -530,6 +576,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -565,6 +612,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -599,6 +647,7 @@ <param name="factorName" value="Genotype"/> <param name="groupNames" value="Mut,Mut,Mut,WT,WT,WT"/> </repeat> + <param name="contrastOpt" value="manual"/> <repeat name="rep_contrast"> <param name="contrast" value="Mut-WT"/> </repeat> @@ -626,6 +675,41 @@ </element> </output_collection> </test> + <!-- Ensure formula and contrast file work --> + <test expect_num_outputs="2"> + <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_factor"> + <param name="factorName" value="Batch"/> + <param name="groupNames" value="b1,b2,b3,b1,b2,b3"/> + </repeat> + <param name="contrastOpt" value="file"/> + <param name="cinfo" value="contrasts_file.txt"/> + <param name="formula" value="~ 0 + Genotype + Batch"/> + <param name="normalisationOption" value="TMM"/> + <output_collection name="outTables" count="3"> + <element name="edgeR_Mut-WT" ftype="tabular"> + <assert_contents> + <has_text_matching expression="GeneID.*logFC.*logCPM.*F.*PValue.*FDR"/> + <has_text_matching expression="11304.*0.4584"/> + </assert_contents> + </element> + <element name="edgeR_WT-Mut" ftype="tabular"> + <assert_contents> + <has_text_matching expression="GeneID.*logFC.*logCPM.*F.*PValue.*FDR"/> + </assert_contents> + </element> + <element name="edgeR_(2*Mut_3*WT)-WT" ftype="tabular"> + <assert_contents> + <has_text_matching expression="GeneID.*logFC.*logCPM.*F.*PValue.*FDR"/> + </assert_contents> + </element> + </output_collection> + </test> </tests> <help><![CDATA[ .. class:: infomark @@ -713,12 +797,20 @@ *Groups:* The names of the groups for the factor. The names should start with a letter, and only contain letters, numbers and underscores, other characters such as spaces and hyphens must not be used. If entered into the tool form above, the order must be the same as the samples (to which the groups correspond) are listed in the columns of the counts matrix, with the values separated by commas. +**Formula:** +By default the tool will construct a formula for modelling counts based on the contents of the factors files or the factors given. +This can be overriden by directly providing the EdgeR formula in section named Formula. + **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 must be entered separately using the Insert Contrast button, spaces must not be used. +Alternatively, you can specify a file with contrasts. The file must contain a header (it's value is irrelevant) +and one contrast per line on the first column (other columns are ignored). If using this option, make sure to +remove any contrast section from the manual part, or the tool will fail. + **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