2. Benchmarking Anomaly Detection Against Ground Truth

Purpose

Vignette 1 showed detectAnomaly() correctly flagging a subset of misannotated cells in one simulated example. A single example only goes so far: how well does this actually work in general, and how sensitive is it to the choices you have to make - which cell type is missing, how much label noise is in the reference, how imbalanced the cell types are, or how large a batch effect separates query from reference?

This vignette works through those questions on the Zeisel mouse brain dataset (scRNAseq), comparing detectAnomaly() (Isolation Forest) against calculateReconstructionError() (PCA reconstruction error), the package’s two cell-type-specific anomaly detection methods.

library(scDiagnostics)
library(SingleCellExperiment)
library(SingleR)
library(ggplot2)
library(dplyr)

Sensitivity to label noise, class imbalance, and batch effects

For the “related” (pyramidal SS) and “rare” (microglia) scenarios, the benchmark also varies three conditions independently: the fraction of reference labels randomly shuffled (label noise), the number of cells retained in the mapped-to reference cluster (class imbalance), and a mean expression shift applied to 20% of query genes (a stand-in for a batch effect).

gradients <- zeisel_benchmark_results$gradients %>%
    filter(TestGroup %in% c("Related", "Rare"))

plot_gradient <- function(df, test_name, x_lab, decreasing_x = FALSE) {
    sub_df <- df %>% filter(Test == test_name)
    sub_df$X_Value <- if (decreasing_x) {
        factor(sub_df$X_Value, levels = sort(as.numeric(unique(sub_df$X_Value)), decreasing = TRUE))
    } else {
        as.numeric(sub_df$X_Value)
    }
    ggplot(sub_df, aes(x = X_Value, y = AUROC, color = Method, group = Method)) +
        geom_line() + geom_point(size = 2) +
        geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
        coord_cartesian(ylim = c(0.4, 1)) +
        facet_wrap(~TestGroup) +
        labs(x = x_lab, y = "AUROC", title = test_name) +
        theme_bw()
}

noise_plot <- plot_gradient(gradients, "Noise", "Fraction of reference labels shuffled")
imbalance_plot <- plot_gradient(gradients, "Imbalance", "Cells in mapped-to reference cluster", decreasing_x = TRUE)
batch_plot <- plot_gradient(gradients, "Batch", "Mean expression shift applied to query")

noise_plot

imbalance_plot

batch_plot

In this benchmark, both methods stay well above the AUROC = 0.5 no-skill baseline across the full range of label noise and class imbalance tested. Under an increasing batch effect, Isolation Forest tends to hold up better than Reconstruction Error in the rare (microglia) scenario - consistent with the idea that a tree-based method partitioning on individual PCs can be more robust to a systematic shift than a global reconstruction-error metric, though this is a pattern observed in this specific benchmark rather than a general guarantee.

Hyperparameter sensitivity

Both methods require choices: how many HVGs or PCs to use, and what threshold marks a cell as anomalous. The benchmark also grid-searches these choices for a single scenario (pyramidal SS withheld, mapped to pyramidal CA1):

Each point is one hyperparameter configuration (a choice of PCs or HVGs, and a threshold rule); splitting the grid into one panel per feature space (for detectAnomaly()) or per MAD threshold (for calculateReconstructionError()) keeps each panel to a handful of points:

