Enable complex contrast strings in DESeq2

modules/nf-core/deseq2/differential/environment.yml

@@ -5,3 +5,4 @@ channels:
   - bioconda
 dependencies:
   - bioconda::bioconductor-deseq2=1.34.0
+  - bioconda::bioconductor-limma=3.50.0

modules/nf-core/deseq2/differential/main.nf

@@ -4,8 +4,8 @@ process DESEQ2_DIFFERENTIAL {
 
     conda "${moduleDir}/environment.yml"
     container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
-        'https://depot.galaxyproject.org/singularity/bioconductor-deseq2:1.34.0--r41hc247a5b_3' :
-        'biocontainers/bioconductor-deseq2:1.34.0--r41hc247a5b_3' }"
+        'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/a1/a15f5d61792b60b6179afd885db27d3fe60eb4c42e805c8887ed0416d88cb484/data' :
+        'community.wave.seqera.io/library/bioconductor-deseq2_bioconductor-limma:b56a0c9ddc3e87e1' }"
 
     input:
     tuple val(meta), val(contrast_variable), val(reference), val(target), val(formula), val(comparison)

modules/nf-core/deseq2/differential/templates/deseq2_differential.R

@@ -231,6 +231,7 @@ for (file_input in c('count_file', 'sample_file')){
 
 library(DESeq2)
 library(BiocParallel)
+library(limma)
 
 ################################################
 ################################################
@@ -425,17 +426,25 @@ dds <- DESeq(
 
 if (!is.null(opt\$contrast_string)) {
   coef_names <- resultsNames(dds)
-  if (!(opt\$contrast_string %in% coef_names)) {
-    stop(sprintf(
-      "Contrast '%s' not in design. Available coefficients: %s",
-      opt\$contrast_string,
-      paste(coef_names, collapse = ", ")
-    ))
-  }
+  if (opt\$contrast_string %in% coef_names) {
+    # Direct coefficient name
+    comp.results <- run_results(name = opt\$contrast_string)
+    if (opt\$shrink_lfc) {
+      comp.results <- run_shrink(coef = opt\$contrast_string)
+    }
+  } else {
+    # Parse as limma-style contrast expression
+    design_mat <- model.matrix(as.formula(model), data = as.data.frame(colData(dds)))
+    colnames(design_mat) <- make.names(colnames(design_mat))
+    numeric_contrast <- as.numeric(
+      limma::makeContrasts(contrasts = opt\$contrast_string, levels = colnames(design_mat))
+    )
 
-  comp.results <- run_results(name = opt\$contrast_string)
-  if (opt\$shrink_lfc) {
-    comp.results <- run_shrink(coef = opt\$contrast_string)
+    # Run DESeq2 results with numeric contrast
+    comp.results <- run_results(contrast = numeric_contrast)
+    if (opt\$shrink_lfc) {
+      comp.results <- run_shrink(contrast = numeric_contrast)
+    }
   }
 } else {
   contrast_var_tg_ref <- c(contrast_variable,
@@ -452,7 +461,7 @@ if (!is.null(opt\$contrast_string)) {
 
 if (!is.null(opt\$transcript_lengths_file)){
   size_factors = estimateSizeFactorsForMatrix(counts(dds) / assays(dds)[["avgTxLength"]])
-}else {
+} else {
   size_factors = sizeFactors(dds)
 }
 

modules/nf-core/deseq2/differential/tests/contrasts_shrink.config

@@ -0,0 +1,5 @@
+process {
+    withName: 'DESEQ2_DIFFERENTIAL' {
+        ext.args = { "--round_digits 2 --sample_id_col sample --vs_method rlog --shrink_lfc TRUE --seed 1" }
+    }
+}

modules/nf-core/deseq2/differential/tests/main.nf.test

@@ -57,6 +57,102 @@ nextflow_process {
         }
     }
 
+    test("RNAseq - Feature Counts - formula + limma contrast string - interaction") {
+        config './contrasts_interaction.config'
+
+        when {
+            process {
+                """
+                input[0] = Channel.of([
+                    'id': 'genotype_WT_KO',
+                    'formula': '~ 0 + genotype',
+                    'comparison': 'genotypeWT - genotypeKO'
+                ])
+                .map {
+                    tuple(it, it.variable, it.reference, it.target, it.formula, it.comparison)
+                }
+
+                input[1] = Channel.of([
+                    [ id: 'test' ],
+                    file("https://github.com/nf-core/test-datasets/raw/differentialabundance/modules_testdata/variancepartition_dream/metadata.tsv", checkIfExists: true),
+                    file("https://github.com/nf-core/test-datasets/raw/differentialabundance/modules_testdata/variancepartition_dream/counts.tsv", checkIfExists: true)
+                ])
+
+                ch_spikes = [[],[]]
+                ch_lengths = [[],[]]
+                input[2] = ch_spikes
+                input[3] = ch_lengths
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.success },
+                { assert snapshot(
+                        process.out.results,
+                        process.out.size_factors,
+                        process.out.normalised_counts,
+                        process.out.rlog_counts,
+                        process.out.vst_counts,
+                        process.out.model,
+                        process.out.versions,
+                        file(process.out.dispersion_plot[0][1]).name,
+                        file(process.out.rdata[0][1]).name
+                    ).match() },
+                { assert path(process.out.session_info[0][1]).text.contains("DESeq2_1.34.0") }
+            )
+        }
+    }
+
+    test("RNAseq - Feature Counts - formula + limma contrast string - shrunken") {
+        config './contrasts_shrink.config'
+
+        when {
+            process {
+                """
+                input[0] = Channel.of([
+                    'id': 'genotype_WT_KO',
+                    'formula': '~ 0 + genotype',
+                    'comparison': 'genotypeWT - genotypeKO'
+                ])
+                .map {
+                    tuple(it, it.variable, it.reference, it.target, it.formula, it.comparison)
+                }
+
+                input[1] = Channel.of([
+                    [ id: 'test' ],
+                    file("https://github.com/nf-core/test-datasets/raw/differentialabundance/modules_testdata/variancepartition_dream/metadata.tsv", checkIfExists: true),
+                    file("https://github.com/nf-core/test-datasets/raw/differentialabundance/modules_testdata/variancepartition_dream/counts.tsv", checkIfExists: true)
+                ])
+
+                ch_spikes = [[],[]]
+                ch_lengths = [[],[]]
+                input[2] = ch_spikes
+                input[3] = ch_lengths
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.success },
+                { assert snapshot(
+                        process.out.results,
+                        process.out.size_factors,
+                        process.out.normalised_counts,
+                        process.out.rlog_counts,
+                        process.out.vst_counts,
+                        process.out.model,
+                        process.out.versions,
+                        file(process.out.dispersion_plot[0][1]).name,
+                        file(process.out.rdata[0][1]).name
+                    ).match() },
+                { assert path(process.out.session_info[0][1]).text.contains("DESeq2_1.34.0") }
+            )
+        }
+    }
+
     test("mouse - contrasts - matrix") {
 
         config './contrasts_matrix.config'