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

## ----message=FALSE------------------------------------------------------------
library(scDiagnostics)
library(SingleCellExperiment)
library(SingleR)
library(ggplot2)
library(dplyr)

## ----message=FALSE------------------------------------------------------------
data("zeisel_reference_data")
data("zeisel_query_data")

table(zeisel_reference_data$true_cell_type)

## ----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"])

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

## ----fig.height=5, fig.width=10-----------------------------------------------
plot(anomaly_output, cell_type = target, data_type = "reference", pc_subset = 1:3)

## ----fig.height=5, fig.width=10-----------------------------------------------
plot(anomaly_output, cell_type = target, data_type = "query", pc_subset = 1:3)

## ----message=FALSE------------------------------------------------------------
data("zeisel_benchmark_results")
names(zeisel_benchmark_results)

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

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

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

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

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

## ----SessionInfo, echo=FALSE, message=FALSE, warning=FALSE, comment=NA--------
options(width = 80)
sessionInfo()