Add limma for rnaseq

.github/workflows/ci.yml

@@ -35,6 +35,7 @@ jobs:
           - "test_affy"
           - "test_maxquant"
           - "test_soft"
+          - "test_rnaseq_limma"
         compute_profile:
           - "conda"
           - "docker"

CHANGELOG.md

@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- [[#286](https://github.com/nf-core/differentialabundance/pull/286)] - Integration of limma voom for rnaseq data ([@KamilMaliszArdigen](https://github.com/KamilMaliszArdigen), review by [@pinin4fjords](https://github.com/pinin4fjords))
+
 ### Fixed
 
 - [[#304](https://github.com/nf-core/differentialabundance/pull/304)] - Removed TXT file options from nextflow_schema where they are equivalent to TSV to make the input files clearer ([@WackerO](https://github.com/WackerO), review by [@pinin4fjords](https://github.com/pinin4fjords))

README.md

@@ -42,7 +42,7 @@ On release, automated continuous integration tests run the pipeline on a full-si
 > [!NOTE]
 > If you are new to Nextflow and nf-core, please refer to [this page](https://nf-co.re/docs/usage/installation) on how to set-up Nextflow. Make sure to [test your setup](https://nf-co.re/docs/usage/introduction#how-to-run-a-pipeline) with `-profile test` before running the workflow on actual data.
 
-RNA-seq:
+RNA-seq with deseq2:
 
 ```bash
  nextflow run nf-core/differentialabundance \
@@ -60,6 +60,27 @@ If you are using the outputs of the nf-core rnaseq workflow as input here **eith
 - supply the raw count matrices (file names like **gene_counts.tsv**) alongide the transcript length matrix via `--transcript_length_matrix` (rnaseq versions >=3.12.0, preferred)
 - **or** supply the **gene_counts_length_scaled.tsv** or **gene_counts_scaled.tsv** matrices.
 
+RNA-seq limma+voom:
+
+```bash
+ nextflow run nf-core/differentialabundance \
+     --input samplesheet.csv \
+     --contrasts contrasts.csv \
+     --matrix assay_matrix.tsv \
+     --gtf mouse.gtf \
+     --outdir <OUTDIR>  \
+     -profile rnaseq_limma,<docker/singularity/podman/shifter/charliecloud/conda/institute>
+```
+
+:::note
+If you are using the outputs of the nf-core rnaseq workflow as input here **either**:
+
+Provide either the **gene_counts_length_scaled.tsv** or **gene_counts_scaled.tsv** matrices. This follows the [recommendation from the tximport documentation](https://bioconductor.org/packages/devel/bioc/vignettes/tximport/inst/doc/tximport.html#limma-voom):
+
+> "Because limma-voom does not use the offset matrix stored in `y$offset`, we recommend using scaled counts generated from abundances, either 'scaledTPM' or 'lengthScaledTPM'."
+
+These matrices, **gene_counts_length_scaled.tsv** or **gene_counts_scaled.tsv**, are generated in the RNA-seq workflow and meet this recommendation by providing appropriately scaled counts for analysis.
+
 See the [usage documentation](https://nf-co.re/differentialabundance/usage) for more information.
 :::
 

conf/modules.config

@@ -274,6 +274,7 @@ process {
             "--lfc ${params.limma_lfc}",
             "--confint ${params.limma_confint}",
             "--subset_to_contrast_samples $params.differential_subset_to_contrast_samples",
+            "--use_voom \"${params.limma_use_voom}\"",
             ((meta.blocking == null) ? '' : "--blocking_variables $meta.blocking"),
             ((meta.exclude_samples_col == null) ? '' : "--exclude_samples_col $meta.exclude_samples_col"),
             ((meta.exclude_samples_values == null) ? '' : "--exclude_samples_values $meta.exclude_samples_values")
@@ -437,6 +438,7 @@ process {
                 pattern: '*.html'
             ]
         ]
+        memory = { 12.GB * task.attempt }
     }
 
     withName: MAKE_REPORT_BUNDLE {

conf/rnaseq_limma.config

@@ -0,0 +1,36 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Nextflow config file for running RNA-seq analysis
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Defines settings specific to RNA-seq analysis
+
+    Use as follows:
+        nextflow run nf-core/differentialabundance -profile rnaseq_limma,<docker/singularity> --outdir <OUTDIR>
+
+----------------------------------------------------------------------------------------
+*/
+
+params {
+
+    config_profile_name        = 'RNA-seq profile with limma'
+    config_profile_description = 'Settings for RNA-seq analysis with limma'
+
+    // Study
+    study_type = 'rnaseq'
+    study_abundance_type = 'counts'
+
+    // Observations
+    observations_id_col        = 'sampleId'
+    observations_name_col      = 'sampleId'
+
+    // Differential options
+    differential_use_limma           = true
+    differential_file_suffix         = ".limma.results.tsv"
+    differential_fc_column           = "logFC"
+    differential_pval_column         = "P.Value"
+    differential_qval_column         = "adj.P.Val"
+    differential_feature_id_column   = "Geneid"
+    differential_feature_name_column = "Geneid"
+    limma_use_voom                   = true
+
+}

conf/test_rnaseq_limma.config

@@ -0,0 +1,67 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Nextflow config file for running RNA-seq analysis
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Defines settings specific to RNA-seq analysis
+
+    Use as follows:
+        nextflow run nf-core/differentialabundance -profile test_rnaseq_limma,<docker/singularity> --outdir <OUTDIR>
+
+----------------------------------------------------------------------------------------
+*/
+
+includeConfig 'rnaseq_limma.config'
+
+
+params {
+    study_name = 'SRP254919'
+    config_profile_name        = 'Test profile'
+    config_profile_description = 'Minimal test dataset to check pipeline function'
+
+    // Input data
+    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 = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.contrasts.csv'
+
+    // 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'
+
+    // Observations
+    observations_id_col        = 'sample'
+    observations_name_col      = 'sample'
+
+    // Features
+    features_type               = 'gene'
+    features_id_col             = 'gene_id'
+    features_name_col           = 'gene_name'
+    features_metadata_cols      = 'gene_id,gene_name,gene_biotype'
+
+    // Differential options
+    differential_feature_id_column   = "gene_id"
+    differential_feature_name_column = "gene_name"
+    differential_fc_column           = "logFC"
+
+    // Apply a higher filter to check that the filtering works
+    filtering_min_abundance=10
+
+    // Exploratory
+    exploratory_assay_names = "raw,normalised"
+    exploratory_final_assay = "normalised"
+    exploratory_log2_assays = 'raw,normalised'
+    exploratory_main_variable = 'contrasts'
+
+    // Test dataset is too small for the nsub default value
+    deseq2_vst_nsub                   = 500
+
+    // Activate GSEA
+    gsea_run = true
+    gene_sets_files = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/gene_set_analysis/mh.all.v2022.1.Mm.symbols.gmt'
+
+    // Report options
+    report_round_digits = 3
+    report_contributors = 'Jane Doe\nDirector of Institute of Microbiology\nUniversity of Smallville;John Smith\nPhD student\nInstitute of Microbiology\nUniversity of Smallville'
+
+}
+
+