---
title: "2. Benchmarking Anomaly Detection Against Ground Truth"
author:
    - name: Anthony Christidis
      affiliation:
        - &core_affiliation Core for Computational Biomedicine, Harvard Medical School
    - name: Andrew Ghazi
      affiliation:
        - *core_affiliation
    - name: Smriti Chawla
      affiliation:
        - *core_affiliation
    - name: Nitesh Turaga
      affiliation:
        - *core_affiliation
    - name: Ludwig Geistlinger
      affiliation:
        - *core_affiliation
    - name: Robert Gentleman
      affiliation:
        - *core_affiliation
package: scDiagnostics
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{2. Benchmarking Anomaly Detection Against Ground Truth}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options:
  markdown:
    wrap: 72
---

```{r setup, include = FALSE}
knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "ragg_png",
  dpi = 72,
  fig.retina = 2,
  fig.align = "center",
  out.width = "100%",
  pngquant = "--speed=1 --quality=1-5"
)
```

# Purpose

[Vignette 1](https://ccb-hms.github.io/scDiagnostics/articles/Introduction.html)
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 (`r BiocStyle::Biocpkg("scRNAseq")`), comparing `detectAnomaly()`
(Isolation Forest) against `calculateReconstructionError()` (PCA
reconstruction error), the package's two cell-type-specific anomaly
detection methods.

```{r, message=FALSE}
library(scDiagnostics)
library(SingleCellExperiment)
library(SingleR)
library(ggplot2)
library(dplyr)
```

# A live example: a withheld cell type related to a retained one

`zeisel_reference_data` and `zeisel_query_data` are a 70/30 split of the
Zeisel dataset (log-normalized, top 250 HVGs, PCA precomputed); see
`?zeisel_reference_data` for details.

```{r, message=FALSE}
data("zeisel_reference_data")
data("zeisel_query_data")

table(zeisel_reference_data$true_cell_type)
```

We withhold "pyramidal SS" from the reference entirely - a harder
detection problem than withholding a rare-but-distinct type, since
pyramidal SS is transcriptionally similar to the retained "pyramidal
CA1" type. We use `SingleR` to see where the query's true pyramidal SS
cells end up being mapped:

```{r, message=FALSE}
set.seed(1)
reference_missing <- zeisel_reference_data[, zeisel_reference_data$true_cell_type != "pyramidal SS"]
reference_missing <- scater::runPCA(reference_missing, ncomponents = 10)

pred <- SingleR(test = zeisel_query_data, ref = reference_missing,
                labels = reference_missing$true_cell_type)
zeisel_query_data$SingleR_annotation <- pred$labels

table(zeisel_query_data$SingleR_annotation[zeisel_query_data$true_cell_type == "pyramidal SS"])
```

Nearly all the true pyramidal SS cells get mapped to "pyramidal CA1".
Since we know which query cells are truly pyramidal SS, we can check
how well `detectAnomaly()` and `calculateReconstructionError()`
separate them from the correctly-labeled pyramidal CA1 cells they were
mapped alongside:

```{r, fig.height=5, fig.width=10}
target <- "pyramidal CA1"

anomaly_output <- detectAnomaly(
    reference_data = reference_missing, query_data = zeisel_query_data,
    ref_cell_type_col = "true_cell_type", query_cell_type_col = "SingleR_annotation",
    cell_types = target, n_hvgs = 30, pc_subset = 1:8, n_tree = 500)

reconstruction_output <- calculateReconstructionError(
    reference_data = reference_missing, query_data = zeisel_query_data,
    ref_cell_type_col = "true_cell_type", query_cell_type_col = "SingleR_annotation",
    cell_types = target, n_hvgs = 30, pc_subset = 1:8)

labels_target <- zeisel_query_data$true_cell_type[zeisel_query_data$SingleR_annotation == target]

data.frame(
    Method = c("detectAnomaly (Isolation Forest)", "calculateReconstructionError"),
    `True pyramidal SS flagged` = c(
        mean(anomaly_output[[target]]$query_anomaly[labels_target == "pyramidal SS"]),
        mean(reconstruction_output[[target]]$query_anomaly[labels_target == "pyramidal SS"])),
    `True pyramidal CA1 flagged` = c(
        mean(anomaly_output[[target]]$query_anomaly[labels_target == target]),
        mean(reconstruction_output[[target]]$query_anomaly[labels_target == target])),
    check.names = FALSE)
```

In this run, both methods flag most of the true pyramidal SS cells
while flagging a smaller fraction of the correctly-labeled pyramidal
CA1 cells - neither is perfect, which is exactly the situation where
combining them (see below) is worth considering. We can visualize the
Isolation Forest result; `data_type` shows one dataset per plot, so we
look at the reference (defining what "normal" looks like) and the query
(colored by anomaly status) side by side:

```{r, fig.height=5, fig.width=10}
plot(anomaly_output, cell_type = target, data_type = "reference", pc_subset = 1:3)
```

```{r, fig.height=5, fig.width=10}
plot(anomaly_output, cell_type = target, data_type = "query", pc_subset = 1:3)
```

This is one scenario (one withheld cell type, one reference/query
split, default settings beyond `n_hvgs`/`pc_subset`). The rest of this
vignette summarizes a systematic evaluation across many such scenarios,
computed once offline on the full (non-downsampled) Zeisel dataset; see
`?zeisel_benchmark_results` and `inst/script/ZeiselBenchmarkResults.R`
for the exact procedure.

```{r, message=FALSE}
data("zeisel_benchmark_results")
names(zeisel_benchmark_results)
```

# Baseline: distinct, related, and rare withheld cell types

Withholding a cell type that is transcriptionally distinct from
everything else (astrocytes) is an easier detection problem than
withholding one that closely resembles a retained cell type (pyramidal
SS, which is related to the retained pyramidal CA1), or one that is
simply rare (microglia):

```{r, fig.height=5, fig.width=8}
baseline <- zeisel_benchmark_results$gradients %>%
    filter(TestGroup == "Baseline") %>%
    mutate(Test = factor(Test, levels = c("Distinct (Astrocytes)",
                                          "Related (Pyramidal SS)",
                                          "Rare (Microglia)")))

ggplot(baseline, aes(x = Test, y = AUROC, fill = Method)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.7, color = "black") +
    geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
    coord_cartesian(ylim = c(0.4, 1)) +
    labs(x = "Missing cell type", y = "AUROC",
        title = "Baseline detection accuracy by missing cell type") +
    theme_bw()
```

Across these three scenarios, both methods reach a high AUROC, with
Isolation Forest at or above 0.94 in all three and Reconstruction Error
weakest on the rare (microglia) scenario in this particular run.

# 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).

```{r, fig.height=10, fig.width=10}
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:

```{r, fig.height=5, fig.width=9}
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()
```

```{r, fig.height=5, fig.width=9}
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:

```{r, message=FALSE}
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)
```

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 SessionInfo, echo=FALSE, message=FALSE, warning=FALSE, comment=NA}
options(width = 80)
sessionInfo()
```