Mercurial > repos > iuc > deseq2

comparison deseq2.R @ 31:9a882d108833 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/deseq2 commit 469558ddf5bc6249874fe5826637fd6ee81588cf
author iuc
date Tue, 18 Jul 2023 14:58:52 +0000
parents 8fe98f7094de
children
comparison
equal deleted inserted replaced
30:8fe98f7094de 31:9a882d108833
33 # setup R error handling to go to stderr 33 # setup R error handling to go to stderr
34 options(show.error.messages = FALSE, error = function() { 34 options(show.error.messages = FALSE, error = function() {
35 cat(geterrmessage(), file = stderr()) 35 cat(geterrmessage(), file = stderr())
36 q("no", 1, FALSE) 36 q("no", 1, FALSE)
37 }) 37 })
38
39 # we need that to not crash galaxy with an UTF8 error on German LC settings.
40 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
41 38
42 library("getopt") 39 library("getopt")
43 library("tools") 40 library("tools")
44 options(stringAsFactors = FALSE, useFancyQuotes = FALSE) 41 options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
45 args <- commandArgs(trailingOnly = TRUE) 42 args <- commandArgs(trailingOnly = TRUE)
67 "fit_type", "t", 1, "integer", 64 "fit_type", "t", 1, "integer",
68 "many_contrasts", "m", 0, "logical", 65 "many_contrasts", "m", 0, "logical",
69 "outlier_replace_off", "a", 0, "logical", 66 "outlier_replace_off", "a", 0, "logical",
70 "outlier_filter_off", "b", 0, "logical", 67 "outlier_filter_off", "b", 0, "logical",
71 "auto_mean_filter_off", "c", 0, "logical", 68 "auto_mean_filter_off", "c", 0, "logical",
72 "beta_prior_off", "d", 0, "logical", 69 "use_beta_priors", "d", 0, "logical",
73 "alpha_ma", "A", 1, "numeric", 70 "alpha_ma", "A", 1, "numeric",
74 "prefilter", "P", 0, "logical", 71 "prefilter", "P", 0, "logical",
75 "prefilter_value", "V", 1, "numeric" 72 "prefilter_value", "V", 1, "numeric"
76 ), byrow = TRUE, ncol = 4) 73 ), byrow = TRUE, ncol = 4)
77 opt <- getopt(spec) 74 opt <- getopt(spec)
215 } 212 }
216 cat("\n---------------------\n") 213 cat("\n---------------------\n")
217 } 214 }
218 215
219 dds <- get_deseq_dataset(sample_table, header = opt$header, design_formula = design_formula, tximport = opt$tximport, txtype = opt$txtype, tx2gene = opt$tx2gene) 216 dds <- get_deseq_dataset(sample_table, header = opt$header, design_formula = design_formula, tximport = opt$tximport, txtype = opt$txtype, tx2gene = opt$tx2gene)
220 # estimate size factors for the chosen method 217
218 # use/estimate size factors with the chosen method
221 if (!is.null(opt$esf)) { 219 if (!is.null(opt$esf)) {
222 dds <- estimateSizeFactors(dds, type = opt$esf) 220 if (opt$esf %in% list("ratio", "poscounts", "iterate")) {
223 } 221 cat("Calculating size factors de novo\n")
222 dds <- estimateSizeFactors(dds, type = opt$esf)
223 } else {
224 sf_table <- read.table(opt$esf)
225 # Sort the provided size factors just in case the order differs from the input file order.
226 merged_table <- merge(sample_table, sf_table, by.x = 0, by.y = 1, sort = FALSE)
227 sf_values <- as.numeric(unlist(merged_table[5]))
228 "sizeFactors"(dds) <- sf_values
229
230 cat("Using user-provided size factors:\n")
231 print(sf_values)
232 }
233 } else {
234 cat("No size factor was used\n")
235 }
236
224 237
225 # estimate size factors for each sample 238 # estimate size factors for each sample
226 # - https://support.bioconductor.org/p/97676/ 239 # - https://support.bioconductor.org/p/97676/
227 if (!is.null(opt$sizefactorsfile)) { 240 if (!is.null(opt$sizefactorsfile)) {
228 nm <- assays(dds)[["avgTxLength"]] 241 nm <- assays(dds)[["avgTxLength"]]
310 independent_filtering <- FALSE 323 independent_filtering <- FALSE
311 if (verbose) cat("automatic filtering on the mean off\n") 324 if (verbose) cat("automatic filtering on the mean off\n")
312 } 325 }
313 326
314 # shrinkage of LFCs 327 # shrinkage of LFCs
315 if (is.null(opt$beta_prior_off)) { 328 if (is.null(opt$use_beta_priors)) {
316 beta_prior <- TRUE
317 } else {
318 beta_prior <- FALSE 329 beta_prior <- FALSE
319 if (verbose) cat("beta prior off\n") 330 if (verbose)
320 } 331 cat("Applied default - beta prior off\n")
332 } else {
333 beta_prior <- opt$use_beta_priors
334 }
335 sprintf("use_beta_prior is set to %s", beta_prior)
336
321 337
322 # dispersion fit type 338 # dispersion fit type
323 if (is.null(opt$fit_type)) { 339 if (is.null(opt$fit_type)) {
324 fit_type <- "parametric" 340 fit_type <- "parametric"
325 } else { 341 } else {