Merge pull request #450 from atrigila/add_complex_contrasts

feat: allow usage of strings for `makeContrasts` in `DREAM`

Author: atrigila

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- [[#450](https://github.com/nf-core/differentialabundance/pull/441)] - Allow usage of strings for makeContrasts in DREAM. ([@atrigila](https://github.com/atrigila), review by [@pinin4fjords](https://github.com/pinin4fjords), [@suzannejin](https://github.com/suzannejin) and [@grst](https://github.com/grst)).
 - [[#441](https://github.com/nf-core/differentialabundance/pull/441)] - Add dream differential module. ([@nschcolnicov](https://github.com/nschcolnicov) and [@alanmmobbs93](https://github.com/alanmobbs93), review by [@pinin4fjords](https://github.com/pinin4fjords), [@suzannejin](https://github.com/suzannejin) and [@grst](https://github.com/grst)).
 - [[#440](https://github.com/nf-core/differentialabundance/pull/440)] - Add handling for formula based contrasts. ([@nschcolnicov](https://github.com/nschcolnicov), review by [@pinin4fjords](https://github.com/pinin4fjords))
 - [[#437](https://github.com/nf-core/differentialabundance/pull/437)] - Add nf-tests to deseq2/differential, rmarkdownnotebook, and proteus modules. ([@nschcolnicov](https://github.com/nschcolnicov), review by [@TODO](TODO)).

assets/schema_contrasts.json

@@ -15,7 +15,7 @@
                     },
                     "formula": {
                         "type": "string",
-                        "pattern": "^~\\s*[a-zA-Z_][a-zA-Z0-9_]*(\\s*\\+\\s*[a-zA-Z_][a-zA-Z0-9_]*)*$"
+                        "pattern": "^~\\s*[a-zA-Z_][a-zA-Z0-9_]*(\\s*([:+*])\\s*[a-zA-Z_][a-zA-Z0-9_]*)*$"
                     },
                     "comparison": {
                         "type": "array",
@@ -31,9 +31,12 @@
                         },
                         "minItems": 1,
                         "uniqueItems": true
+                    },
+                    "make_contrasts_str": {
+                        "type": "string"
                     }
                 },
-                "required": ["id", "comparison"],
+                "required": ["id"],
                 "additionalProperties": false
             },
             "minItems": 1

conf/test.config

@@ -20,7 +20,7 @@ params {
     input                    = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.samplesheet.csv'
     matrix                   = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.salmon.merged.gene_counts.top1000cov.tsv'
     transcript_length_matrix = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.spoofed_lengths.tsv'
-    contrasts_yml            = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/SRP254919.contrasts.yaml'
+    contrasts_yml            = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/SRP254919.contrasts.yaml'
 
     // To do: replace this with a cut-down mouse GTF matching the matrix for testing
     gtf = 'https://ftp.ensembl.org/pub/release-81/gtf/mus_musculus/Mus_musculus.GRCm38.81.gtf.gz'

conf/test_affy.config

@@ -20,7 +20,7 @@ params {
 
     // Input data
     input                  = 'https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/GSE50790.csv'
-    contrasts_yml          = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/GSE50790_contrasts.yaml'
+    contrasts_yml          = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/GSE50790_contrasts.yaml'
     affy_cel_files_archive = 'https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/GSE50790_RAW.tar'
 
     // Observations

conf/test_full.config

@@ -16,7 +16,7 @@ params {
 
     // Input data
     input         = 'https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/rnaseq_featurecounts_sample_preparations.tsv'
-    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/rnaseq_featurecounts_contrast_file.yaml'
+    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/rnaseq_featurecounts_contrast_file.yaml'
     matrix        = 'https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/rnaseq_featurecounts_merged_gene_counts.tsv'
     gtf           = 'https://ftp.ensembl.org/pub/release-81/gtf/mus_musculus/Mus_musculus.GRCm38.81.gtf.gz'
 

conf/test_maxquant.config

@@ -21,7 +21,7 @@ params {
     // Input data
     input         = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/proteomics/maxquant/MaxQuant_samplesheet.tsv'
     matrix        = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/proteomics/maxquant/MaxQuant_proteinGroups.txt'
-    contrasts_yml = 'https://github.com/nf-core/test-datasets/raw/refs/heads/differentialabundance/testdata/formula_contrasts/MaxQuant_contrasts.yaml'
+    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/MaxQuant_contrasts.yaml'
 
     // Observations
     observations_id_col   = 'Experiment'

conf/test_nogtf.config

@@ -22,7 +22,7 @@ params {
 
     input         = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.samplesheet.csv'
     matrix        = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.salmon.merged.gene_counts.top1000cov.tsv'
-    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/SRP254919.contrasts.yaml'
+    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/SRP254919.contrasts.yaml'
 
     //Features
     features_metadata_cols = 'gene_id,gene_name'

conf/test_rnaseq_limma.config

@@ -22,7 +22,7 @@ params {
     input                    = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.samplesheet.csv'
     matrix                   = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.salmon.merged.gene_counts.top1000cov.tsv'
     transcript_length_matrix = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.spoofed_lengths.tsv'
-    contrasts_yml            = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/SRP254919.contrasts.yaml'
+    contrasts_yml            = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/SRP254919.contrasts.yaml'
 
     // To do: replace this with a cut-down mouse GTF matching the matrix for testing
     gtf = 'https://ftp.ensembl.org/pub/release-81/gtf/mus_musculus/Mus_musculus.GRCm38.81.gtf.gz'

conf/test_soft.config

@@ -20,7 +20,7 @@ params {
 
     // Input
     input         = 'https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/GSE50790.csv'
-    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/formula_contrasts/GSE50790_contrasts.yaml'
+    contrasts_yml = 'https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/differentialabundance/testdata/GSE50790_contrasts.yaml'
     querygse      = 'GSE50790'
 
 }

docs/usage.md

@@ -175,30 +175,58 @@ contrasts:
     blocking_factors: ["replicate"]
 ```
 
-Alternatively, the YAML contrasts also supports formula based model definitions:
+The necessary fields in order are:
+
+- `id` - an arbitrary identifier, will be used to name contrast-wise output files
+- `comparison`(respectively):
+- - `variable` - which column from the observations information will be used to define groups
+- - `reference` - the base/ reference level for the comparison. If features have higher values in this group than target they will generate negative fold changes
+- - `target` - the target/ non-reference level for the comparison. If features have higher values in this group than the reference they will generate positive fold changes
+- `blocking_factors` - Any additional variables (also observation columns) that should be modelled alongside the contrast variable
+- `exclude_samples_col` and `exclude_samples_values` - the former being a valid column in the samples sheet, the latter a list of values in that column which should be used to select samples prior to differential modelling. This is helpful where certain samples need to be excluded prior to analysis of a given contrast.
+
+Alternatively, the YAML contrasts also supports formula based model definitions for tools such as `VARIANCEPARTITION_DREAM`:
 
 ```yaml
 contrasts:
   - id: condition_control_treated
     formula: "~ condition"
-    comparison: ["condition", "control", "treated"]
+    make_contrasts_str: "conditiontreated"
   - id: condition_control_treated_blockrep
     formula: "~ condition + replicate"
-    comparison: ["condition", "control", "treated"]
+    make_contrasts_str: "conditiontreated"
 ```
 
 The necessary fields in order are:
 
-- `id` - an arbitrary identifier, will be used to name contrast-wise output files
-- `comparison`(respectively):
-- - `variable` - which column from the observations information will be used to define groups
-- - `reference` - the base/ reference level for the comparison. If features have higher values in this group than target they will generate negative fold changes
-- - `target` - the target/ non-reference level for the comparison. If features have higher values in this group than the reference they will generate positive fold changes
+- `formula` - A string representation of the model formula. It is used to build the design matrix.
+- `make_contrasts_str` - An explicit literal contrast string (e.g., "treatmenthND6 - treatmentmCherry") that is passed directly to [`limma::makeContrasts()`](https://rdrr.io/bioc/limma/man/makeContrasts.html) in `VARIANCEPARTITION_DREAM`. The parameter names must be syntactically valid variable names in R (see [`make.names`](https://stat.ethz.ch/R-manual/R-devel/library/base/html/make.names.html)). This field provides full control for complex designs. Requires `formula`.
 
-You can optionally supply:
+> [!WARNING]
+> Formula-based contrasts are currently only supported by `VARIANCEPARTITION_DREAM`. They **do not work** with tools like `DESEQ2` or base `LIMMA` implementations yet.
 
-- `blocking_factors` - Any additional variables (also observation columns) that should be modelled alongside the contrast variable
-- `exclude_samples_col` and `exclude_samples_values` - the former being a valid column in the samples sheet, the latter a list of values in that column which should be used to select samples prior to differential modelling. This is helpful where certain samples need to be excluded prior to analysis of a given contrast.
+> [!NOTE]
+>
+> #### Notes on `make_contrasts_str`
+>
+> This string must match exactly the name of the coefficient in the model matrix as generated by the specified `formula`. It is passed to `limma::makeContrasts()` without modification. For example:
+>
+> - `formula: "~ condition"` will generate model coefficients like `conditiontreated` (if `control` is the reference).
+> - Then, `make_contrasts_str: "conditiontreated"` selects that coefficient for testing.
+>
+> This gives full control over the contrast definition but requires understanding of the model matrix.
+
+Beyond the basic one-factor comparison, the YAML contrasts format supports advanced experimental designs through the use of interaction terms and custom contrast strings. These are particularly useful in multifactorial experiments where the effect of one variable may depend on the level of another (e.g. genotype × treatment). To model an interaction between genotype and treatment, use a formula like `~ genotype * treatment`, which expands the yaml to:
+
+```yaml
+contrasts:
+  - id: genotype_WT_KO_treatment_Control_Treated
+    formula: "~ genotype * treatment"
+    comparison: ["genotype", "WT", "KO"]
+    make_contrasts_str: "genotypeKO.treatmentTreated"
+```
+
+To facilitate constructing and validating such models and contrast strings, consider using the [`ExploreModelMatrix`](https://www.bioconductor.org/packages/release/bioc/html/ExploreModelMatrix.html) Shiny app to have visual inspection of the design matrix and interactive contrast building. Another helpful resource is the [guide to creating design matrices for gene expression experiments](https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/designmatrices.html).
 
 ## Feature annotations