Purpose

Automated cell type annotation - transferring labels from a reference dataset onto a new query dataset - is now a routine step in single-cell RNA-seq (scRNA-seq) analysis. It is fast and reproducible, but it is only as trustworthy as the alignment between reference and query: batch effects, cell states missing from the reference, or systematic differences in sequencing depth can all produce confidently-labeled cells that are, in fact, misannotated.

scDiagnostics provides diagnostics for exactly this problem. Rather than another annotation method, it is a toolkit for auditing annotations you already have: is the query well-aligned with the reference in PCA space? Are there cells that look anomalous relative to their assigned cell type? If so, which genes distinguish them from the reference? Answering these questions helps decide whether an annotation transfer can be trusted, and where to look if it cannot.

The package operates on SingleCellExperiment objects and is not specific to scRNA-seq: the same diagnostics apply directly to spatial data stored as a r BiocStyle::Biocpkg("SpatialExperiment") or r BiocStyle::Biocpkg("SpatialFeatureExperiment") object, without any modification, since both extend SingleCellExperiment. Vignette 4 demonstrates this on MERFISH spatial data.

Overview of the workflow

Across the case-study vignettes in this package (see below), the same three-step diagnostic pattern recurs:

1. Project the query data onto the reference’s PCA space and visualize how each cell type compares (e.g. plotCellTypePCA(), projectPCA()).
2. Detect cells whose projection looks anomalous relative to their assigned reference cell type (detectAnomaly()).
3. Characterize what makes the anomalous cells different, by testing for expression shifts in genes that drive the relevant principal components (calculateGeneShifts()).

The four panels below illustrate this on the COVID-19 case study described in vignette 3: CD14 monocytes from a severe COVID-19 query project further along PC1 and PC2 than the healthy reference (panel A); a subset of those cells is flagged as anomalous by detectAnomaly() (panel B); and calculateGeneShifts() shows that the anomalous cells specifically over-express a panel of interferon-response genes, both as a heatmap of per-gene z-scores (panel C) and as fold-changes relative to the reference (panel D).

These specific results are from the COVID-19 case study and should not be read as a general property of every dataset - see vignette 2 for how detection accuracy was benchmarked against known ground truth more broadly.

Installation

BiocManager::install("ccb-hms/scDiagnostics")

BiocManager::install("ccb-hms/scDiagnostics",
                     build_vignettes = TRUE,
                     dependencies = TRUE)

library(scDiagnostics)

A worked example with known ground truth

Before applying detectAnomaly() to real data (as in the case-study vignettes), it helps to see it work on data where the “right answer” is known. We simulate two batches of cells with the splatter package, three cell types each, and treat one batch as the reference and the other as the query.

library(splatter)
library(scuttle)
library(scater)
library(SingleR)

set.seed(100)

# Simulate two batches of 500 cells, 3 balanced cell types
sce_ref <- mockSCE()
params <- splatEstimate(sce_ref)
params <- setParams(
    params,
    batchCells = c(500, 500), batch.facLoc = 0, batch.facScale = 0,
    group.prob = c(1 / 3, 1 / 3, 1 / 3),
    de.prob = c(0.1, 0.2, 0.2),
    de.facLoc = c(0.250, 0.375, 0.375),
    de.facScale = c(0.2, 0.3, 0.4),
    out.prob = 0, out.facLoc = 4, out.facScale = 0.5)
simulated_data <- splatSimulate(params, method = "groups", verbose = FALSE)

# Treat Batch1 as reference, Batch2 as query
reference_data <- simulated_data[, simulated_data$Batch == "Batch1"]
query_data <- simulated_data[, simulated_data$Batch == "Batch2"]

reference_data$Cell_Type <- factor(reference_data$Group)
levels(reference_data$Cell_Type) <- c("Cell Type A", "Cell Type B", "Cell Type C")
query_data$Cell_Type <- factor(query_data$Group)
levels(query_data$Cell_Type) <- c("Cell Type A", "Cell Type B", "Cell Type C")

reference_data <- logNormCounts(reference_data)
#> Warning in .library_size_factors(assay(x, assay.type), ...): 'librarySizeFactors' is deprecated.
#> Use 'scrapper::centerSizeFactors' instead.
#> See help("Deprecated")
#> Warning in .local(x, ...): 'normalizeCounts' is deprecated.
#> Use 'scrapper::normalizeCounts' instead.
#> See help("Deprecated")
query_data <- logNormCounts(query_data)
#> Warning in .library_size_factors(assay(x, assay.type), ...): 'librarySizeFactors' is deprecated.
#> Use 'scrapper::centerSizeFactors' instead.
#> See help("Deprecated")
#> Warning in .library_size_factors(assay(x, assay.type), ...): 'normalizeCounts' is deprecated.
#> Use 'scrapper::normalizeCounts' instead.
#> See help("Deprecated")

