Merge pull request #139 from jr-leary7/dev

updated tests with better error handling

.github/workflows/R-CMD-check.yaml

@@ -2,7 +2,7 @@
 
 on:
   push:
-    branches: [main, master]
+    branches: main
 
 name: R-CMD-check
 

.github/workflows/render-README.yaml

@@ -3,7 +3,7 @@
 
 on:
   push:
-    branches: [main, master]
+    branches: main
     paths: README.Rmd
 
 name: render-README
@@ -16,15 +16,17 @@ jobs:
     steps:
       - uses: actions/checkout@v3
       - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
       - uses: r-lib/actions/setup-r-dependencies@v2
       - uses: r-lib/actions/setup-pandoc@v2
 
       - name: install CRAN packages
         run: Rscript -e 'install.packages(c("rmarkdown","ggplot2", "dplyr", "purrr", "remotes"))'
       - name: install BioConductor packages
-        run: Rscript -e 'install.packages("BiocManager"); BiocManager::install(c("SingleCellExperiment", "scater", "scran", "scaffold"))'
+        run: Rscript -e 'install.packages("BiocManager"); BiocManager::install(c("SingleCellExperiment", "scater", "scran")'
       - name: install GitHub packages
-        run: Rscript -e 'remotes::install_github("jr-leary7/scLANE")'
+        run: Rscript -e 'remotes::install_github("jr-leary7/scLANE"); remotes::install_github("rhondabacher/scaffold")'
 
       - name: render README
         run: Rscript -e 'rmarkdown::render("README.Rmd", output_format = "md_document")'
@@ -33,4 +35,4 @@ jobs:
         run: |
           git add README.md man/figures/README-*
           git commit -m "Recompiled README.Rmd" || echo "No changes to commit"
-          git push origin main || echo "No changes to commit"
+          git push origin || echo "No changes to commit"

.github/workflows/test-coverage.yaml

@@ -2,7 +2,7 @@
 
 on:
   push:
-    branches: [main, master]
+    branches: main
 
 name: test-coverage
 

R/clusterGenes.R

