bug fixes and added heatmap sorting fn -- closes #109

NAMESPACE

@@ -13,6 +13,7 @@ export(nbGAM)
 export(plotClusteredGenes)
 export(plotModels)
 export(smoothedCountsMatrix)
+export(sortGenesHeatmap)
 export(stripGLM)
 export(testDynamic)
 export(testSlope)
@@ -48,6 +49,7 @@ 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%")

R/fitGLMM.R

@@ -120,7 +120,7 @@ fitGLMM <- function(X_pred = NULL,
     } else {
       glmm_mod <- glmmTMB::glmmTMB(Y ~ X1 + X2 + X3 + X4 + (1 + X1 + X2 + X3 + X4 | subject),
                                    data = glmm_basis_df,
-                                   offset = log(Y.offset),
+                                   offset = log(1 / Y.offset),
                                    family = glmmTMB::nbinom2(link = "log"),
                                    se = TRUE,
                                    REML = FALSE)

R/sortGenesHeatmap.R

@@ -0,0 +1,48 @@
+#' Generate a table of fitted values and celltype metadata for genes of interest.
+#'
+#' @name sortGenesHeatmap
+#' @author Jack Leary
+#' @import magrittr
+#' @importFrom purrr map reduce
+#' @importFrom dplyr filter distinct select with_groups summarise arrange pull
+#' @description Sort genes such that genes with peak expression occurring earlier in pseudotime are first, and vice versa for genes with late peak expression. Useful for ordering genes in order to create heatmaps of expression cascades.
+#' @param heatmap.mat A matrix of raw or smoothed expression values with genes as columns and cells as rows. Defaults to NULL.
+#' @param pt.vec A numeric vector of pseudotime values for each cell i.e., for each row in the heatmap matrix. Defaults to NULL.
+#' @return A character vector of genes sorted by their peak expression values over pseudotime.
+#' @export
+#' @examples
+#' \dontrun{
+#' smoothed_counts <- smoothedCountsMatrix(gene_stats, pt = pt_df)
+#' sortGenesheatmap(heatmap.mat = smoothed_counts$Lineage_A,
+#'                  pt.vec = pt_df[!is.na(pt_df$Lineage_A), ]$Lineage_A)
+#' }
+
+sortGenesHeatmap <- function(heatmap.mat = NULL, pt.vec = NULL) {
+  # check inputs
+  if (!inherits(heatmap.mat, "matrix")) {
+    heatmap.mat <- try(as.matrix(heatmap.mat), silent = TRUE)
+    if (inherits(heatmap.mat, "try-error")) {
+      stop("heatmap.mat must be coerceable to a matrix.")
+    }
+  }
+  if (!is.numeric(pt.vec) || any(is.na(pt.vec))) { stop("pt.vec must be a numeric vector with no NA values.") }
+
+  # identify point at which peak expression occurs for each gene across pseudotime
+  gene_peak_order <- purrr::map(seq(ncol(heatmap.mat)), \(x) {
+    data.frame(gene = colnames(heatmap.mat)[x],
+               pt = pt.vec,
+               mRNA = heatmap.mat[, x]) %>%
+      dplyr::filter(mRNA == max(mRNA)) %>%
+      dplyr::distinct() %>%
+      dplyr::select(gene,
+                    pt,
+                    mRNA)
+  }) %>%
+    purrr::reduce(rbind) %>%
+    dplyr::with_groups(gene,
+                       dplyr::summarise,
+                       mu = mean(pt)) %>%
+    dplyr::arrange(mu) %>%
+    dplyr::pull(gene)
+  return(gene_peak_order)
+}

R/testDynamic.R

@@ -512,7 +512,7 @@ testDynamic <- function(expr.mat = NULL,
           }, silent = TRUE)
         } else {
           null_sumy_df <- try({
-            as.data.frame(broom.mixed::tidy(null_mod)) # saves a few bytes by converting from tibble
+            as.data.frame(broom.mixed::tidy(null_mod))  # saves a few bytes by converting from tibble
           }, silent = TRUE)
           null_pred_df <- try({
             data.frame(stats::predict(null_mod, type = "link", se.fit = TRUE)[1:2]) %>%
@@ -631,11 +631,11 @@ testDynamic <- function(expr.mat = NULL,
     total_time_numeric <- as.numeric(total_time)
     print(paste0("testDynamic evaluated ",
                  length(genes),
-                 " genes with ",
+                 " genes across ",
                   n_lineages,
                  " ",
                  ifelse(n_lineages == 1, "lineage ", "lineages "),
-                 "apiece in ",
+                 "in ",
                  round(total_time_numeric, 3),
                  " ",
                  total_time_units))

tests/testthat/test_scLANE.R

@@ -147,6 +147,8 @@ withr::with_output_sink(tempfile(), {
                                           size.factor.offset = cell_offset,
                                           parallel.exec = TRUE,
                                           n.cores = 2)
+  sorted_genes <- sortGenesHeatmap(heatmap.mat = smoothed_counts$Lineage_A,
+                                   pt.vec = pt_test$PT)
   fitted_values_table <- getFittedValues(test.dyn.res = glm_gene_stats,
                                          genes = names(glm_gene_stats),
                                          pt = pt_test,
@@ -276,6 +278,11 @@ test_that("smoothedCountsMatrix() output", {
   expect_type(smoothed_counts$Lineage_A, "double")
 })
 
+test_that("sortGenesHeatmap() output", {
+  expect_type(sorted_genes, "character")
+  expect_length(sorted_genes, ncol(smoothed_counts$Lineage_A))
+})
+
 test_that("getFittedValues() output", {
   expect_s3_class(fitted_values_table, "data.frame")
   expect_equal(ncol(fitted_values_table), 17)