When the reference contains all three cell types, SingleR recovers the ground-truth labels essentially perfectly:

reference_data <- runPCA(reference_data, ncomponents = 10)
query_data <- runPCA(query_data, ncomponents = 10)

pred <- SingleR(query_data, reference_data, labels = reference_data$Cell_Type)
query_data$SingleR_annotation <- pred$labels

mean(query_data$SingleR_annotation == query_data$Cell_Type)
#> [1] 1

Now suppose “Cell Type C” is missing from the reference entirely - a common real-world scenario where a cell state present in the query simply was not sampled in the reference. SingleR is forced to assign those cells to the closest remaining type:

reference_missing <- reference_data[, reference_data$Cell_Type != "Cell Type C"]
reference_missing <- runPCA(reference_missing, ncomponents = 10)

pred_missing <- SingleR(query_data, reference_missing,
                        labels = reference_missing$Cell_Type)
query_data$SingleR_annotation_missing <- pred_missing$labels

# Where do the true Cell Type C cells get misannotated to?
table(query_data$SingleR_annotation_missing[query_data$Cell_Type == "Cell Type C"])
#> 
#> Cell Type A Cell Type B 
#>         173           5

Most of the true “Cell Type C” cells are misannotated as “Cell Type A”. Because we know which cells are truly misannotated in this simulation, we can check whether detectAnomaly() actually flags them as such, relative to the correctly-annotated “Cell Type A” cells:

anomaly_output <- detectAnomaly(
    reference_data = reference_missing,
    query_data = query_data,
    ref_cell_type_col = "Cell_Type",
    query_cell_type_col = "SingleR_annotation_missing",
    cell_types = "Cell Type A",
    pc_subset = 1:2,
    n_tree = 1000,
    threshold_method = "absolute",
    anomaly_threshold = 0.5)

is_anomalous <- anomaly_output[["Cell Type A"]]$query_anomaly
labels_a <- query_data$Cell_Type[query_data$SingleR_annotation_missing == "Cell Type A"]

# Fraction flagged as anomalous, split by true identity
tapply(is_anomalous, labels_a, mean)
#> Cell Type A Cell Type B Cell Type C 
#>  0.08695652          NA  0.47398844

In this run, roughly half of the truly-misannotated “Cell Type C” cells are flagged as anomalous, compared to a small fraction of the correctly labeled “Cell Type A” cells - detectAnomaly() is picking up a real signal here, though far from a perfect separation. We can visualize the same result:

plot(anomaly_output, cell_type = "Cell Type A", data_type = "query", pc_subset = 1:2)

For a systematic evaluation of how detection accuracy holds up across label noise, class imbalance, and batch effects - rather than this one simulated example - see vignette 2.

Finding the function you need

scDiagnostics groups its functions into five broad categories. Functions marked with a vignette link are walked through in more depth there; the rest are documented in the reference manual (e.g. ?detectAnomaly). The full, finer-grained reference index is also available if you’d rather browse by a more specific task.

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] SingleR_2.15.4              scater_1.41.2              
 [3] ggplot2_4.0.3               scuttle_1.23.2             
 [5] splatter_1.37.1             SingleCellExperiment_1.35.2
 [7] SummarizedExperiment_1.43.0 Biobase_2.73.2             
 [9] GenomicRanges_1.65.4        Seqinfo_1.3.2              
[11] IRanges_2.47.5              S4Vectors_0.51.10          
[13] BiocGenerics_0.59.12        generics_0.1.4             
[15] MatrixGenerics_1.25.0       matrixStats_1.5.0          
[17] 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 dplyr_1.2.1          
 [91] gtable_0.3.6          isotree_0.6.1-5       sass_0.4.10          
 [94] digest_0.6.39         SparseArray_1.13.2    ggrepel_0.9.8        
 [97] dqrng_0.4.1           rjson_0.2.23          farver_2.1.2         
[100] htmltools_0.5.9       lifecycle_1.0.5       GlobalOptions_0.1.4  
[103] statmod_1.5.2         MASS_7.3-66