guard against small inputs in nCountRnaStratification

R/PseudoBulk.R

@@ -60,9 +60,21 @@ PseudobulkSeurat <- function(seuratObj,
                                                                              replace = F))
       }
       #drop NA column formed when initializing matrix
-      cluster_nCount_RNA_matrix <- cluster_nCount_RNA_matrix[,-1]
+      cluster_nCount_RNA_matrix <- cluster_nCount_RNA_matrix[, -1, drop = FALSE]
+      #guard against small input dataset sizes
+      if (ncol(cluster_nCount_RNA_matrix) < 2L || nrow(cluster_nCount_RNA_matrix) < 2L) {
+        print(paste0("Skipping nCount_RNA stratification for ", stratificationGroupingField, " due to small input dataset size. Minimum cluster size: ", minimum_cluster_size, ". Cluster nCount_RNA matrix dimensions: ", nrow(cluster_nCount_RNA_matrix), "x", ncol(cluster_nCount_RNA_matrix)))
+        next
+      }
       #compute KL divergences
-      KL_divergences <- flexmix::KLdiv(cluster_nCount_RNA_matrix)
+      KL_divergences <- tryCatch(
+        flexmix::KLdiv(cluster_nCount_RNA_matrix),
+        error = function(e) NULL
+      )
+      if (is.null(KL_divergences)) {
+        print(paste0("Skipping nCount_RNA stratification for ", stratificationGroupingField, " due to error in KL divergence calculation. Cluster nCount_RNA matrix dimensions: ", nrow(cluster_nCount_RNA_matrix), "x", ncol(cluster_nCount_RNA_matrix)))
+        next
+      }
 
       #perform rudimentary outlier detection on column sums of KL divergences assuming a half-normal distribution. 
       abnormal_RNA_clusters_index <- which(colSums(KL_divergences) >= mean(colSums(KL_divergences)) + 2*sd(colSums(KL_divergences)))

R/Seurat_III.R

@@ -512,6 +512,10 @@ ScoreCellCycle <- function(seuratObj, min.genes = 10, useAlternateG2M = FALSE) {
   s.genes <- s.genes[which(s.genes %in% rownames(seuratObj))]
   g2m.genes <- g2m.genes[which(g2m.genes %in% rownames(seuratObj))]
 
+  data_genes <- rownames(Seurat::GetAssayData(seuratObj, layer = 'data'))
+  s.genes <- intersect(s.genes, data_genes)
+  g2m.genes <- intersect(g2m.genes, data_genes)
+
   print(paste0("Genes present in seurat object: g2m (", length(g2m.genes), ") and s (", length(s.genes), ")"))
 
   if (length(g2m.genes) < min.genes || length(s.genes) < min.genes) {

tests/testthat/test-pseudobulk.R

@@ -32,6 +32,22 @@ test_that("Pseudobulk works", {
 
 })
 
+test_that("nCount RNA stratification skips when data is too small for KL divergence", {
+  seuratObj <- suppressWarnings(Seurat::UpdateSeuratObject(readRDS('../testdata/seuratOutput.rds')))
+  tiny <- seuratObj[, 1:2]
+  tiny$kl_test_cluster <- 0L
+  pseudo_tiny <- testthat::expect_no_error(
+    PseudobulkSeurat(
+      tiny,
+      groupFields = c('ClusterNames_0.2'),
+      nCountRnaStratification = TRUE,
+      stratificationGroupingFields = 'kl_test_cluster'
+    )
+  )
+  testthat::expect_false(is.null(pseudo_tiny))
+  testthat::expect_false('nCount_RNA_Stratification' %in% colnames([email protected]))
+})
+
 test_that("Pseudobulk-based differential expression works", {
   testthat::expect_equal(length(colnames(pbmc_small)), expected = 80) #weak insurance that pbmc_small doesn't change.
   [email protected][,"random_cohort"] <- base::rep(c(1,2), length.out = length(colnames(pbmc_small)))