Merge branch 'master' into new_tag_strategy

modules/nf-core/dragmap/hashtable/tests/main.nf.test

@@ -16,7 +16,7 @@ nextflow_process {
                 """
                 input[0] = [
                     [id:'test'],
-                    file(params.test_data['sarscov2']['genome']['genome_fasta'], checkIfExists: true)
+                    file(params.modules_testdata_base_path + 'genomics/sarscov2/genome/genome.fasta', checkIfExists: true)
                 ]
                 """
             }
@@ -28,7 +28,7 @@ nextflow_process {
                 { assert snapshot(
                         file(process.out.hashmap[0][1]).name,
                         file(process.out.versions[0]).name
-                    ).match() 
+                    ).match()
                 }
             )
         }
@@ -44,7 +44,7 @@ nextflow_process {
                 """
                 input[0] = [
                     [id:'test'],
-                    file(params.test_data['sarscov2']['genome']['genome_fasta'], checkIfExists: true)
+                    file(params.modules_testdata_base_path + 'genomics/sarscov2/genome/genome.fasta', checkIfExists: true)
                 ]
                 """
             }
@@ -56,7 +56,7 @@ nextflow_process {
                 { assert snapshot(
                         file(process.out.hashmap[0][1]).name,
                         file(process.out.versions[0]).name
-                    ).match() 
+                    ).match()
                 }
             )
         }

modules/nf-core/propr/grea/meta.yml

@@ -0,0 +1,53 @@
+---
+# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json
+name: "propr_grea"
+description: Perform Gene Ratio Enrichment Analysis
+keywords:
+  - logratio
+  - differential
+  - propr
+  - grea
+  - enrichment
+  - expression
+tools:
+  - "grea":
+      description: "Gene Ratio Enrichment Analysis"
+      homepage: "https://github.com/tpq/propr"
+      documentation: "https://rdrr.io/cran/propr/man/propr.html"
+      tool_dev_url: "https://github.com/tpq/propr"
+      doi: "10.2202/1544-6115.1175"
+      licence: ["GPL-2"]
+input:
+  - meta:
+      type: map
+      description: |
+        Groovy Map containing sample information.
+        This can be used at the workflow level to pass optional parameters to the module.
+        [id: 'test', ...]
+  - A:
+      type: file
+      description: adjacency matrix for gene ratio proportionality/differential proportionality
+      pattern: "*.{csv,tsv}"
+  - K:
+      type: file
+      description: relational database containing genes and GO terms (generated by mygene module)
+      pattern: "*.{gmt}"
+output:
+  - meta:
+      type: map
+      description: |
+        Groovy Map containing sample information.
+        This can be used at the workflow level to pass optional parameters to the module.
+        [id: 'test', ...]
+  - enriched:
+      type: file
+      description: File containing GO terms and their enrichment values
+      pattern: "*.{csv}"
+  - versions:
+      type: file
+      description: File containing software versions
+      pattern: "versions.yml"
+authors:
+  - "@caraiz2001"
+maintainers:
+  - "@caraiz2001"

modules/nf-core/propr/propr/main.nf

