diff --git a/DESCRIPTION b/DESCRIPTION
index 8d91eab..d04390d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: monocle3
 Title: Clustering, Differential Expression, and Trajectory Analysis for
     Single-Cell RNA-Seq
-Version: 1.4.11
+Version: 1.4.12
 Authors@R: c(
     person(given = "Hannah",
            family = "Pliner",
@@ -99,6 +99,7 @@ Imports:
     RcppHNSW (>= 0.3.0),
     Rtsne (>= 0.15),
     S4Vectors,
+    sf,
     shiny,
     slam (>= 0.1-45),
     spdep (>= 1.1-2),
@@ -116,7 +117,6 @@ Suggests:
     pryr (>= 0.1.4),
     rmarkdown,
     scran,
-    sf,
     spelling,
     testthat (>= 2.1.0)
 Remotes:
diff --git a/R/cluster_genes.R b/R/cluster_genes.R
index f4a3c95..efdaefb 100644
--- a/R/cluster_genes.R
+++ b/R/cluster_genes.R
@@ -536,5 +536,10 @@ aggregate_gene_expression <- function(cds,
   if (exclude.na){
     agg_mat <- agg_mat[row.names(agg_mat) != "NA", colnames(agg_mat) != "NA",drop=FALSE]
   }
+
+  if(is(agg_mat, 'IterableMatrix')) {
+    agg_mat <- as(agg_mat, 'dgCMatrix')
+  }
+
   return(agg_mat)
 }
diff --git a/R/find_markers.R b/R/find_markers.R
index 4a5f16a..65e15f9 100644
--- a/R/find_markers.R
+++ b/R/find_markers.R
@@ -74,7 +74,7 @@ top_markers <- function(cds,
                         cores=1,
                         verbose=FALSE) {
 
-  if(is(counts(cds), 'IterableMatrix') && is.null(counts_row_order)) {
+  if(is(counts(cds), 'IterableMatrix') && is.null(counts_row_order(cds))) {
     stop(paste('This CDS has a BPCells counts matrix but no counts_row_order matrix, which',
                'top_markers() requires. Use the command',
                 '  cds <- set_cds_row_order_matrix(cds=cds)',
@@ -121,6 +121,8 @@ top_markers <- function(cds,
                                                            norm_method="size_only",
                                                            scale_agg_values=FALSE))
 
+# bge the cluster_binary_exprs and cluster_mean_exprs pairs appear to be the same when run with dgCMatrix vs BPCells matrix
+
   if (verbose)
     message("Computing Jensen-Shannon specificities")
 
@@ -364,7 +366,7 @@ test_marker_for_cell_group = function(gene_id, cell_group, cell_group_df, cds,
     }
     else {
       f_expression <-
-        log(as.numeric(as(counts_row_order(cds)[gene_id,], 'dgCMatrix')) / size_factors(cds) + 0.1)
+        log(as.numeric(as(monocle3::counts_row_order(cds)[gene_id,], 'dgCMatrix')) / size_factors(cds) + 0.1)
     }
 
     #print(sum(SingleCellExperiment::counts(cds)[gene_id,] > 0))
diff --git a/R/generics.R b/R/generics.R
index 7326646..557f6b3 100644
--- a/R/generics.R
+++ b/R/generics.R
@@ -652,24 +652,24 @@ setMethod("counts<-", signature(object="SingleCellExperiment"),
 )
 
 
-#' Generic to access cds row order BPCells count matrix.
+#' Generic to access cds row order BPCells counts matrix.
 #' @param x A cell_data_set object.
 #' 
 #' @examples
 #'  \donttest{
 #'    cds <- load_a549()
-#'    exprs(cds)
+#'    mat_row_order <- counts_row_order(cds)
 #'  }         
 #'          
-#' @return BPCells row order ount matrix.
+#' @return BPCells row order counts matrix.
 #'
 #' @export
 setGeneric("counts_row_order", function(x) standardGeneric("counts_row_order"))
   
-#' Method to access cds row order BPCells count matrix
+#' Method to access cds row order BPCells counts matrix
 #' @param x A cell_data_set object.
 #'
-#' @return BPCells row order count matrix.
+#' @return BPCells row order counts matrix.
 #' 
 #' @export
 setMethod("counts_row_order", "cell_data_set", function(x) {
@@ -681,6 +681,11 @@ setMethod("counts_row_order", "cell_data_set", function(x) {
       stop('CDS has no BPCells row order counts matrix')
     }
   }
+  if(get_global_variable('bpcells_matrix_pair_check')) {
+    if(!check_bpcells_counts_matrix_pair(x)) {
+      stop('')
+    }
+  }
   value <- assay(x, 'counts_row_order')
   return(value)
 })                                 
diff --git a/R/graph_test.R b/R/graph_test.R
index 2f4ad5d..a65a524 100644
--- a/R/graph_test.R
+++ b/R/graph_test.R
@@ -154,15 +154,17 @@ graph_test <- function(cds,
     message("Performing Moran's I test: ...")
   }
 
-  exprs_mat <- SingleCellExperiment::counts(cds)
+   # Use row major order BPCells count matrix.
+  if(is(counts(cds), 'IterableMatrix')) {
+    exprs_mat <- monocle3::counts_row_order(cds)
+  }
+  else {
+    exprs_mat <- SingleCellExperiment::counts(cds)
+  }
+
   exprs_mat <- exprs_mat[, attr(lw, "region.id"), drop=FALSE]
   sz <- size_factors(cds)[attr(lw, "region.id")]
 
-  # Use row major order BPCells count matrix.
-  if(is(exprs_mat, 'IterableMatrix')) {
-    exprs_mat <- counts_row_order(cds)
-  }
-
   wc <- spdep::spweights.constants(lw, zero.policy = TRUE, adjust.n = TRUE)
   test_res <- pbmcapply::pbmclapply(row.names(exprs_mat),
                                     FUN = function(x, sz, alternative,
diff --git a/R/matrix.R b/R/matrix.R
index 7575739..9878fdc 100644
--- a/R/matrix.R
+++ b/R/matrix.R
@@ -878,6 +878,38 @@ set_cds_row_order_matrix <- function(cds) {
 }
 
 
+# Test (partially) that the BPCells counts() and
+# counts_row_order() matrices are consistent.
+# The test checks that the two matrices exist and
+# that row sums are the same for the two.
+check_bpcells_counts_matrix_pair <- function(cds) {
+  if(!is(assays(cds)[['counts']], 'IterableMatrix')) {
+    message('Error: the cds does not have a BPCells counts matrix.')
+    return(FALSE)
+  }
+
+  if(!is(assays(cds)[['counts_row_order']], 'IterableMatrix')) {
+    message('Error: the cds does not have a BPCells counts_row_order matrix.')
+    return(FALSE)
+  }
+
+  counts_rowsums <- BPCells::rowSums(counts(cds))
+  counts_row_order_rowsums <- BPCells::rowSums(assays(cds)[['counts_row_order']])
+
+  if(length(counts_rowsums) != length(counts_row_order_rowsums)) {
+    message('Error: the cds counts and counts_row_order matrix row sums have different lengths')
+    return(FALSE)
+  }
+
+  if(any(counts_rowsums != counts_row_order_rowsums)) {
+    message('Error: the cds counts and counts_row_order matrix row sums have different values')
+    return(FALSE)
+  }
+
+  return(TRUE)
+}
+
+
 #' Convert the counts matrix class in the given CDS.
 #'
 #' @description Converts the counts matrix that is in the
diff --git a/R/utils.R b/R/utils.R
index 503580d..e6363a6 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -221,6 +221,7 @@ sparse_par_c_apply <- function (cl = NULL, x, FUN, convert_to_dense, ...) {
 mc_es_apply <- function(cds, MARGIN, FUN, required_packages, cores=1,
                         convert_to_dense=TRUE,
                         reduction_method="UMAP", ...) {
+# message('mc_es_apply: start')
   parent <- environment(FUN)
   if (is.null(parent))
     parent <- emptyenv()
@@ -287,12 +288,22 @@ mc_es_apply <- function(cds, MARGIN, FUN, required_packages, cores=1,
   }
 
   #
-  # 20230424 bge: These two following calls using counts(cds) appear to work based on my simple test.
+  # 20230424 bge: The following sparse_par_?_apply calls using counts(cds) appear to work based on my simple test.
   #
   if (MARGIN == 1){
-    suppressWarnings(res <- sparse_par_r_apply(cl=cl, x=SingleCellExperiment::counts(cds), FUN=FUN,
-                                               convert_to_dense=convert_to_dense, ...))
-  }else{
+# message('mc_es_apply: MARGIN 1')
+    if( is(counts(cds), 'IterableMatrix')) {
+# message('mc_es_apply: BPCells matrix')
+      suppressWarnings(res <- sparse_par_r_apply(cl=cl, x=monocle3::counts_row_order(cds), FUN=FUN,
+                                                 convert_to_dense=convert_to_dense, ...))
+    }
+    else {
+# message('mc_es_apply: dgCMatrix matrix')
+      suppressWarnings(res <- sparse_par_r_apply(cl=cl, x=SingleCellExperiment::counts(cds), FUN=FUN,
+                                                 convert_to_dense=convert_to_dense, ...))
+    }
+  }
+  else {
     suppressWarnings(res <- sparse_par_c_apply(cl=cl, x=SingleCellExperiment::counts(cds), FUN=FUN,
                                                convert_to_dense=convert_to_dense, ...))
   }
@@ -303,6 +314,7 @@ mc_es_apply <- function(cds, MARGIN, FUN, required_packages, cores=1,
 #' @importFrom Biobase multiassign
 smart_es_apply <- function(cds, MARGIN, FUN, convert_to_dense,
                            reduction_method="UMAP", ...) {
+# message('smart_es_apply: start')
   parent <- environment(FUN)
   if (is.null(parent))
     parent <- emptyenv()
@@ -317,9 +329,24 @@ smart_es_apply <- function(cds, MARGIN, FUN, convert_to_dense,
                        as.data.frame(coldata_df), envir=e1)
   environment(FUN) <- e1
 
-  if (is_sparse_matrix(SingleCellExperiment::counts(cds))){
+  if (is(SingleCellExperiment::counts(cds), 'IterableMatrix')) {
+# message('smart_es_apply: BPCells matrix')
+    if(MARGIN == 1) {
+# message('smart_es_apply: MARGIN 1')
+      res <- sparse_apply(monocle3::counts_row_order(cds), MARGIN, FUN, convert_to_dense, ...)
+    }
+    else {
+# message('smart_es_apply: MARGIN 2')
+      res <- sparse_apply(SingleCellExperiment::counts(cds), MARGIN, FUN, convert_to_dense, ...)
+    }
+  }
+  else
+  if (is_sparse_matrix(SingleCellExperiment::counts(cds))) {
+# message('smart_es_apply: dgCMatrix matrix')
     res <- sparse_apply(SingleCellExperiment::counts(cds), MARGIN, FUN, convert_to_dense, ...)
-  } else {
+  }
+  else {
+# message('smart_es_apply: dense matrix')
     res <- pbapply::pbapply(SingleCellExperiment::counts(cds), MARGIN, FUN, ...)
   }
 
@@ -596,7 +623,7 @@ combine_cds <- function(cds_list,
     # up to this pass through the loop.
     exp <- counts(cds_list[[i]])
     if(bpcells_matrix_flag && !is(exp, 'IterableMatrix')) {
-      exp <- as(exp, 'IterableMatrix')
+      exp <- as(exp, 'IterableMatrix')       # wraps dgCMatrix in IterableMatrix
     }
     exp <- exp[intersect(row.names(exp), gene_list),, drop=FALSE]
 
@@ -661,7 +688,7 @@ combine_cds <- function(cds_list,
 
       # Append additional rows.
       if(bpcells_matrix_flag) {
-        exp <- rbind2(exp, as(extra_rows, 'IterableMatrix'))
+        exp <- rbind2(exp, as(extra_rows, 'IterableMatrix'))       # wraps dgCMatrix in IterableMatrix
       }
       else {
         exp <- rbind(exp, extra_rows)
@@ -889,7 +916,7 @@ combine_cds_for_maddy <- function(cds_list,
     # up to this pass through the loop.
     exp <- counts(cds_list[[i]])
     if(bpcells_matrix_flag && !is(exp, 'IterableMatrix')) {
-      exp <- as(exp, 'IterableMatrix')
+      exp <- as(exp, 'IterableMatrix')                               # wraps dgCMatrix in IterableMatrix
     }
     exp <- exp[intersect(row.names(exp), gene_list),, drop=FALSE]
 
@@ -954,7 +981,7 @@ combine_cds_for_maddy <- function(cds_list,
 
       # Append additional rows.
       if(bpcells_matrix_flag) {
-        exp <- rbind2(exp, as(extra_rows, 'IterableMatrix'))
+        exp <- rbind2(exp, as(extra_rows, 'IterableMatrix'))    # wraps dgCMatrix in IterableMatrix
       }
       else {
         exp <- rbind(exp, extra_rows)
diff --git a/R/zzz.R b/R/zzz.R
index b2ea304..0ae6bde 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -94,6 +94,7 @@ get_global_variable <- function(variable_name=NULL) {
   set_global_variable('monocle3_timer_t0', 0)
   set_global_variable('monocle3_timer_msg', "")
   set_global_variable('monocle_gc_matrix_path', list())
+  set_global_variable('bpcells_matrix_pair_check', TRUE)
 
   # Default nn_control list for functions that do not need
   # an index, which is all but the label transfer functions.