diff deseq2.R @ 30:8fe98f7094de draft

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/deseq2 commit 6868b66f73ddbe947986d1a20b546873cbd515a9

--- a/deseq2.R	Mon Nov 29 18:16:48 2021 +0000
+++ b/deseq2.R	Fri Aug 26 11:16:15 2022 +0000
@@ -31,8 +31,9 @@
 #   3 "mean"
 
 # setup R error handling to go to stderr
-options(show.error.messages = F, error = function() {
-  cat(geterrmessage(), file = stderr()); q("no", 1, F)
+options(show.error.messages = FALSE, error = function() {
+  cat(geterrmessage(), file = stderr())
+  q("no", 1, FALSE)
 })
 
 # we need that to not crash galaxy with an UTF8 error on German LC settings.
@@ -69,7 +70,9 @@
   "outlier_filter_off", "b", 0, "logical",
   "auto_mean_filter_off", "c", 0, "logical",
   "beta_prior_off", "d", 0, "logical",
-  "alpha_ma", "A", 1, "numeric"
+  "alpha_ma", "A", 1, "numeric",
+  "prefilter", "P", 0, "logical",
+  "prefilter_value", "V", 1, "numeric"
 ), byrow = TRUE, ncol = 4)
 opt <- getopt(spec)
 
@@ -239,7 +242,7 @@
             size_factors <- estimateSizeFactorsForMatrix(counts(dds))
         }
     }
-    write.table(size_factors, file = opt$sizefactorsfile, sep = "\t", col.names = F, quote = FALSE)
+    write.table(size_factors, file = opt$sizefactorsfile, sep = "\t", col.names = FALSE, quote = FALSE)
 }
 
 apply_batch_factors <- function(dds, batch_factors) {
@@ -253,7 +256,7 @@
   dds_data <- colData(dds)
   # Merge dds_data with batch_factors using indexes, which are sample names
   # Set sort to False, which maintains the order in dds_data
-  reordered_batch <- merge(dds_data, batch_factors, by.x = 0, by.y = 0, sort = F)
+  reordered_batch <- merge(dds_data, batch_factors, by.x = 0, by.y = 0, sort = FALSE)
   batch_factors <- reordered_batch[, ncol(dds_data):ncol(reordered_batch)]
   for (factor in colnames(batch_factors)) {
     dds[[factor]] <- batch_factors[[factor]]
@@ -263,7 +266,7 @@
 }
 
 if (!is.null(opt$batch_factors)) {
-  batch_factors <- read.table(opt$batch_factors, sep = "\t", header = T)
+  batch_factors <- read.table(opt$batch_factors, sep = "\t", header = TRUE)
   dds <- apply_batch_factors(dds = dds, batch_factors = batch_factors)
   batch_design <- colnames(batch_factors)[-c(1, 2)]
   design_formula <- as.formula(paste("~", paste(c(batch_design, rev(factors)), collapse = " + ")))
@@ -280,6 +283,12 @@
   cat(paste(ncol(dds), "samples with counts over", nrow(dds), "genes\n"))
 }
 
+# minimal pre-filtering
+if (!is.null(opt$prefilter)) {
+    keep <- rowSums(counts(dds)) >= opt$prefilter_value
+    dds <- dds[keep, ]
+}
+
 # optional outlier behavior
 if (is.null(opt$outlier_replace_off)) {
   min_rep <- 7