Allow dream to process non-formula contrasts

[delfiterradas]

PR checklist

Closes nf-core/differentialabundance#578

• This comment contains a description of changes (with reason).
• If you've fixed a bug or added code that should be tested, add tests!
• If you've added a new tool - have you followed the module conventions in the contribution docs
• If necessary, include test data in your PR.
• Remove all TODO statements.
• Emit the versions.yml file.
• Follow the naming conventions.
• Follow the parameters requirements.
• Follow the input/output options guidelines.
• Add a resource label
• Use BioConda and BioContainers if possible to fulfil software requirements.
• Ensure that the test works with either Docker / Singularity. Conda CI tests can be quite flaky:
  • For modules:
    • nf-core modules test <MODULE> --profile docker
    • nf-core modules test <MODULE> --profile singularity
    • nf-core modules test <MODULE> --profile conda
  • For subworkflows:
    • nf-core subworkflows test <SUBWORKFLOW> --profile docker
    • nf-core subworkflows test <SUBWORKFLOW> --profile singularity
    • nf-core subworkflows test <SUBWORKFLOW> --profile conda

[grst]

I think you still need to take the baseline level into account when building contrasts, as done here for DESeq2:

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

Lines 450 to 457 in 8846461

	contrast_var_tg_ref <- c(contrast_variable,
	opt\$target_level,
	opt\$reference_level)
	comp.results <- run_results(contrast = contrast_var_tg_ref)
	if (opt\$shrink_lfc) {
	comp.results <- run_shrink(contrast = contrast_var_tg_ref)
	}

DESeq2 allows to test contrasts with a c(variable, target, reference) tuple. I don't know of DREAM offers the same. If it doesn't, you'll need to manually build the contrast or choose appropriate baseline levels for the factor with the contrast variable.

[Copilot]

Pull request overview

This pull request enables the dream module to handle contrasts without requiring an explicit formula. Previously, dream would fail if no formula was provided; now it automatically generates a formula from the contrast variable and optional blocking variables.

Key Changes:

• Modified the workflow filter to allow dream to process contrasts without formulas
• Added automatic formula generation logic when no formula is provided
• Added validation for contrast variables and blocking variables
• Updated tests to cover the no-formula scenario

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
subworkflows/nf-core/abundance_differential_filter/main.nf	Removed the formula requirement from the dream inputs filter, allowing processing of non-formula contrasts
modules/nf-core/variancepartition/dream/templates/dream.R	Added logic to automatically generate formulas from contrast and blocking variables when no formula is provided, including validation and variable checking
modules/nf-core/variancepartition/dream/tests/main.nf.test	Reorganized tests and added a new test case for matrix-based contrasts without formulas
modules/nf-core/variancepartition/dream/tests/main.nf.test.snap	Added snapshot for the new no-formula test case
modules/nf-core/variancepartition/dream/meta.yml	Updated documentation to indicate formula parameter is now optional instead of mandatory

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[grst]

Getting close!

Can you please add a check that (variable, reference, target) contrasts are not combined with a custom formula? I'd rather have the module fail that doing something incorrect.

[grst]

I think we can error in this case?

Perfect colinearity of the factors would show that it is an incorrect model.
Reference level as baseline can't happen if we always create a zero-intercept model.

[grst]

@pinin4fjords, I'm inclined to keep things simple and avoid possible error sources by making the following options mutually exclusive:

• specify a formula and a constrast_string
• specify a contrast tuple (variable, baseline, target) and blocking_factors.

I.e. I would raise an error if someone specifies a formula together with a contrast tuple. Does that sound alright to you?

[grst]

So regarding the eBayes, it is necessary to call it after contrasts.fit, according to this example from the limma user guide:

https://bioconductor.org/packages//2.11/bioc/vignettes/limma/inst/doc/usersguide.pdf

It wouldn't be necessary to call it before contrasts.fit, but it doesn't harm. (we still need it if we don't call contrasts.fit at all, but results() directly on the first fit).

[pinin4fjords]

I'll defer to @grst on approval for this

[grst]

ok, then here are my last two requests:

• Can you please add a check that (variable, reference, target) contrasts are not combined with a custom formula?
• remove the special case when factors are missing from the design matrix, just fail instead.

Like that we focus on the usecases we currently need for differentialabundance and keep it simple.

[delfiterradas]

@grst

ok, then here are my last two requests:

• Can you please add a check that (variable, reference, target) contrasts are not combined with a custom formula?
• remove the special case when factors are missing from the design matrix, just fail instead.

Like that we focus on the usecases we currently need for differentialabundance and keep it simple.

I will remove the special case now, but I am not sure how to approach the check for the contrast tuple since the script currently uses the contrast_variable and contrast_reference to relevel the metadata:

if (!is.null(opt\$contrast_reference) && opt\$contrast_reference != "") {
    metadata[[opt\$contrast_variable]] <- factor(metadata[[opt\$contrast_variable]])
    metadata[[opt\$contrast_variable]] <- relevel(metadata[[opt\$contrast_variable]], ref = opt\$contrast_reference)
}

I could make the check only for blocking_variables and contrast_reference but then that wouldn't allow subsetting to contrast samples.

Would the check be necessary if we already have a validation where contrast_string must be provided if using formula?

[grst]

I don't think the relevel or subset to contrast samples should work together with a custom contrast/formula. At the pipeline level, when specifying a formula and a contrast string, contrast variable, target and reference won't be set.

You can write the check such that this possibility is excluded.

[delfiterradas]

Okay, I have added the validation so that the module fails when a formula is used together with the contrast tuple.

Some of the tests are failing because no explicit reference level is set and therefore no releveling is performed, so sometimes the "target" level is absorbed into the intercept by the default:

    >   ~genotype * treatment
    >   Warning message:
    >   Zero sample variances detected, have been offset away from zero 
    >           (Intercept) genotypeWT treatmentTreated genotypeWT:treatmentTreated
    >   Sample1           1          1                0                           0
    >   Sample2           1          1                0                           0
    >   Sample3           1          1                0                           0
    >   [1] "(Intercept)"                 "genotypeWT"                 
    >   [3] "treatmentTreated"            "genotypeWT:treatmentTreated"
    >   Using contrast string: genotypeKO - treatmentTreated 
    >   Error in eval(ej, envir = levelsenv) : object 'genotypeKO' not found

In this case, genotype (KO) is encoded in the intercept by default and therefore does not exist as a standalone coefficient in the design matrix and cannot be directly referenced in a contrast expression.

@grst Do you think this is behavior that users are generally expected to understand and handle themselves?

[grst]

yes, users are expected to understand how this works if they use formulas and contrast strings.

One thing we should consider is how to explicitly set the baseline level of the factors (I believe currently it's the default, i.e. alphabetically). But that's a different topic that we should discuss at the pipeline level first.

The flagged issues have been fixed, but the "requested changes" status prevents us from merging. Using the "dimiss" feature to keep things going during the Christmas break.