Merge pull request #443 from nf-core/add_toolsheet

Add toolsheet-related implementations

Author: suzannejin

.github/workflows/ci.yml

@@ -44,6 +44,8 @@ jobs:
           - "maxquant"
           - "soft"
           - "rnaseq_limma"
+          - "rnaseq_dream"
+          - "custom_paramsheet"
         compute_profile:
           - "docker"
           - "singularity"

.nf-core.yml

@@ -1,6 +1,8 @@
 lint:
   files_exist:
     - assets/multiqc_config.yml
+  files_unchanged:
+    - .github/CONTRIBUTING.md
   multiqc_config: false
   nextflow_config:
     - config_defaults:

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- [[#443](https://github.com/nf-core/differentialabundance/pull/443)] - Add toolsheet-related implementations. ([@suzannejin](https://github.com/suzannejin), review by [@pinin4fjords](https://github.com/pinin4fjords), [@mirpedrol](https://github.com/mirpedrol) and [@JoseEspinosa](https://github.com/JoseEspinosa))
 - [[#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))
@@ -38,6 +39,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+- [[#448](https://github.com/nf-core/differentialabundance/pull/448)] - Simplify toolsheet handling and restructure workflow to use paramset in meta. ([@pinin4fjords](https://github.com/pinin4fjords), review by [@suzannejin](https://github.com/suzannejin) and [@grst](https://github.com/grst))
 - [[#431](https://github.com/nf-core/differentialabundance/pull/431)] - Replace the calls to differential and functional analysis modules by subworkflows. ([@suzannejin](https://github.com/suzannejin), review by [@pinin4fjords](https://github.com/pinin4fjords))
 - [[#410](https://github.com/nf-core/differentialabundance/pull/410)] - Update contrasts file format to allow yaml ([@nschcolnicov](https://github.com/nschcolnicov), review by [@pinin4fjords](https://github.com/pinin4fjords)).
 - [[#374](https://github.com/nf-core/differentialabundance/pull/374)] - Update all modules and subworkflows ([@nschcolnicov](https://github.com/nschcolnicov), review by [@pinin4fjords](https://github.com/pinin4fjords)).

README.md

@@ -95,6 +95,21 @@ Affymetrix microarray:
 > [!WARNING]
 > Please provide pipeline parameters via the CLI or Nextflow `-params-file` option. Custom config files including those provided by the `-c` Nextflow option can be used to provide any configuration _**except for parameters**_; see [docs](https://nf-co.re/docs/usage/getting_started/configuration#custom-configuration-files).
 
+The paramsheet file (ie. `paramsheet.csv`) stored in the `assets` directory defines the combination of tools and parameters that make sense to run for a given study type. You can use the flag `--paramset_name` to specify which set of tools to run. For example:
+
+```bash
+ nextflow run nf-core/differentialabundance \
+     --input samplesheet.csv \
+     --contrasts contrasts.yaml \
+     --matrix assay_matrix.tsv \
+     --gtf mouse.gtf \
+     --outdir <OUTDIR>  \
+     -profile rnaseq,<docker/singularity/podman/shifter/charliecloud/conda/institute> \
+     --paramset_name deseq2_rnaseq_gprofiler2
+```
+
+You could also provide your own paramsheet through the `--paramsheet` parameter.
+
 For more details and further functionality, please refer to the [usage documentation](https://nf-co.re/differentialabundance/usage) and the [parameter documentation](https://nf-co.re/differentialabundance/parameters).
 
 ### Reporting

assets/differentialabundance_report.Rmd

@@ -26,7 +26,7 @@ params:
     report_author: NULL,
     report_description: NULL,
     report_scree: NULL
-    report_round_digits: NULL
+    round_digits: NULL
     observations_type: NULL
     observations: NULL                                          # GSE156533.samplesheet.csv
     observations_id_col: NULL
@@ -90,6 +90,7 @@ params:
     filtering_min_abundance: 1
     filtering_min_proportion: NULL
     filtering_grouping_var: NULL
+    differential_method: NULL
     differential_file_suffix: NULL
     differential_feature_id_column: NULL
     differential_feature_name_column: NULL
@@ -118,7 +119,7 @@ params:
     deseq2_cores: NULL
     deseq2_vs_blind: NULL
     deseq2_vst_nsub: NULL
-    gsea_run: false
+    functional_method: NULL
     gsea_nperm: NULL
     gsea_permute: NULL
     gsea_scoring_scheme: NULL
@@ -137,7 +138,6 @@ params:
     gsea_save_rnd_lists: NULL
     gsea_zip_report: NULL
     gsea_chip_file: NULL
-    gprofiler2_run: false
     gprofiler2_organism: NULL
     gprofiler2_significant: NULL
     gprofiler2_measure_underrepresentation: NULL
@@ -176,7 +176,7 @@ round_dataframe_columns <- function(df, columns = NULL, digits = -1) {
     if (is.null(columns)) {
         columns <- colnames(df)[(unlist(lapply(df, is.numeric), use.names=F))]    # extract only numeric columns for rounding
     }
-    df[,columns] <- format(data.frame(df[, columns], check.names = FALSE), scientific=T, digits=params$report_round_digits)
+    df[,columns] <- format(data.frame(df[, columns], check.names = FALSE), scientific=T, digits=params$round_digits)
     # Convert columns back to numeric
 
     for (c in columns) {
@@ -364,14 +364,12 @@ treatment-mCherry-hND6-batcheffect.deseq2.results.tsv
 ```{r, echo=FALSE}
 differential_file_suffix <- params$differential_file_suffix
 if (is.null(differential_file_suffix)) {
-    differential_file_suffix <- ifelse(params$study_type %in% c('rnaseq'), ".deseq2.results.tsv", ".limma.results.tsv")
+    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))
 })
-differential_names <- lapply(contrasts$id, function(d){
-    paste0(gsub(' |;', '_', d), '_', params$study_name)
-})
+differential_names <- paste0(contrasts$id, '_', params$study_name)
 
 # Initialize vector to store warning messages before merging tables
 warnings_list <- c()
@@ -535,7 +533,7 @@ contrasts_to_print <- contrasts
 colnames(contrasts_to_print) <- prettifyVariablename(colnames(contrasts_to_print))
 
 # Add design/model formulae to report
-de_tool <- ifelse(params$study_type %in% c('rnaseq'), "deseq2", "limma")
+de_tool <- params$differential_method
 contrasts_to_print$model <- sapply(contrasts_to_print$Id, function(id) {
     model_file <- paste0(id, ".", de_tool, ".model.txt")
     if (file.exists(model_file)) {
@@ -934,7 +932,7 @@ for (i in 1:nrow(contrasts)){
         colnames(contrast_de) <- prettifyVariablename(colnames(contrast_de))
 
         if (nrow(contrast_de) > 0){
-            contrast_de <- round_dataframe_columns(contrast_de, digits=params$report_round_digits)
+            contrast_de <- round_dataframe_columns(contrast_de, digits=params$round_digits)
             print( htmltools::tagList(datatable(contrast_de, caption = paste('Differential genes', dir, 'in', contrast_descriptions[i], " (check", differential_files[[i]], "for more detail)"), rownames = FALSE) ))
 
         if ("Gene biotype" %in% colnames(contrast_de)) {
@@ -969,56 +967,51 @@ for (i in 1:nrow(contrasts)){
 <!-- Gene set analysis results -->
 
 ```{r, echo=FALSE, results='asis'}
-possible_gene_set_methods <- c('gsea', 'gprofiler2')
-if (any(unlist(params[paste0(possible_gene_set_methods, '_run')]))){
+if (!is.null(params$functional_method)){
     cat("\n### Gene set analysis\n")
+    cat("\n#### ", toupper(params$functional_method) ," {.tabset}\n")
+
+    if (params$functional_method == 'gsea') {
+        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) ))
+        }
+    }
 
-    for (gene_set_method in possible_gene_set_methods){
-        if (unlist(params[paste0(gene_set_method, '_run')])){
-        cat("\n#### ", toupper(gene_set_method) ," {.tabset}\n")
-        if (gene_set_method == 'gsea') {
-            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$report_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$report_round_digits)
-                print( htmltools::tagList(datatable(ref_gsea_results, caption = paste0("\nReference (", contrasts$reference[i], ")\n"), rownames = FALSE) ))
-            }
-            }
-
-        } else if (gene_set_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"))
-
-            # Make sure to grab only non-empty files
-            for (name in differential_names) {
-                cat(paste0("\n##### ", name, "\n"))
-
-                table <- paste0(name, ".gprofiler2.all_enriched_pathways.tsv")
-                table_path <- file.path(params$input_dir, table)
-                if (!file.exists(table_path) || file.size(table_path) == 0){
-                    cat(paste0("No ", ifelse(params$gprofiler2_significant, "significantly", ""), " enriched pathways were found for this contrast."))
-                } else {
-                    all_enriched <- read.table(table_path, header=T, sep="\t", quote="\"")
-                    all_enriched <- data.frame("Pathway name" = all_enriched$term_name, "Pathway code" = all_enriched$term_id,
-                    "Differential features" = all_enriched$intersection_size, "Pathway size" = all_enriched$term_size,
-                    "Differential fraction" = (all_enriched$intersection_size/all_enriched$term_size),
-                    "Adjusted p value" = all_enriched$p_value, check.names = FALSE)
-                    all_enriched <- round_dataframe_columns(all_enriched, digits=params$report_round_digits)
-                    print(htmltools::tagList(datatable(all_enriched, caption = paste('Enriched pathways in', name, " (check", table, "for more detail)"), rownames = FALSE)))
-                }
-                cat("\n")
+    } 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"))
+
+        # Make sure to grab only non-empty files
+        for (name in differential_names) {
+            cat(paste0("\n##### ", name, "\n"))
+
+            table <- paste0(name, ".gprofiler2.all_enriched_pathways.tsv")
+            table_path <- file.path(params$input_dir, table)
+            if (!file.exists(table_path) || file.size(table_path) == 0){
+                cat(paste0("No ", ifelse(params$gprofiler2_significant, "significantly", ""), " enriched pathways were found for this contrast."))
+            } else {
+                all_enriched <- read.table(table_path, header=T, sep="\t", quote="\"")
+                all_enriched <- data.frame("Pathway name" = all_enriched$term_name, "Pathway code" = all_enriched$term_id,
+                "Differential features" = all_enriched$intersection_size, "Pathway size" = all_enriched$term_size,
+                "Differential fraction" = (all_enriched$intersection_size/all_enriched$term_size),
+                "Adjusted p value" = all_enriched$p_value, check.names = FALSE)
+                all_enriched <- round_dataframe_columns(all_enriched, digits=params$round_digits)
+                print(htmltools::tagList(datatable(all_enriched, caption = paste('Enriched pathways in', name, " (check", table, "for more detail)"), rownames = FALSE)))
             }
-        }
+            cat("\n")
         }
     }
 }
@@ -1066,7 +1059,7 @@ make_params_table('exploratory analysis', 'exploratory_', remove_pattern = TRUE)
 
 
 ```{r, echo=FALSE, results='asis'}
-if (params$study_type == 'rnaseq'){
+if (params$differential_method == 'deseq2'){
     make_params_table('DESeq2', 'deseq2_', remove_pattern = TRUE)
 }
 make_params_table('downstream differential analysis', 'differential_', remove_pattern = TRUE)
@@ -1075,18 +1068,10 @@ make_params_table('downstream differential analysis', 'differential_', remove_pa
 <!-- If any gene set methods have been activated show their params -->
 
 ```{r, echo=FALSE, results='asis'}
-possible_gene_set_methods <- c('gsea', 'gprofiler2')
-
-if (any(unlist(params[paste0(possible_gene_set_methods, '_run')]))){
+if (!is.null(params$functional_method)){
     cat("\n### Gene set analysis\n")
-
-    for (gene_set_method in possible_gene_set_methods){
-        if (unlist(params[paste0(gene_set_method, '_run')])){
-        cat("\n#### ", toupper(gene_set_method) ,"  {.tabset}\n")
-        make_params_table(toupper(gene_set_method), paste0(gene_set_method, '_'), remove_pattern = TRUE)
-        }
-    }
-
+    cat("\n#### ", toupper(params$functional_method) ,"  {.tabset}\n")
+    make_params_table(toupper(params$functional_method), paste0(params$functional_method, '_'), remove_pattern = TRUE)
 }
 ```
 

assets/paramsheet.csv

@@ -0,0 +1,15 @@
+paramset_name,study_type,differential_method,limma_use_voom,functional_method
+deseq2_rnaseq,rnaseq,deseq2,,
+deseq2_rnaseq_gsea,rnaseq,deseq2,,gsea
+deseq2_rnaseq_gprofiler2,rnaseq,deseq2,,gprofiler2
+limma_rnaseq,rnaseq,limma,true,
+limma_rnaseq_gsea,rnaseq,limma,true,gsea
+limma_rnaseq_gprofiler2,rnaseq,limma,true,gprofiler2
+dream_rnaseq,rnaseq,dream,,
+dream_rnaseq_gsea,rnaseq,dream,,gsea
+dream_rnaseq_gprofiler2,rnaseq,dream,,gprofiler2
+limma_affy,affy_array,limma,,
+limma_affy_gsea,affy_array,limma,,gsea
+limma_affy_gprofiler2,affy_array,limma,,gprofiler2
+limma_soft,geo_soft_file,limma,,
+limma_maxquant,maxquant,limma,,

conf/affy.config

@@ -30,6 +30,7 @@ params {
     exploratory_log2_assays = null
 
     // Differential options
+    differential_method              = "limma"
     differential_file_suffix         = ".limma.results.tsv"
     differential_fc_column           = "logFC"
     differential_pval_column         = "P.Value"

conf/maxquant.config

@@ -31,6 +31,7 @@ params {
     exploratory_log2_assays = null
 
     // Differential options
+    differential_method              = "limma"
     differential_file_suffix         = ".limma.results.tsv"
     differential_fc_column           = "logFC"
     differential_pval_column         = "P.Value"