Merge pull request #136 from jr-leary7/dev

Dev

DESCRIPTION

@@ -2,8 +2,8 @@ Package: scLANE
 Type: Package
 Title: Model gene expression dynamics with spline-based NB GLMs, GEEs, & GLMMs
 Version: 0.7.3
-Authors@R: c(person(given = "Jack", family = "Leary", email = "j.leary@ufl.edu", role = c("aut", "cre")), 
-             person(given = "Rhonda", family = "Bacher", email = "rbacher@ufl.edu", role = c("ctb", "fnd")))
+Authors@R: c(person(given = "Jack", family = "Leary", email = "j.leary@ufl.edu", role = c("aut", "cre"), comment = c(ORCID = "0009-0004-8821-3269")), 
+             person(given = "Rhonda", family = "Bacher", email = "rbacher@ufl.edu", role = c("ctb", "fnd"), comment = c(ORCID = "0000-0001-5787-476X")))
 Description: This package uses truncated power basis spline models to build flexible, interpretable models of single cell gene expression over pseudotime or latent time. 
              The modeling architectures currently supported are negative binomial GLMs, GEEs, & GLMMs. 
              Downstream analysis functionalities include model comparison, dynamic gene clustering, smoothed counts generation, gene set enrichment testing, & visualization. 

NAMESPACE

@@ -7,6 +7,7 @@ export(enrichDynamicGenes)
 export(extractBreakpoints)
 export(fitGLMM)
 export(getFittedValues)
+export(getKnotDist)
 export(getResultsDE)
 export(marge2)
 export(modelLRT)
@@ -48,9 +49,11 @@ importFrom(dplyr,ntile)
 importFrom(dplyr,pull)
 importFrom(dplyr,relocate)
 importFrom(dplyr,rename)
+importFrom(dplyr,rowwise)
 importFrom(dplyr,sample_frac)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
+importFrom(dplyr,ungroup)
 importFrom(dplyr,with_groups)
 importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
@@ -84,10 +87,10 @@ importFrom(glmmTMB,nbinom2)
 importFrom(parallel,clusterEvalQ)
 importFrom(parallel,clusterExport)
 importFrom(parallel,clusterSetRNGStream)
-importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
 importFrom(purrr,discard)
+importFrom(purrr,imap)
 importFrom(purrr,map)
 importFrom(purrr,map2)
 importFrom(purrr,map_chr)

R/createCellOffset.R

@@ -14,6 +14,7 @@
 #' @export
 #' @examples
 #' \dontrun{
+#' createCellOffset(expr.mat = sce_obj)
 #' createCellOffset(expr.mat = counts(sce_obj))
 #' createCellOffset(expr.mat = seu_obj, scale.factor = 1e5)
 #' }

R/getKnotDist.R

@@ -0,0 +1,34 @@
+#' Pull the set of knots for dynamic genes across each lineage.
+#'
+#' @name getKnotDist
+#' @author Jack Leary
+#' @import magrittr
+#' @importFrom purrr imap reduce
+#' @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.
+#' @return A data.frame containing gene name, lineage ID, and knot location in pseudotime.
+#' @export
+#' @examples
+#' \dontrun{
+#' getKnotDist(gene_stats)
+#' }
+
+getKnotDist <- function(test.dyn.res = NULL, dyn.genes = NULL) {
+  # check inputs
+  if (is.null(test.dyn.res)) { stop("You forgot one of the arguments to getKnotDist().") }
+  if (is.null(dyn.genes)) {
+    dyn.genes <- names(test.dyn.res)
+  }
+  # 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)
+  return(knot_df)
+}

R/summarizeModel.R

@@ -23,9 +23,15 @@ summarizeModel <- function(marge.model = NULL, pt = NULL) {
                      Slope.Segment = NA_real_,
                      Trend.Segment = NA_real_)
   } else {
-    # extract model equation & slopes
+    # extract model equation, slopes, & covariances
+    if (marge.model$model_type == "GEE") {
+      coef_vcov <- as.matrix(marge.model$final_mod$var)
+    } else {
+      coef_vcov <- vcov(marge.model$final_mod)
+    }
     coef_df <- data.frame(coef_name = names(coef(marge.model$final_mod)),
-                          coef_value = unname(coef(marge.model$final_mod)))
+                          coef_value = unname(coef(marge.model$final_mod)),
+                          coef_var = unname(diag(coef_vcov)))
 
     coef_df <- coef_df[-which(coef_df$coef_name == "Intercept"), ]
 

R/testDynamic.R

@@ -8,7 +8,7 @@
 #' @importFrom bigstatsr as_FBM
 #' @importFrom foreach foreach %dopar% registerDoSEQ
 #' @importFrom doParallel registerDoParallel
-#' @importFrom parallel makeCluster detectCores stopCluster clusterEvalQ clusterExport clusterSetRNGStream
+#' @importFrom parallel makeCluster stopCluster clusterEvalQ clusterExport clusterSetRNGStream
 #' @importFrom withr with_output_sink
 #' @importFrom MASS glm.nb negative.binomial theta.mm
 #' @importFrom dplyr rename mutate relocate
