Vignette

Set up the environment

This vignette uses the following additional R packages. If any are not installed, run:

# In R
install.packages(
  c("Seurat", "ggplot2", "patchwork", "reticulate", "BiocManager"),
  repos = "https://cloud.r-project.org"
)
BiocManager::install("anndataR")

The Python environment is needed for ranking integration methods in step 3. To use scib-metrics from R, use reticulate to work with Python. If Conda is not installed, first install Miniconda following the official installation guide. Run the following commands in a terminal to create the Python environment and install the required packages:

# In a terminal
conda create -n bcmp-scib python=3.12
conda activate bcmp-scib
python -m pip install "pandas>=2.2.3,<3" anndata scib-metrics

1. Load a Seurat object

This bundled example Seurat object¹ contains raw counts, batch metadata, and a precomputed UMAP.

# In R
library(bcmp)
library(ggplot2)
library(patchwork)

seurat_path <- system.file(
  "extdata", "diabetic_kidney_lite_seurat.rds", package = "bcmp"
)
seurat <- readRDS(seurat_path)
seurat

An object of class Seurat
27980 features across 1500 samples within 1 assay
Active assay: RNA (27980 features, 0 variable features)
 1 layer present: counts
 12 dimensional reductions calculated: X_demo_umap, X_CSS, X_ComBat, X_Harmony,
 X_LIGER, X_Scanorama, X_Seurat_cca, X_Seurat_rpca, X_Uncorrected, X_fastMNN,
 X_scANVI, X_scVI

2. Run BCMP and visualize domains

Run BCMP

Simply run BCMP with the Seurat object and its batch key:

r_result <- bcmp(
  seurat,
  batch_key = "batch",
  verbose = FALSE
)

colnames(r_result$object[[]])

[1] "orig.ident"   "nCount_RNA"   "nFeature_RNA" "batch"        "bcmp_domain"

Visualize domains

umap <- as.data.frame(
  SeuratObject::Embeddings(r_result$object[["X_demo_umap"]])
)
colnames(umap) <- c("UMAP1", "UMAP2")

plot_umap <- function(labels, title, palette) {
  ggplot(transform(umap, label = as.character(labels)), aes(UMAP1, UMAP2, colour = label)) +
    geom_point() +
    scale_colour_brewer(palette = palette) +
    theme_void() + theme(legend.position = "bottom") +
    labs(title = title, colour = NULL)
}
plot_umap(r_result$object$batch, "Batch", "Accent") +
  plot_umap(r_result$object$bcmp_domain, "BCMP domain", "Paired")

3. Rank candidate integration methods

The eleven bundled candidate embeddings are stored as Seurat reductions. Their keys and display names are recorded in seurat@misc$integration_benchmark.

benchmark <- seurat@misc$integration_benchmark
candidate_reductions <- benchmark$embedding_keys
method_names <- benchmark$method_names

Prepare the scIB input

Because scib-metrics requires a Python AnnData object, use anndataR and reticulate to convert the Seurat object directly in memory. See the anndataR manual for conversion details.

reticulate::use_condaenv("bcmp-scib", required = TRUE)

reduction_mapping <- stats::setNames(
  as.list(candidate_reductions), candidate_reductions
)
adata <- reticulate::r_to_py(anndataR::as_AnnData(
  r_result$object,
  assay_name = "RNA",
  x_mapping = "counts",
  layers_mapping = FALSE,
  obs_mapping = c(batch = "batch", bcmp_domain = "bcmp_domain"),
  obsm_mapping = reduction_mapping,
  output_class = "ReticulateAnnData"
))

Calculate the ranking

Use the BCMP domain labels as the population labels for every candidate reduction.

benchmark_module <- reticulate::import("scib_metrics.benchmark", convert = FALSE)
benchmarker <- benchmark_module$Benchmarker(
  adata = adata,
  batch_key = "batch",
  label_key = "bcmp_domain",
  embedding_obsm_keys = candidate_reductions,
  bio_conservation_metrics = NULL,
  batch_correction_metrics = benchmark_module$BatchCorrection(),
  n_jobs = 1L,
  progress_bar = FALSE
)
benchmarker$benchmark()

scores <- reticulate::py_to_r(
  benchmarker$get_results(min_max_scale = TRUE)
)
ranking <- unlist(
  scores[candidate_reductions, "Batch correction"], use.names = FALSE
)
names(ranking) <- candidate_reductions
ranking <- sort(ranking, decreasing = TRUE)
names(ranking) <- unname(method_names[names(ranking)])
print(round(ranking, 2L))

     Harmony         scVI       scANVI   Seurat CCA  Seurat RPCA      fastMNN
        0.91         0.88         0.88         0.85         0.82         0.77
      ComBat    Scanorama        LIGER          CSS  Uncorrected
        0.70         0.65         0.58         0.57         0.22

A bar plot provides another view of the ranking:

ranking_data <- data.frame(
  method = factor(names(ranking), levels = rev(names(ranking))),
  score = unname(ranking)
)
ggplot(ranking_data, aes(method, score)) +
  geom_col(fill = "grey55") +
  coord_flip() +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = NULL, y = "scIB batch-correction score")