ggplot(zeisel_benchmark_results$if_tuning,
      aes(x = Specificity, y = Sensitivity, color = Threshold)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", color = "gray70") +
    geom_vline(xintercept = 0.8, linetype = "dashed", color = "gray70") +
    geom_point(size = 3, alpha = 0.85) +
    facet_wrap(~Mode) +
    coord_cartesian(xlim = c(0.6, 1), ylim = c(0.3, 1)) +
    labs(title = "detectAnomaly(): sensitivity vs. specificity across hyperparameters") +
    theme_bw()

ggplot(zeisel_benchmark_results$re_tuning,
      aes(x = Specificity, y = Sensitivity, color = HVGs)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", color = "gray70") +
    geom_vline(xintercept = 0.8, linetype = "dashed", color = "gray70") +
    geom_point(size = 3, alpha = 0.85) +
    facet_wrap(~MAD_Threshold) +
    coord_cartesian(xlim = c(0.6, 1), ylim = c(0.3, 1)) +
    labs(title = "calculateReconstructionError(): sensitivity vs. specificity across hyperparameters") +
    theme_bw()

(The dashed lines mark 80% sensitivity/specificity as a rough visual reference, not a formal threshold.) Within each panel, points still vary by PC subset or HVG count - the full per-configuration breakdown is in zeisel_benchmark_results$if_tuning/re_tuning if you want to identify a specific one.

In this grid, no single configuration dominates on both sensitivity and specificity simultaneously (the usual precision/recall trade-off). For this particular scenario, configurations using a small, cell-type-targeted set of HVGs with a MAD-based threshold tend to land closer to the top-right (high sensitivity and specificity) corner - but that is a property of this benchmark, not a universal ranking of hyperparameters, and a different dataset could favor a different configuration. n_hvgs = 30 with a MAD-based threshold (the defaults used earlier in this vignette) is a reasonable starting point rather than a claim that it is optimal in general; it’s worth re-checking against your own data if detection accuracy matters a lot for your use case.

Combining both methods

detectAnomaly() and calculateReconstructionError() look at different parts of the data and can fail in different ways, which makes them candidates for use together rather than as competitors:

  • detectAnomaly() partitions cells directly along the retained principal components - it flags cells that sit in an unusual location within that low-dimensional PC subspace.
  • calculateReconstructionError() does the opposite in a sense: it compresses each cell to that same low-dimensional subspace and back, and flags cells whose original expression profile isn’t well reconstructed - i.e. it is sensitive to signal in the subspace orthogonal to the retained PCs (loosely, the null space of the PCA projection), which Isolation Forest never looks at directly.

Because they emphasize different subspaces, flagging a cell as anomalous whenever either method flags it (the union of the two) can catch cells that one method misses but the other doesn’t - raising sensitivity beyond what either method achieves alone, at the cost of more false positives. In the pyramidal SS example above, neither method alone is perfect, and their errors don’t fully overlap:

if_flag <- anomaly_output[[target]]$query_anomaly
re_flag <- reconstruction_output[[target]]$query_anomaly
union_flag <- if_flag | re_flag

data.frame(
    Rule = c("Isolation Forest only", "Reconstruction Error only",
            "Either flags (union)"),
    `True pyramidal SS flagged` = c(
        mean(if_flag[labels_target == "pyramidal SS"]),
        mean(re_flag[labels_target == "pyramidal SS"]),
        mean(union_flag[labels_target == "pyramidal SS"])),
    `True pyramidal CA1 flagged` = c(
        mean(if_flag[labels_target == target]),
        mean(re_flag[labels_target == target]),
        mean(union_flag[labels_target == target])),
    check.names = FALSE)
#>                        Rule True pyramidal SS flagged
#> 1     Isolation Forest only                 0.7830189
#> 2 Reconstruction Error only                 0.8962264
#> 3      Either flags (union)                 0.9528302
#>   True pyramidal CA1 flagged
#> 1                  0.2332016
#> 2                  0.1264822
#> 3                  0.3003953

In this run, the union flags more true pyramidal SS cells than either method alone - each method catches some cells the other misses. That gain isn’t free: the union also flags more of the correctly-labeled pyramidal CA1 cells than either method alone, since it inherits every false positive from both. Whether that trade-off is worth it (versus requiring both methods to agree, which pushes the other way - fewer false positives, but only the anomalies both methods happen to catch) depends on whether missing a real anomaly or chasing a false one is more costly for your analysis. Neither combination rule is “correct” in general, and this result is specific to this scenario, not a claim that the union always beats each method individually.

R Session Info

R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 26.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] dplyr_1.2.1                 SpatialExperiment_1.23.0   
 [3] SingleR_2.15.4              scater_1.41.2              
 [5] ggplot2_4.0.3               scuttle_1.23.2             
 [7] splatter_1.37.1             SingleCellExperiment_1.35.2
 [9] SummarizedExperiment_1.43.0 Biobase_2.73.2             
[11] GenomicRanges_1.65.4        Seqinfo_1.3.2              
[13] IRanges_2.47.5              S4Vectors_0.51.10          
[15] BiocGenerics_0.59.12        generics_0.1.4             
[17] MatrixGenerics_1.25.0       matrixStats_1.5.0          
[19] scDiagnostics_1.7.15        BiocStyle_2.41.0           

loaded via a namespace (and not attached):
  [1] gridExtra_2.3.1       rlang_1.3.0           magrittr_2.0.5       
  [4] clue_0.3-68           GetoptLong_1.1.1      otel_0.2.0           
  [7] ggridges_0.5.7        compiler_4.6.1        png_0.1-9            
 [10] systemfonts_1.3.2     vctrs_0.7.3           shape_1.4.6.1        
 [13] crayon_1.5.3          pkgconfig_2.0.3       fastmap_1.2.0        
 [16] backports_1.5.1       magick_2.9.1          XVector_0.53.0       
 [19] labeling_0.4.3        rmarkdown_2.32        ggbeeswarm_0.7.3     
 [22] ragg_1.5.2            purrr_1.2.2           xfun_0.61            
 [25] bluster_1.23.1        cachem_1.1.0          beachmat_2.29.2      
 [28] jsonlite_2.0.0        DelayedArray_0.39.6   BiocParallel_1.47.0  
 [31] irlba_2.3.7           parallel_4.6.1        cluster_2.1.8.3      
 [34] R6_2.6.1              bslib_0.12.0          RColorBrewer_1.1-3   
 [37] limma_3.99.0          GGally_2.4.0          jquerylib_0.1.4      
 [40] Rcpp_1.1.2            iterators_1.0.14      knitr_1.52           
 [43] splines_4.6.1         Matrix_1.7-6          igraph_2.3.3         
 [46] tidyselect_1.2.1      abind_1.4-8           yaml_2.3.12          
 [49] viridis_0.6.5         doParallel_1.0.17     codetools_0.2-20     
 [52] lattice_0.23-1        tibble_3.3.1          withr_3.0.3          
 [55] S7_0.2.2              evaluate_1.0.5        survival_3.8-12      
 [58] ggstats_0.14.0        fitdistrplus_1.2-6    circlize_0.4.18      
 [61] pillar_1.11.1         BiocManager_1.30.27   checkmate_2.3.4      
 [64] foreach_1.5.2         scales_1.4.0          RhpcBLASctl_0.23-42  
 [67] glue_1.8.1            metapod_1.21.0        maketools_1.3.2      
 [70] tools_4.6.1           BiocNeighbors_2.7.3   sys_3.4.3            
 [73] ScaledMatrix_1.21.0   locfit_1.5-9.12       buildtools_1.0.0     
 [76] scran_1.41.1          grid_4.6.1            tidyr_1.3.2          
 [79] colorspace_2.1-3      edgeR_4.99.6          beeswarm_0.4.0       
 [82] BiocSingular_1.29.1   vipor_0.4.7           cli_3.6.6            
 [85] rsvd_1.0.5            textshaping_1.0.5     S4Arrays_1.13.1      
 [88] viridisLite_0.4.3     ComplexHeatmap_2.29.0 gtable_0.3.6         
 [91] isotree_0.6.1-5       sass_0.4.10           digest_0.6.39        
 [94] SparseArray_1.13.2    ggrepel_0.9.8         dqrng_0.4.1          
 [97] rjson_0.2.23          farver_2.1.2          htmltools_0.5.9      
[100] lifecycle_1.0.5       GlobalOptions_0.1.4   statmod_1.5.2        
[103] MASS_7.3-66