@@ -93,6 +93,9 @@ testDynamic <- function(expr.mat = NULL,
   if (is.null(expr.mat) || is.null(pt)) { stop("You forgot some inputs to testDynamic().") }
 
   # get raw counts from SingleCellExperiment or Seurat object & transpose to cell x gene dense matrix
+  if (is.null(genes)) {
+    genes <- rownames(expr.mat)
+  }
   if (inherits(expr.mat, "SingleCellExperiment")) {
     expr.mat <- BiocGenerics::counts(expr.mat)[genes, ]
     expr.mat <- as.matrix(expr.mat)
@@ -102,15 +105,12 @@ testDynamic <- function(expr.mat = NULL,
                                      assay = Seurat::DefaultAssay(expr.mat))
     expr.mat <- as.matrix(expr.mat[genes, ])
   } else if (inherits(expr.mat, "dgCMatrix")) {
-    expr.mat <- as.matrix(expr.mat)
+    expr.mat <- as.matrix(expr.mat[genes, ])
+  } else {
+    expr.mat <- expr.mat[genes, ]
   }
   if (!(inherits(expr.mat, "matrix") || inherits(expr.mat, "array"))) { stop("Input expr.mat must be coerceable to a matrix of integer counts.") }
   expr.mat <- t(expr.mat)  # transpose to cell x gene matrix
-  if (is.null(genes)) {
-    genes <- colnames(expr.mat)
-  } else {
-    expr.mat <- expr.mat[, genes]
-  }
 
   # extract pseudotime dataframe if input is results from Slingshot
   if (inherits(pt, "SlingshotDataSet")) {
@@ -165,8 +165,7 @@ testDynamic <- function(expr.mat = NULL,
   package_list <- c("glm2", "scLANE", "MASS",  "bigstatsr", "broom.mixed", "dplyr", "stats")
   if (is.gee) {
     package_list <- c(package_list, "geeM")
-  }
-  if (is.glmm) {
+  } else if (is.glmm) {
     package_list <- c(package_list, "glmmTMB")
   }
 

inst/rmarkdown/templates/Bacher_Group_HTML/skeleton/skeleton.Rmd

@@ -1,6 +1,6 @@
 ---
 title: "Title"
-subtitle: "University of Florida - Dept. of Biostatistics - Bacher Group"
+subtitle: "UF Dept. of Biostatistics - Bacher Group"
 author: "Name"
 date: "`r Sys.Date()`"
 output:
@@ -20,29 +20,40 @@ knitr::opts_chunk$set(echo = TRUE,
                       comment = NA, 
                       message = FALSE, 
                       warning = FALSE, 
-                      fig.align = "center")
+                      fig.align = "center", 
+                      dev = "png", 
+                      dpi = 300)
 set.seed(312)  # lucky seed
 ```
 
 # Libraries
 
 ```{r}
-library(glm2)
-library(mgcv)
 library(dplyr)
-library(scran)
 library(scLANE)
-library(scater)
 library(ggplot2)
-library(tradeSeq)
-library(slingshot)
-library(doParallel)
-library(kableExtra)
-library(SingleCellExperiment)
 ```
 
 # Data 
 
+```{r}
+
+```
+
 # Analysis 
 
+```{r}
+
+```
+
 # Conclusions 
+
+```{r}
+
+```
+
+# Session info 
+
+```{r}
+sessioninfo::session_info()
+```

tests/testthat/test_scLANE.R

@@ -34,7 +34,7 @@ withr::with_output_sink(tempfile(), {
                                 size.factor.offset = cell_offset,
                                 n.cores = 2,
                                 track.time = TRUE)
-  gee_gene_stats <- testDynamic(sim_data,
+  gee_gene_stats <- testDynamic(expr.mat = sim_data,
                                 pt = pt_test,
                                 genes = genes_to_test,
                                 n.potential.basis.fns = 5,
@@ -53,7 +53,7 @@ withr::with_output_sink(tempfile(), {
                                  glmm.adaptive = TRUE,
                                  id.vec = sim_data$subject,
                                  n.cores = 2,
-                                 track.time = TRUE)
+                                 track.time = FALSE)
   # get results tables overall
   glm_test_results <- getResultsDE(glm_gene_stats)
   gee_test_results <- getResultsDE(gee_gene_stats)
@@ -191,6 +191,7 @@ withr::with_output_sink(tempfile(), {
   gsea_res <- enrichDynamicGenes(glm_test_results, species = "hsapiens")
   coef_summary_glm <- summarizeModel(marge_mod_offset, pt = pt_test)
   coef_summary_gee <- summarizeModel(marge_mod_GEE_offset, pt = pt_test)
+  knot_df <- getKnotDist(glm_gene_stats)
 })
 
 # run tests
@@ -359,3 +360,8 @@ test_that("summarizeModels() output", {
   expect_type(coef_summary_glm$Slope.Segment, "double")
   expect_type(coef_summary_gee$Slope.Segment, "double")
 })
+
+test_that("getKnotDist() output", {
+  expect_s3_class(knot_df, "data.frame")
+  expect_equal(ncol(knot_df), 3)
+})