Merge pull request #182 from jr-leary7/dev

added parallel processing to geneProgramDrivers

Author: jr-leary7

DESCRIPTION

@@ -9,7 +9,7 @@ Description: This package uses truncated power basis spline models to build flex
              Downstream analysis functionalities include model comparison, dynamic gene clustering, smoothed counts generation, gene set enrichment testing, & visualization. 
 License: MIT + file LICENSE
 Encoding: UTF-8
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Depends: 
   glm2, 
   magrittr, 

R/geneProgramDrivers.R

@@ -4,6 +4,11 @@
 #' @author Jack Leary
 #' @importFrom Matrix Matrix
 #' @importFrom purrr map reduce
+#' @importFrom foreach foreach registerDoSEQ
+#' @importFrom parallel makeCluster clusterEvalQ stopCluster
+#' @importFrom withr with_output_sink
+#' @importFrom utils txtProgressBar setTxtProgressBar
+#' @importFrom doSNOW registerDoSNOW
 #' @importFrom stats cor.test p.adjust
 #' @importFrom dplyr arrange desc mutate filter
 #' @description This function computes the correlation
@@ -13,6 +18,8 @@
 #' @param cor.method (Optional) The correlation method to be used. Defaults to "spearman".
 #' @param fdr.cutoff (Optional) The FDR threshold for determining statistical significance. Defaults to 0.01.
 #' @param p.adj.method (Optional) The method used to adjust \emph{p}-values for multiple hypothesis testing. Defaults to "holm".
+#' @param n.cores (Optional) The number of cores used when iterating over genes to perform testing. Defaults to 2.
+#' @param verbose (Optional) Should a progress bar be printed to the console during processing? Defaults to TRUE.
 #' @return Either a \code{Seurat} or \code{SingleCellExperiment} object if \code{expr.mat} is in either form, or a data.frame containing per-cell program scores if \code{expr.mat} is a matrix.
 #' @seealso \code{\link{geneProgramScoring}}
 #' @seealso \code{\link[stats]{cor.test}}
@@ -39,7 +46,9 @@ geneProgramDrivers <- function(expr.mat = NULL,
                                gene.program = NULL,
                                cor.method = "spearman",
                                fdr.cutoff = 0.01,
-                               p.adj.method = "holm") {
+                               p.adj.method = "holm",
+                               n.cores = 2L,
+                               verbose = TRUE) {
   # check inputs
   if (is.null(expr.mat) || is.null(genes) || is.null(gene.program)) { stop("Arguments to geneProgramDrivers() are missing.") }
   cor.method <- tolower(cor.method)
@@ -54,16 +63,50 @@ geneProgramDrivers <- function(expr.mat = NULL,
   } else if (inherits(expr.mat, "dgCMatrix")) {
     counts_matrix <- Matrix::Matrix(expr.mat, sparse = FALSE)
   }
+  # set up parallel processing
+  if (n.cores > 1L) {
+    cl <- parallel::makeCluster(n.cores)
+    doSNOW::registerDoSNOW(cl)
+  } else {
+    cl <- foreach::registerDoSEQ()
+  }
+  # set up progress bar
+  if (verbose) {
+    withr::with_output_sink(tempfile(), {
+      pb <- utils::txtProgressBar(0, length(genes), style = 3)
+    })
+    progress_fun <- function(n) utils::setTxtProgressBar(pb, n)
+    snow_opts <- list(progress = progress_fun)
+  } else {
+    snow_opts <- list()
+  }
   # iteratively compute correlations & p-values
-  cor_tests <- purrr::map(genes, \(g) {
-    cor_res <- stats::cor.test(counts_matrix[g, ],
+  cor_tests <- foreach::foreach(g = seq(genes),
+                                .combine = "list",
+                                .multicombine = ifelse(length(genes) > 1, TRUE, FALSE),
+                                .maxcombine = ifelse(length(genes) > 1, length(genes), 2),
+                                .packages = c("stats"),
+                                .errorhandling = "pass",
+                                .inorder = TRUE,
+                                .verbose = FALSE,
+                                .options.snow = snow_opts) %dopar% {
+    cor_res <- stats::cor.test(counts_matrix[genes[g], ],
                                gene.program,
                                method = cor.method,
                                exact = FALSE)
-    cor_df <- data.frame(gene = g,
+    cor_df <- data.frame(gene = genes[g],
                          corr = unname(cor_res$estimate),
                          pvalue = cor_res$p.value)
-    return(cor_df)
+    cor_df
+  }
+  # end parallelization & clean up each worker node
+  withr::with_output_sink(tempfile(), {
+    if (n.cores > 1L) {
+      parallel::clusterEvalQ(cl, expr = {
+        gc(verbose = FALSE, full = TRUE)
+      })
+      parallel::stopCluster(cl)
+    }
   })
   cor_tests <- purrr::reduce(cor_tests, rbind)
   rownames(cor_tests) <- genes

R/geneProgramScoring.R

@@ -8,6 +8,7 @@
 #' @param genes A character vector of gene IDs. Defaults to NULL.
 #' @param gene.clusters A factor containing the cluster assignment of each gene in \code{genes}. Defaults to NULL.
 #' @param program.labels (Optional) A character vector specifying a label for each gene cluster. Defaults to NULL.
+#' @param minmax.norm (Optional) Should each program's score be min-max normalized to be on [0, 1]? Defaults to FALSE.
 #' @param n.cores (Optional) The number of cores used under the hood in \code{\link[UCell]{ScoreSignatures_UCell}}. Defaults to 2.
 #' @return Either a \code{Seurat} or \code{SingleCellExperiment} object if \code{expr.mat} is in either form, or a data.frame containing per-cell program scores if \code{expr.mat} is a matrix.
 #' @seealso \code{\link[UCell]{ScoreSignatures_UCell}}
@@ -30,6 +31,7 @@ geneProgramScoring <- function(expr.mat = NULL,
                                genes = NULL,
                                gene.clusters = NULL,
                                program.labels = NULL,
+                               minmax.norm = FALSE,
                                n.cores = 2L) {
   # check inputs
   if (is.null(expr.mat) || is.null(genes) || is.null(gene.clusters)) { stop("Arguments to geneProgramScoring() are missing.") }
@@ -42,6 +44,9 @@ geneProgramScoring <- function(expr.mat = NULL,
     program.labels <- paste0("cluster_", cluster.labels)
   } else {
     program.labels <- gsub(" ", "_", program.labels)
+    if (length(program.labels) != length(levels(gene.clusters))) {
+      stop("Each cluster must have a label.")
+    }
   }
   # set up query matrix
   if (inherits(expr.mat, "SingleCellExperiment")) {
@@ -60,6 +65,15 @@ geneProgramScoring <- function(expr.mat = NULL,
   program_scores <- UCell::ScoreSignatures_UCell(counts_matrix,
                                                  features = program_list,
                                                  ncores = n.cores)
+  # min-max normalize if desired
+  if (minmax.norm) {
+    program_scores <- purrr::map(seq(ncol(program_scores)), \(i) {
+      scores <- program_scores[, i]
+      normed_scores <- (scores - min(scores)) / (max(scores) - min(scores))
+      return(normed_scores)
+    })
+    program_scores <- purrr::reduce(program_scores, cbind)
+  }
   # reformat program scores depending on input format
   if (inherits(expr.mat, "matrix") || inherits(expr.mat, "array") || inherits(expr.mat, "dgCMatrix")) {
     colnames(program_scores) <- program.labels