Merge pull request #320 from suzannejin/dev-ratio-propd-loop

[dev-ratio] update propr container and propd.R

assets/tools_samplesheet.csv

@@ -1,8 +1,10 @@
 pathway_name,diff_method,args_diff,cor_method,args_cor,enr_method,args_enr
 propd,propd,,,,,
 propd_fdr,propd,--permutation 100,,,,
+propd_grea,propd,,,,grea,
+propd_ora,propd,,,,gprofiler2,
+deseq2_gsea,deseq2,,,,gsea,
+deseq2_ora,deseq2,,,,gprofiler2,
 pcorbshrink,,,propr,--metric pcor.bshrink,,
 propr,,,propr,--metric rho,,
 cor,,,propr,--metric cor,,
-propd_grea,propd,,,,grea,
-deseq2,deseq2,,,,gsea,

conf/test_experimental.config

@@ -29,7 +29,11 @@ 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 =  'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/mus_musculus/rnaseq_expression/SRP254919.contrasts.csv'
-    pathway = "propd,propd_fdr,pcorbshrink,propr,cor,propd_grea"
+    pathway = "propd,propd_fdr,propd_grea,propr"
+
+    // some pathways are not included in the test, because they don't make sense to test for this dataset
+    // for example propd_fdr don't give significant results on this dataset, and the permutation will be scaped
+    // and pcor.bshrink is taking too long for this test (15min), and it should be tested on matrices with n > p, or n close to p. (eg. specific pathways, or on DE genes only, etc)
 
     //Features
     features_metadata_cols      = 'gene_id,gene_name'

modules/local/propr/grea/main.nf

@@ -4,8 +4,8 @@ process PROPR_GREA {
 
     // conda "${moduleDir}/environment.yml"
     container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
-        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:75e6dfd62a2313a9':
-        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:fbd569ab00953cb0' }"
+        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:209490acb0e524e3' :
+        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:17abd3f137436739' }"
 
     input:
     tuple val(meta), path(adj)

modules/local/propr/propd/main.nf

@@ -4,8 +4,8 @@ process PROPR_PROPD {
 
     // conda "${moduleDir}/environment.yml"
     container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
-        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:75e6dfd62a2313a9':
-        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:fbd569ab00953cb0' }"
+        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:209490acb0e524e3' :
+        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:17abd3f137436739' }"
 
     input:
     tuple val(meta), path(count), path(samplesheet), val(contrast_variable), val(reference), val(target)

modules/local/propr/propd/templates/propd.R

@@ -465,21 +465,50 @@ if (opt\$permutation == 0) {
 
     warning('Permutation tests are used to compute FDR values.')
 
-    # update FDR values using permutation tests
+    # calculate FDR values using permutation tests
+    # this part will call the updateCutoffs function iteratively
+    # as far as it does not find a meaningful theta value
+    # and does not reach the maximum number of iterations
+
+    theta_cutoff <- FALSE
+    max_cutoff <- 1
+    ntry <- 0
+    while (!theta_cutoff & max_cutoff > 0 & ntry < 10) {
+        ntry <- ntry + 1
+
+        # get theta cutoffs to test the FDR
+
+        if (ntry > 1) {
+            part <- pd@fdr[which(pd@fdr\$truecounts > 0),]
+            if (nrow(part) > 1) {
+                max_cutoff <- min(part\$cutoff)
+            } else {
+                break
+            }
+        }
 
-    pd <- updateCutoffs(
-        pd,
-        number_of_cutoffs = 100,
-        ncores = opt\$ncores
-    )
+        cutoffs <- as.numeric(quantile(
+            pd@results[pd@results\$theta < max_cutoff, 'theta'],
+            seq(0, 1, length.out = opt\$number_of_cutoffs)
+        ))
 
-    # get theta cutoff
+        # update FDR values
+
+        pd <- updateCutoffs(
+            pd,
+            custom_cutoffs = cutoffs,
+            ncores = opt\$ncores
+        )
+
+        # check if any theta value has FDR below desired threshold
+
+        theta_cutoff <- getCutoffFDR(
+            pd,
+            fdr=opt\$fdr,
+            window_size=1
+        )
+    }
 
-    theta_cutoff <- getCutoffFDR(
-        pd,
-        fdr=opt\$fdr,
-        window_size=1
-    )
     if (theta_cutoff) {
 
         warning('Significant theta value found: ', theta_cutoff)

modules/local/propr/propr/main.nf

@@ -4,8 +4,8 @@ process PROPR_PROPR {
 
     // conda "${moduleDir}/environment.yml"
     container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
-        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:75e6dfd62a2313a9':
-        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:fbd569ab00953cb0' }"
+        'oras://community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:209490acb0e524e3' :
+        'community.wave.seqera.io/library/bioconductor-limma_r-ggplot2_r-propr:17abd3f137436739' }"
 
     input:
     tuple val(meta), path(count)

subworkflows/local/differential/main.nf

@@ -39,6 +39,9 @@ workflow DIFFERENTIAL {
     // Perform differential analysis with propd
     // ----------------------------------------------------
 
+    // TODO propd currently don't support blocking, so we should not run propd with same contrast_variable, reference and target,
+    // but different blocking variable, since it will simply run the same analysis again.
+
     ch_counts
         .join(ch_samplesheet)
         .combine(ch_contrasts)