@@ -42,7 +42,7 @@ clusterGenes <- function(test.dyn.res = NULL,
                          size.factor.offset = NULL,
                          clust.algo = "leiden",
                          use.pca = FALSE,
-                         n.PC = 15,
+                         n.PC = 15L,
                          lineages = NULL) {
   # check inputs
   if (is.null(test.dyn.res) || is.null(pt)) { stop("test.dyn.res & pt must be supplied to clusterGenes().") }
@@ -59,6 +59,7 @@ clusterGenes <- function(test.dyn.res = NULL,
     fitted_vals_mat <- purrr::map(test.dyn.res, \(x) x[[lineage_name]]$MARGE_Preds) %>%
                        stats::setNames(names(test.dyn.res)) %>%
                        purrr::discard(rlang::is_na) %>%
+                       purrr::discard(rlang::is_null) %>%
                        purrr::discard(\(p) rlang::inherits_only(p, "try-error")) %>%
                        purrr::map2(.y = names(.), function(x, y) {
                          if (is.null(size.factor.offset)) {

R/createSlopeTestData.R

@@ -48,10 +48,16 @@ createSlopeTestData <- function(marge.model = NULL,
       brkpts <- purrr::map_dbl(rounded_brkpts, \(x) pt[, 1][which.min(abs(pt[, 1] - x))])
       # extract p-values for coefficients other than intercept
       if (is.gee) {
+        est_betas <- summary(marge.model$final_mod)$beta[-1]
+        test_stats <- summary(marge.model$final_mod)$wald.test[-1]
         p_vals <- summary(marge.model$final_mod)$p[-1]
       } else if (is.glmm) {
+        est_betas <- unname(summary(marge.model$final_mod)$coefficients$cond[, "Estimate"][-1])
+        test_stats <- unname(summary(marge.model$final_mod)$coefficients$cond[, "z value"][-1])
         p_vals <- unname(summary(marge.model$final_mod)$coefficients$cond[, "Pr(>|z|)"][-1])
       } else {
+        est_betas <- unname(summary(marge.model$final_mod)$coefficients[, "Estimate"][-1])
+        test_stats <- unname(summary(marge.model$final_mod)$coefficients[, "z value"][-1])
         p_vals <- unname(summary(marge.model$final_mod)$coefficients[, "Pr(>|z|)"][-1])
       }
       mod_notes <- rep(NA_character_, length(p_vals))
@@ -61,6 +67,8 @@ createSlopeTestData <- function(marge.model = NULL,
   slope_df <- data.frame(Breakpoint = brkpts,
                          Rounded_Breakpoint = rounded_brkpts,
                          Direction = brkpt_dirs,
+                         Beta = est_betas,
+                         Test_Stat = test_stats,
                          P_Val = p_vals,
                          Notes = mod_notes)
   return(slope_df)

R/getFittedValues.R

@@ -83,7 +83,7 @@ getFittedValues <- function(test.dyn.res = NULL,
                              data.frame(scLANE_pred_link = test.dyn.res[[g]][[paste0("Lineage_", y)]]$MARGE_Preds$marge_link_fit,
                                         scLANE_se_link = test.dyn.res[[g]][[paste0("Lineage_", y)]]$MARGE_Preds$marge_link_se)
                            }, silent = TRUE)
-                           if (inherits(pred_df, "try-error")) {
+                           if (inherits(pred_df, "try-error") || is.null(pred_df) || all(is.na(pred_df))) {
                              gene_df <- dplyr::mutate(gene_df,
                                                       scLANE_pred_link = NA_real_,
                                                       scLANE_se_link = NA_real_,

R/getKnotDist.R

@@ -4,6 +4,7 @@
 #' @author Jack Leary
 #' @import magrittr
 #' @importFrom purrr imap reduce
+#' @importFrom dplyr filter
 #' @description Pulls knot locations for dynamic genes across each lineage, allowing comparisons of where transcriptional switches occur between lineages.
 #' @param test.dyn.res The output from \code{\link{testDynamic}}. Defaults to NULL.
 #' @param dyn.genes The set of genes to pull knots for. If unspecified, pulls knots for all modeled genes. Defaults to NULL.
@@ -22,13 +23,19 @@ getKnotDist <- function(test.dyn.res = NULL, dyn.genes = NULL) {
   }
   # pull knots per-lineage
   knot_df <- purrr::imap(test.dyn.res[dyn.genes], \(x, y) {
-    purrr::imap(x, \(z, w) {
-      data.frame(gene = y,
-                 lineage = w,
-                 knot = z$MARGE_Slope_Data$Breakpoint)
-    }) %>%
-      purrr::reduce(rbind)
-  }) %>%
-    purrr::reduce(rbind)
+              purrr::imap(x, \(z, w) {
+                if (is.null(z$MARGE_Slope_Data$Breakpoint)) {
+                  knots <- NA_real_
+                } else {
+                  knots <- z$MARGE_Slope_Data$Breakpoint
+                }
+                data.frame(gene = y,
+                           lineage = w,
+                           knot = knots)
+              }) %>%
+                 purrr::reduce(rbind)
+             }) %>%
+               purrr::reduce(rbind) %>%
+               dplyr::filter(!is.na(knot))
   return(knot_df)
 }

R/smoothedCountsMatrix.R

@@ -54,6 +54,7 @@ smoothedCountsMatrix <- function(test.dyn.res = NULL,
     fitted_vals_list <- furrr::future_map(test.dyn.res, \(y) y[[lineage_name]]$MARGE_Preds) %>%
                         stats::setNames(names(test.dyn.res)) %>%
                         purrr::discard(rlang::is_na) %>%
+                        purrr::discard(rlang::is_null) %>%
                         purrr::discard(\(p) rlang::inherits_only(p, "try-error"))
     fitted_vals_mat <- purrr::map(fitted_vals_list, \(z) {
                          if (is.null(size.factor.offset)) {

R/testDynamic.R

@@ -140,10 +140,11 @@ testDynamic <- function(expr.mat = NULL,
   if (parallel.exec) {
     cl <- parallel::makeCluster(n.cores)
     doParallel::registerDoParallel(cl)
+    parallel::clusterSetRNGStream(cl, iseed = random.seed)
   } else {
     cl <- foreach::registerDoSEQ()
+    set.seed(random.seed)
   }
-  parallel::clusterSetRNGStream(cl, iseed = random.seed)
 
   # convert dense counts matrix to file-backed matrix
   expr.mat <- bigstatsr::as_FBM(expr.mat,
@@ -342,10 +343,12 @@ testDynamic <- function(expr.mat = NULL,
 
   # end parallelization & clean up each worker node
   withr::with_output_sink(tempfile(), {
-    parallel::clusterEvalQ(cl, expr = {
-      gc(verbose = FALSE, full = TRUE)
-    })
-    parallel::stopCluster(cl)
+    if (parallel.exec) {
+      parallel::clusterEvalQ(cl, expr = {
+        gc(verbose = FALSE, full = TRUE)
+      })
+      parallel::stopCluster(cl)
+    }
   })
   # clean up errors w/ proper formatting
   names(test_stats) <- genes
@@ -384,7 +387,7 @@ testDynamic <- function(expr.mat = NULL,
     total_time <- end_time - start_time
     total_time_units <- attributes(total_time)$units
     total_time_numeric <- as.numeric(total_time)
-    time_message <- paste0("testDynamic evaluated ",
+    time_message <- paste0("scLANE testing completed for ",
                            length(genes),
                            " genes across ",
                            n_lineages,

R/testSlope.R

@@ -29,7 +29,7 @@ testSlope <- function(test.dyn.res = NULL,
   if (is.null(test.dyn.res)) { stop("You forgot to provide results from testDynamic() to testSlope().") }
   # create table of results
   slope_df <- purrr::map_dfr(test.dyn.res, \(x) purrr::map_dfr(x, \(y) data.frame(y["MARGE_Slope_Data"][[1]]))) %>%
-              dplyr::arrange(P_Val, dplyr::desc(Test_Stat)) %>%
+              dplyr::arrange(P_Val, dplyr::desc(abs(Test_Stat))) %>%
               dplyr::mutate(P_Val_Adj = stats::p.adjust(P_Val, method = p.adj.method)) %>%
               dplyr::arrange(Gene, Breakpoint) %>%
               dplyr::mutate(Gene_Dynamic_Lineage_Slope = dplyr::if_else(P_Val_Adj < fdr.cutoff, 1, 0, missing = 0)) %>%

tests/testthat/test_scLANE.R

@@ -32,7 +32,7 @@ withr::with_output_sink(tempfile(), {
                                 genes = genes_to_test,
                                 n.potential.basis.fns = 5,
                                 size.factor.offset = cell_offset,
-                                n.cores = 2,
+                                parallel.exec = FALSE,
                                 track.time = TRUE)
   gee_gene_stats <- testDynamic(expr.mat = sim_data,
                                 pt = pt_test,