@@ -14,6 +14,8 @@ process PROPR_PROPR {
     tuple val(meta), path("*.propr.rds"), emit: propr
     tuple val(meta), path("*.propr.tsv"), emit: matrix
     tuple val(meta), path("*.fdr.tsv"),   emit: fdr         , optional:true
+    tuple val(meta), path("*.adj.csv"),   emit: adj         , optional:true
+    path "*.warnings.log",                emit: warnings
     path "*.R_sessionInfo.log",           emit: session_info
     path "versions.yml",                  emit: versions
 

modules/nf-core/propr/propr/meta.yml

@@ -59,6 +59,10 @@ output:
       type: file
       description: (optional) propr fdr table
       pattern: "*.fdr.tsv"
+  - adj:
+      type: file
+      description: (optional) propr adjacency table
+      pattern: "*.adj.csv"
   - session_info:
       type: file
       description: dump of R SessionInfo

modules/nf-core/propr/propr/templates/propr.R

@@ -96,6 +96,79 @@ set_reference <- function(ivar, mat){
     return(ivar)
 }
 
+#' Set the appropiate range for the sequence of cutoffs used in updateCutoffs.
+#' Adjusts the interval to the different metrics.
+#'
+#' @param object propr object. Output from propr function.
+#'
+#' @return sequence of cutoff values.
+seqCutoff <- function(object){
+    matrix <- getMatrix(object)
+    matrix[matrix <= 0] <- NA
+    diag(matrix) <- NA
+    min_cutoff <- round(min(matrix, na.rm = TRUE),3)
+    max_cutoff <- round(max(matrix, na.rm = TRUE),3)
+    step_cutoff <- (max_cutoff - min_cutoff)/ 20
+    seq_cutoff <- seq(min_cutoff, max_cutoff, step_cutoff)
+    return(seq_cutoff)
+}
+
+#' Extract the proportionality cutoff for a specified FDR value.
+#' Gene pairs with a proportionality value higher than the extracted cutoff will be considered significantly proportional.
+#'
+#' @param object propr object. Output from propr function. updateCutoffs function should be applied to the object previous to valCutoff.
+#' @param fdrVal FDR value to extract the cutoff for. Per default 0.05
+#' @param metric Metric used to calculate the proportionality values. Options are 'cor', 'rho', 'phi', 'phs', 'vlr', 'pcor', 'pcor.shrink', 'pcor.bshrink'
+#'
+#' @return cutoff value. Proportionality values higher than this cutoff are considered significant.
+valCutoff  <- function(object, metric, fdrVal = 0.05){
+    fdr_df <- object@fdr
+    print(fdr_df)
+    metric_up <- c("rho", "cor", "pcor", "pcor.shrink", "pcor.bshrink")
+    if (prod(dim(fdr_df) == 0)){
+        warning("Please run updateCutoff on propr first")
+    }else{
+        fdr_vals <- fdr_df\$FDR
+        if (any(is.na(fdr_vals))){
+            stop("FDR not defined. This metric is not appropiate for the given dataset")
+        }
+        threshold <- any(fdr_vals <= fdrVal)
+        if (metric %in% metric_up){
+            if (threshold){
+                fdr_threshold <- fdr_vals[which.max(fdr_vals <= fdrVal)]
+            }else{
+                warning("FDR is higher than the specified threshold for all proportionality values. Using the lowest fdr instead")
+                fdr_threshold <- fdr_vals[length(fdr_vals)]
+            }
+        }else{
+            if (threshold){
+                fdr_threshold <- fdr_vals[which.min(fdr_vals <= fdrVal) - 1]
+            }else{
+                warning("FDR is higher than the specified threshold for all proportionality values. Using the lowest fdr instead")
+                fdr_threshold <- fdr_vals[1]
+            }
+        }
+    cutoff <- fdr_df\$cutoff[fdr_df\$FDR == fdr_threshold]
+    }
+    return(cutoff)
+}
+
+
+#' Convert a proportionality matrix to an adjacency matrix based on a threshold.
+#'
+#' @param matrix proportionality matrix. Can be extracted from propr object with getMatrix().
+#' @param cutoff Significant proportionality value extracted from valCutoff function.
+#'
+#' @return Adjacency matrix. Gene pairs with a proportionality value higher than the threshold will have 1, otherwise 0.
+convert_to_adjacency <- function(matrix, cutoff, metric) {
+    if (metric == 'cor' || metric == 'rho' || metric == 'pcor' || metric == 'pcor.shrink' || metric == 'pcor.bshrink'){
+        adjacency <- ifelse(matrix > cutoff, 1, 0)
+    } else {
+        adjacency <- ifelse(matrix < cutoff, 1, 0)
+    }
+    return(adjacency)
+}
+
 ################################################
 ################################################
 ## Parse arguments                            ##
@@ -115,7 +188,9 @@ opt <- list(
     cutoff_interval  = NA,
     ncores           = as.integer('$task.cpus'),
     features_id_col  = 'gene_id',
-    fixseed          = FALSE
+    fixseed          = FALSE,
+    adjacency        = FALSE,
+    fdrVal           = 0.05
 )
 opt_types <- list(
     count            = 'character',
@@ -130,11 +205,14 @@ opt_types <- list(
     cutoff_interval  = 'numeric',
     ncores           = 'numeric',
     features_id_col  = 'character',
-    fixseed          = 'logical'
+    fixseed          = 'logical',
+    adjacency        = 'logical',
+    fdrVal           = 'numeric'
 )
 
 # Apply parameter overrides
 args_opt <- parse_args('$task.ext.args')
+
 for ( ao in names(args_opt)){
     if (! ao %in% names(opt)){
         stop(paste("Invalid option:", ao))
@@ -235,15 +313,36 @@ pro <- propr(
 )
 
 # update FDR by permutation, if required
+
 if (opt\$permutation > 0) {
     cutoff <- seq(
         opt\$cutoff_min,
         opt\$cutoff_max,
         opt\$cutoff_interval
         )
+    if (is.na(opt\$cutoff_min) ||  is.na(opt\$cutoff_max) || is.na(opt\$cutoff_interval)) {
+        warning("cutoff values were not provided. Using the default cutoff values.")
+        cutoff <- seqCutoff(pro)
+    }
+    m <- getMatrix(pro)
+    diag(m) <- NA
+    print((opt\$cutoff_max - opt\$cutoff_min)/2 + opt\$cutoff_min)
+    print(max(m, na.rm = TRUE))
+    if ((opt\$cutoff_max - opt\$cutoff_min)/2 + opt\$cutoff_min > max(m, na.rm = TRUE)) {
+        warning("The provided cutoff values are out of range. Using the default cutoff values.")
+        cutoff <- seqCutoff(pro)
+    }
     pro <- updateCutoffs(pro, cutoff=cutoff, ncores=opt\$ncores)
 }
 
+# calculate cutoff and adjacency matrix, if required
+
+if (opt\$adjacency == TRUE) {
+    cutoff <- valCutoff(pro, opt\$metric, opt\$fdrVal)
+    matrix <- getMatrix(pro)
+    adj <- convert_to_adjacency(matrix, cutoff, opt\$metric)
+}
+
 ################################################
 ################################################
 ## Generate outputs                           ##
@@ -275,6 +374,27 @@ if (opt\$permutation > 0) {
     )
 }
 
+if (opt\$adjacency == TRUE) {
+    write.table(
+        adj,
+        file      = paste0(opt\$prefix, '.adj.csv'),
+        col.names = TRUE,
+        row.names = TRUE,
+        sep       = ',',
+        quote     = FALSE
+    )
+}
+
+################################################
+################################################
+## WARNINGS                                   ##
+################################################
+################################################
+
+sink(paste0(opt\$prefix, ".warnings.log"))
+print(warnings())
+sink()
+
 ################################################
 ################################################
 ## R SESSION INFO                             ##
@@ -291,13 +411,11 @@ sink()
 ################################################
 ################################################
 
-r.version <- strsplit(version[['version.string']], ' ')[[1]][3]
 propr.version <- as.character(packageVersion('propr'))
 
 writeLines(
     c(
         '"${task.process}":',
-        paste('    r-base:', r.version),
         paste('    r-propr:', propr.version)
     ),
 'versions.yml')

modules/nf-core/propr/propr/tests/adjacency.config

@@ -0,0 +1,3 @@
+process {
+    ext.args = {"--metric rho --permutation 10 --cutoff_min 0.05 --cutoff_max 0.95 --cutoff_interval 0.1 --fixseed true --adjacency true"}
+}

modules/nf-core/propr/propr/tests/adjacency_pcor.config

@@ -0,0 +1,3 @@
+process {
+    ext.args = {"--metric pcor.bshrink --permutation 10 --cutoff_min 0.05 --cutoff_max 0.95 --cutoff_interval 0.1 --fixseed true --adjacency true"}
+}

modules/nf-core/propr/propr/tests/adjacency_phs.config

@@ -0,0 +1,3 @@
+process {
+    ext.args = {"--metric phs --permutation 10 --cutoff_min 0.05 --cutoff_max 0.95 --cutoff_interval 0.1 --fixseed true --adjacency true"}
+}

modules/nf-core/propr/propr/tests/adjacency_rho.config

@@ -0,0 +1,3 @@
+process {
+    ext.args = {"--metric rho --permutation 10 --cutoff_min 0.05 --cutoff_max 0.95 --cutoff_interval 0.1 --fixseed true --adjacency true"}
+}