Fix fold change calculation for SCT assay

satijalab · Jul 6, 2023 · e08e108 · e08e108
1 parent a989986
commit e08e108
Show file tree

Hide file tree

Showing 2 changed files with 68 additions and 5 deletions.
diff --git a/R/differential_expression.R b/R/differential_expression.R
@@ -752,12 +752,18 @@ FindMarkers.SCTAssay <- function(
     'scale.data' = GetAssayData(object = object, slot = "counts"),
     numeric()
   )
-  if (is.null(x = mean.fxn)){
-    mean.fxn <- function(x) {
-      return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
-
-    }
+  # Default assumes the input is log1p(corrected counts)
+  default.mean.fxn <- function(x) {
+    return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
   }
+  mean.fxn <- mean.fxn %||% switch(
+    EXPR = slot,
+    'counts' = function(x) {
+      return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
+    },
+    'scale.data' = rowMeans,
+    default.mean.fxn
+  )
   fc.results <- FoldChange(
     object = object,
     slot = data.slot,
@@ -1102,6 +1108,9 @@ FoldChange.Assay <- function(
       default.mean.fxn
     ),
     'scale.data' = rowMeans,
+    'counts' = function(x) {
+      return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
+    },
     default.mean.fxn
   )
   # Omit the decimal value of e from the column name if base == exp(1)
@@ -1125,6 +1134,59 @@ FoldChange.Assay <- function(
   )
 }
 
+#' @importFrom Matrix rowMeans
+#' @rdname FoldChange
+#' @concept differential_expression
+#' @export
+#' @method FoldChange SCTAssay
+FoldChange.SCTAssay <- function(
+    object,
+    cells.1,
+    cells.2,
+    features = NULL,
+    slot = "data",
+    pseudocount.use = 1,
+    fc.name = NULL,
+    mean.fxn = NULL,
+    base = 2,
+    ...
+) {
+  pseudocount.use <- pseudocount.use %||% 1
+  data <- GetAssayData(object = object, slot = slot)
+  default.mean.fxn <- function(x) {
+    return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
+  }
+  mean.fxn <- mean.fxn %||% switch(
+    EXPR = slot,
+    'data' = default.mean.fxn,
+    'scale.data' = rowMeans,
+    'counts' = function(x) {
+      return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
+    },
+    default.mean.fxn
+  )
+  # Omit the decimal value of e from the column name if base == exp(1)
+  base.text <- ifelse(
+    test = base == exp(1),
+    yes = "",
+    no = base
+  )
+  fc.name <- fc.name %||% ifelse(
+    test = slot == "scale.data",
+    yes = "avg_diff",
+    no = paste0("avg_log", base.text, "FC")
+  )
+  FoldChange(
+    object = data,
+    cells.1 = cells.1,
+    cells.2 = cells.2,
+    features = features,
+    mean.fxn = mean.fxn,
+    fc.name = fc.name
+  )
+}
+
+
 #' @importFrom Matrix rowMeans
 #' @rdname FoldChange
 #' @concept differential_expression

diff --git a/tests/testthat/test_differential_expression.R b/tests/testthat/test_differential_expression.R
@@ -300,6 +300,7 @@ test_that("FindAllMarkers works as expected", {
 
   # CLR normalization
   expect_equal(results.clr[1, "p_val"], 1.209462e-11)
+
   expect_equal(results.clr[1, "avg_log2FC"], -1.079924, tolerance = 1e-6)
   expect_equal(results.clr[1, "pct.1"], 0.083)
   expect_equal(results.clr[1, "pct.2"], 0.909)