Merge pull request #462 from atrigila/limma_formula

Author: atrigila

assets/differentialabundance_report.Rmd

@@ -367,7 +367,7 @@ if (is.null(differential_file_suffix)) {
     differential_file_suffix <- paste0(".", params$differential_method, ".results.tsv")
 }
 differential_files <- lapply(contrasts$id, function(d){
-    file.path(params$input_dir, paste0(gsub(' |;', '_', d), '_', params$study_name, differential_file_suffix))
+    file.path(params$input_dir, paste0(gsub(' |;', '_', d), differential_file_suffix))
 })
 differential_names <- paste0(contrasts$id, '_', params$study_name)
 
@@ -972,24 +972,54 @@ if (!is.null(params$functional_method)){
     cat("\n#### ", toupper(params$functional_method) ," {.tabset}\n")
 
     if (params$functional_method == 'gsea') {
+
+    # Only keep contrasts that have both reference and target
+    gsea_contrasts <- contrasts[
+        !is.na(contrasts$reference) & contrasts$reference != "" &
+        !is.na(contrasts$target) & contrasts$target != "",
+    ]
+
+    if (nrow(gsea_contrasts) == 0) {
+        warning("No contrasts with reference and target defined. Skipping GSEA report section.")
+    } else {
+
         for (gmt_file in simpleSplit(params$gene_sets_files)) {
-        gmt_name <- basename(tools::file_path_sans_ext(gmt_file))
-        cat("\n##### ", gmt_name ," {.tabset}\n")
-
-        reference_gsea_tables <- paste0(differential_names, ".", gmt_name, '.gsea_report_for_', contrasts$reference, '.tsv')
-        target_gsea_tables <- paste0(differential_names, ".", gmt_name, '.gsea_report_for_', contrasts$target, '.tsv')
-        for (i in 1:nrow(contrasts)){
-            cat("\n###### ", contrast_descriptions[i], "\n")
-            target_gsea_results <- read_metadata(target_gsea_tables[i])[,c(-2,-3)]
-            target_gsea_results <- round_dataframe_columns(target_gsea_results, digits=params$round_digits)
-            print( htmltools::tagList(datatable(target_gsea_results, caption = paste0("\nTarget (", contrasts$target[i], ")\n"), rownames = FALSE) ))
-            ref_gsea_results <- read_metadata(reference_gsea_tables[i])[,c(-2,-3)]
-            ref_gsea_results <- round_dataframe_columns(ref_gsea_results, digits=params$round_digits)
-            print( htmltools::tagList(datatable(ref_gsea_results, caption = paste0("\nReference (", contrasts$reference[i], ")\n"), rownames = FALSE) ))
+            gmt_name <- basename(tools::file_path_sans_ext(gmt_file))
+            cat("\n##### ", gmt_name ," {.tabset}\n")
+
+            reference_gsea_tables <- paste0(differential_names, ".", gmt_name, '.gsea_report_for_', gsea_contrasts$reference, '.tsv')
+            target_gsea_tables <- paste0(differential_names, ".", gmt_name, '.gsea_report_for_', gsea_contrasts$target, '.tsv')
+
+            for (i in seq_len(nrow(gsea_contrasts))) {
+                cat("\n###### ", contrast_descriptions[i], "\n")
+
+                if (file.exists(target_gsea_tables[i])) {
+                    target_gsea_results <- read_metadata(target_gsea_tables[i])[,c(-2,-3)]
+                    target_gsea_results <- round_dataframe_columns(target_gsea_results, digits=params$round_digits)
+                    print( htmltools::tagList(
+                        datatable(target_gsea_results,
+                            caption = paste0("\nTarget (", gsea_contrasts$target[i], ")\n"),
+                            rownames = FALSE)
+                    ))
+                } else {
+                    cat("\n*Target GSEA file missing: ", target_gsea_tables[i], "*\n")
+                }
+
+                if (file.exists(reference_gsea_tables[i])) {
+                    ref_gsea_results <- read_metadata(reference_gsea_tables[i])[,c(-2,-3)]
+                    ref_gsea_results <- round_dataframe_columns(ref_gsea_results, digits=params$round_digits)
+                    print( htmltools::tagList(
+                        datatable(ref_gsea_results,
+                            caption = paste0("\nReference (", gsea_contrasts$reference[i], ")\n"),
+                            rownames = FALSE)
+                    ))
+                } else {
+                    cat("\n*Reference GSEA file missing: ", reference_gsea_tables[i], "*\n")
+                }
+            }
         }
     }
-
-    } else if (params$functional_method == 'gprofiler2') {
+    }  else if (params$functional_method == 'gprofiler2') {
 
         cat(paste0("\nThis section contains the results tables of the pathway analysis which was done with the R package gprofiler2. The differential fraction is the number of differential genes in a pathway divided by that pathway's size, i.e. the number of genes annotated for the pathway.",
         ifelse(params$gprofiler2_significant, paste0(" Enrichment was only considered if significant, i.e. adjusted p-value <= ", params$gprofiler2_max_qval, "."), "Enrichment was also considered if not significant."), "\n"))

docs/usage.md

@@ -200,10 +200,10 @@ contrasts:
 The necessary fields in order are:
 
 - `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`.
+- `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` and `LIMMA_DIFFERENTIAL`. 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`.
 
 > [!WARNING]
-> Formula-based contrasts are currently only supported by `VARIANCEPARTITION_DREAM`. They **do not work** with tools like `DESEQ2` or base `LIMMA` implementations yet.
+> Formula-based contrasts are currently only supported by `VARIANCEPARTITION_DREAM` and `LIMMA_DIFFERENTIAL`. They **do not work** with tools like `DESEQ2` yet.
 
 > [!NOTE]
 >
@@ -215,6 +215,7 @@ The necessary fields in order are:
 > - 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.
+> Some downstream applications (e.g. `GSEA_GSEA`, `SHINYNGS_APP`) do not support formula-based contrasts as they require a `meta.variable`.
 
 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:
 
@@ -465,6 +466,8 @@ process {
 
 You will not get the final reporting outcomes of the workflow, but you will get the differential tables produced by DESeq2 or Limma, and the results of any gene seta analysis you have enabled.
 
+We have also added a dedicated pipeline parameter, `--skip_reports` that allows you to skip only the RMarkdown notebook and bundled report while leaving other reporting processes active. The `RMARKDOWNNOTEBOOK` process assumes that every grouping variable you pass to it (from the contrasts file’s variable column or PCA-derived informative_variables) exists as a valid, named column in your sample metadata. If you know your metadata or contrasts might be incomplete or non-standard (such as using formula-based yaml files), the you can use this flag to skip these steps.
+
 #### Restricting samples considered by DESeq2 or Limma
 
 By default, the DESeq2 or Limma differential modules model all samples at once, rather than just the samples involved in the contrast. This is usually the correct thing to do, but when there are are large numbers of samples involved in each contrast it may be unnecessary, and things can be sped up significantly by setting `--differential_subset_to_contrast_samples`. This will remove any samples not relevant to the contrast before the main differential analysis routines are called.

modules.json

@@ -72,7 +72,7 @@
                     },
                     "limma/differential": {
                         "branch": "master",
-                        "git_sha": "6ac0725bddb061e6f8eb3e22d20129f901021d9a",
+                        "git_sha": "41d1a6ad16756e12a424e508cbe3a1b3b8639946",
                         "installed_by": ["abundance_differential_filter"]
                     },
                     "propr/grea": {
@@ -82,7 +82,7 @@
                     },
                     "propr/propd": {
                         "branch": "master",
-                        "git_sha": "25c228dd132d91596e87b636b0e03fd48d4c951a",
+                        "git_sha": "ce6e8e8bc0c00d284733d7c48ddc794873f60ffe",
                         "installed_by": ["abundance_differential_filter"]
                     },
                     "proteus/readproteingroups": {
@@ -136,12 +136,12 @@
                 "nf-core": {
                     "abundance_differential_filter": {
                         "branch": "master",
-                        "git_sha": "6ac0725bddb061e6f8eb3e22d20129f901021d9a",
+                        "git_sha": "4006954b7bfbbc61c7d5fb009319ed536d2ffca4",
                         "installed_by": ["subworkflows"]
                     },
                     "differential_functional_enrichment": {
                         "branch": "master",
-                        "git_sha": "6ac0725bddb061e6f8eb3e22d20129f901021d9a",
+                        "git_sha": "4006954b7bfbbc61c7d5fb009319ed536d2ffca4",
                         "installed_by": ["subworkflows"]
                     },
                     "utils_nextflow_pipeline": {