diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index fc8c8f8..c9749ed 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -2,7 +2,7 @@ on: push: - branches: [main, master] + branches: main name: R-CMD-check diff --git a/.github/workflows/render-README.yaml b/.github/workflows/render-README.yaml index 7d373b9..9e51818 100644 --- a/.github/workflows/render-README.yaml +++ b/.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" diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index a1f48cb..8917e9c 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -2,7 +2,7 @@ on: push: - branches: [main, master] + branches: main name: test-coverage diff --git a/R/clusterGenes.R b/R/clusterGenes.R index 368ed58..e536323 100644 --- a/R/clusterGenes.R +++ b/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)) { diff --git a/R/createSlopeTestData.R b/R/createSlopeTestData.R index 1875d6f..339213d 100644 --- a/R/createSlopeTestData.R +++ b/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) diff --git a/R/getFittedValues.R b/R/getFittedValues.R index 4bcfba2..4d0fc48 100644 --- a/R/getFittedValues.R +++ b/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_, diff --git a/R/getKnotDist.R b/R/getKnotDist.R index 5b1547c..c88fc4b 100644 --- a/R/getKnotDist.R +++ b/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) } diff --git a/R/smoothedCountsMatrix.R b/R/smoothedCountsMatrix.R index 50b971b..bdcbd85 100644 --- a/R/smoothedCountsMatrix.R +++ b/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)) { diff --git a/R/testDynamic.R b/R/testDynamic.R index 4a8b762..0d3a24f 100644 --- a/R/testDynamic.R +++ b/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, diff --git a/R/testSlope.R b/R/testSlope.R index 0cb5d32..0544216 100644 --- a/R/testSlope.R +++ b/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)) %>% diff --git a/man/clusterGenes.Rd b/man/clusterGenes.Rd index 135c081..a0d770d 100644 --- a/man/clusterGenes.Rd +++ b/man/clusterGenes.Rd @@ -10,7 +10,7 @@ clusterGenes( size.factor.offset = NULL, clust.algo = "leiden", use.pca = FALSE, - n.PC = 15, + n.PC = 15L, lineages = NULL ) } diff --git a/tests/testthat/test_scLANE.R b/tests/testthat/test_scLANE.R index f07554b..ccc4508 100644 --- a/tests/testthat/test_scLANE.R +++ b/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,