## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7,
    fig.height = 5,
    message = FALSE,
    warning = FALSE
)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("scBatchQC")

## ----quickstart---------------------------------------------------------------
library(scBatchQC)
library(SingleCellExperiment)

set.seed(123)
n_genes <- 500

# Batch 1: high-depth fresh tissue (10x v3)
counts_b1 <- matrix(
    rpois(n_genes * 100, lambda = 12),
    nrow = n_genes, ncol = 100
)

# Batch 2: lower-depth cryopreserved (10x v2)
# with 5 deliberately damaged cells
counts_b2 <- matrix(
    rpois(n_genes * 100, lambda = 5),
    nrow = n_genes, ncol = 100
)
counts_b2[, 1:5] <- floor(counts_b2[, 1:5] / 10)

# Combine and label
counts <- cbind(counts_b1, counts_b2)
rownames(counts) <- paste0("Gene", seq_len(n_genes))
rownames(counts)[1:30] <- paste0("MT-", seq_len(30))
colnames(counts) <- paste0("Cell", seq_len(200))

sce <- SingleCellExperiment(assays = list(counts = counts))
sce$batch <- rep(c("Batch1_v3", "Batch2_v2"), each = 100)
sce

## ----qc_metrics---------------------------------------------------------------
sce <- batchAwareQCMetrics(
    sce,
    batch = "batch",
    nmads = 3,
    shrink_strength = 0.5
)

# Columns added
grep("^scBatchQC", names(colData(sce)), value = TRUE)

## ----qc_table-----------------------------------------------------------------
# Outlier counts per batch
table(Outlier = sce$scBatchQC_outlier, Batch = sce$batch)

## ----plot_qc, fig.height = 8, fig.cap = "QC distributions per batch. Dashed red lines = harmonised thresholds. Red points = flagged outliers."----
plotBatchQC(sce, batch = "batch")

## ----doublet_rates------------------------------------------------------------
cells_loaded <- c(Batch1_v3 = 8000, Batch2_v2 = 5000)

sce <- estimateBatchDoubletRate(
    sce,
    batch = "batch",
    cells_loaded = cells_loaded,
    protocol = c(Batch1_v3 = "10x_v3", Batch2_v2 = "10x_v2")
)

# Estimated doublet rate per batch
tapply(sce$scBatchQC_doublet_rate, sce$batch, unique)

## ----threshold_sweep----------------------------------------------------------
sweep <- lapply(c(2, 2.5, 3, 3.5, 4), function(n) {
    r <- harmonizeQCThresholds(sce, batch = "batch", nmads = n)
    data.frame(
        nmads = n,
        flagged = sum(rowSums(as.data.frame(r$n_flagged)))
    )
})
do.call(rbind, sweep)

## ----filter-------------------------------------------------------------------
sce_clean <- sce[, !sce$scBatchQC_outlier]
cat(
    "Kept", ncol(sce_clean), "of", ncol(sce),
    "cells (removed", sum(sce$scBatchQC_outlier), ")\n"
)

## ----bqc_result---------------------------------------------------------------
library(S4Vectors)

result <- BQCResult(
    qcFlags = DataFrame(outlier = sce$scBatchQC_outlier),
    doubletScores = sce$scBatchQC_doublet_rate,
    batchSummary = estimateBatchDoubletRate(
        sce,
        batch = "batch",
        cells_loaded = cells_loaded,
        return_sce = FALSE
    )
)
result

## ----accessors----------------------------------------------------------------
head(qcFlags(result))
head(doubletScores(result))
batchSummary(result)

## ----full_workflow, eval = FALSE----------------------------------------------
# # Step 1: batch-aware QC
# sce <- batchAwareQCMetrics(sce, batch = "batch", nmads = 3)
# 
# # Step 2: doublet rates
# sce <- estimateBatchDoubletRate(
#     sce,
#     batch = "batch",
#     cells_loaded = my_cells_loaded
# )
# 
# # Step 3: visualise
# plotBatchQC(sce, batch = "batch")
# 
# # Step 4: filter
# sce_clean <- sce[, !sce$scBatchQC_outlier]
# 
# # Step 5: downstream (scran, Seurat, Harmony, etc.)

## ----session------------------------------------------------------------------
sessionInfo()