--- title: "Immune-inspired prototypes for cytometry and single-cell data with bHIVE" author: name: Nick Borcherding email: ncborch@gmail.com affiliation: Washington University in St. Louis, School of Medicine, St. Louis, MO, USA date: "`r Sys.Date()`" output: BiocStyle::html_document: toc_float: true package: bHIVE vignette: > %\VignetteIndexEntry{Introduction to bHIVE} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(message = FALSE, tidy = FALSE, fig.align = "center") library(BiocStyle) fmt <- function(x) formatC(x, format = "f", digits = 2) ``` # Introduction Cytometry and single-cell RNA-seq measure many markers on many cells. Two questions come up in almost every analysis. Which cell populations are present? And which population does a new cell belong to? Clustering answers the first. Classification, or label transfer from an annotated reference, answers the second. `r Biocpkg("bHIVE")` answers both with one model, an artificial immune network. A set of prototypes, called antibodies, evolves the way B cells do in an immune response. Antibodies that bind the data well are cloned, the clones mutate, and near-duplicates are suppressed. What remains is a small set of points in marker space that summarize the cells. What this gives an analyst: - **Prototypes you can read.** Each antibody is a point in the original marker units, so a heatmap of the antibodies reads like a table of population phenotypes. - **One model for clustering and label transfer.** The same fitted object that clusters a sample can label new cells with `$predict()`. - **Rare populations.** In supervised mode every class is seeded with its own antibodies, so a population of a few cells is not lost. In unsupervised mode `coverageBoost` seeds antibodies where coverage is poor. - **A cluster count you choose or let emerge.** By default the number of clusters comes out of the suppression dynamics. Set `targetK` to fix it. - **Composable immune modules.** The R6 engine, `AINet`, accepts modules for somatic hypermutation, germinal center selection, microenvironment and more. The second vignette, `vignette("immune-modules", package = "bHIVE")`, covers them. bHIVE works with Bioconductor data structures. `bHIVE()` accepts a `SummarizedExperiment` or `SingleCellExperiment` and stores its labels in `colData()`. `bHIVEParam()` plugs bHIVE into `r Biocpkg("bluster")`, so it also runs through `bluster::clusterRows()` and `bluster::clusterSweep()`. bHIVE does not replace graph-based clustering or `r Biocpkg("FlowSOM")`. The comparison below shows that a shared-nearest-neighbor graph recovers manually gated populations better on this dataset. bHIVE is a prototype model to inspect and reuse. # Installation ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("bHIVE") ``` # A mass cytometry dataset We use bone marrow mass cytometry data from Levine et al. (2015), available through `r Biocpkg("HDCytoData")`. Each cell was gated by hand into one of 14 populations or left unassigned. The gates give a reference to score clusterings against. HDCytoData stores cells in rows and markers in columns. Bioconductor single-cell containers use the opposite orientation, so we transpose. We keep the 32 lineage markers, take 5,000 gated cells, and apply the usual `asinh(x / 5)` transform for mass cytometry. ```{r data} library(bHIVE) library(SingleCellExperiment) library(HDCytoData) library(bluster) library(ggplot2) se <- Levine_32dim_SE() type_markers <- colData(se)$marker_class == "type" gated <- which(rowData(se)$population_id != "unassigned") set.seed(1) keep <- sample(gated, 5000) sce <- SingleCellExperiment( assays = list(exprs = t(asinh(assay(se, "exprs")[keep, type_markers] / 5))), colData = droplevels(rowData(se)[keep, ])) colnames(sce) <- paste0("cell", seq_len(ncol(sce))) sce sort(table(sce$population_id), decreasing = TRUE) ``` The smallest populations have only a few dozen cells in this sample. # Clustering a SingleCellExperiment `bHIVE()` clusters the cells (the columns) of a `SingleCellExperiment` using an assay, here `exprs`, or a `reducedDim()` embedding through `use.dimred`. For scRNA-seq, cluster a PCA embedding rather than a sparse assay. Two settings matter for high-dimensional data. `affinityParams = list(alpha = "auto")` sets the bandwidth of the Gaussian affinity from the data. A fixed bandwidth of 1 makes almost every affinity zero in 32 dimensions. `targetK` asks for a fixed number of clusters, here the 14 gated populations. ```{r cluster} set.seed(42) sce <- bHIVE(sce, assay.type = "exprs", targetK = 14, affinityParams = list(alpha = "auto"), maxIter = 20, verbose = FALSE) table(sce$bHIVE) ari_42 <- pairwiseRand(sce$population_id, sce$bHIVE, mode = "index") ari_42 ``` The labels are in `colData(sce)$bHIVE`. The full result, including the fitted model, is in `metadata(sce)$bHIVE`. The last line is the adjusted Rand index (ARI) against the manual gates. An ARI of 1 means perfect agreement, and 0 means chance. This run reaches `r fmt(ari_42)`. The comparison below shows how much that changes with the random seed. ## Reading the prototypes Each cluster has one prototype in marker space. A heatmap of the prototypes shows each cluster's phenotype. We label each cluster with the gated population most of its cells came from. ```{r prototypes, fig.width = 9, fig.height = 4.5} res <- metadata(sce)$bHIVE protos <- res$antibodies_unscaled majority <- tapply(as.character(sce$population_id), sce$bHIVE, function(p) names(which.max(table(p)))) rownames(protos) <- paste0(seq_len(nrow(protos)), ": ", majority) long <- data.frame( cluster = factor(rep(rownames(protos), ncol(protos)), levels = rev(rownames(protos))), marker = factor(rep(colnames(protos), each = nrow(protos)), levels = colnames(protos)), value = as.vector(protos)) ggplot(long, aes(marker, cluster, fill = value)) + geom_tile() + scale_fill_viridis_c(name = "asinh") + labs(x = NULL, y = NULL) + theme_minimal() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) ``` ```{r prototypes-summary, include = FALSE} n_rows <- table(majority) split_pops <- names(n_rows)[n_rows > 1] merged_pops <- setdiff(levels(sce$population_id), majority) list_or_none <- function(x) if (length(x)) paste(x, collapse = ", ") else "none" ``` Reading the rows shows what a run separates and what it merges. A population on more than one row was split. This run splits `r list_or_none(split_pops)`. A population with no row of its own was merged into its neighbors. Here that is `r list_or_none(merged_pops)`. That explains most of the gap between this run and the gates. ## Plotting the cells with their prototypes `visualizeHIVE()` fits a PCA to the cells and projects the prototypes into the same coordinates. ```{r visualize, fig.width = 7, fig.height = 5} X <- t(assay(sce, "exprs")) visualizeHIVE(res, X = X, title = "bHIVE prototypes on the PCA of the cells") ``` # bHIVE through bluster `bHIVEParam()` describes a bHIVE run for `r Biocpkg("bluster")`. It works anywhere a `BlusterParam` does. ```{r bluster} param <- bHIVEParam(targetK = 14, affinityParams = list(alpha = "auto"), maxIter = 20) param set.seed(42) clusters <- clusterRows(X, param) table(clusters) ``` The rows of `X` are the cells of `sce`, so the labels go straight into `colData()`. ```{r collabels} colLabels(sce) <- clusters nlevels(colLabels(sce)) ``` `bluster::clusterSweep()` varies any setting of the parameter object. Here it tries three cluster counts. ```{r sweep} set.seed(42) sweep <- clusterSweep(X, bHIVEParam(affinityParams = list(alpha = "auto"), maxIter = 10), targetK = c(8, 14, 20)) vapply(as.list(sweep$clusters), function(cl) pairwiseRand(sce$population_id, cl, mode = "index"), numeric(1)) ``` # Comparison with FlowSOM and bluster methods We compare bHIVE with three methods used for cytometry and single-cell data: - `r Biocpkg("FlowSOM")` builds a self-organizing map and merges its nodes into metaclusters. It is the standard prototype method for cytometry. We ask for 14 metaclusters. - k-means with 14 centers, through `bluster::KmeansParam()`. - A shared-nearest-neighbor graph with Louvain community detection, through `bluster::SNNGraphParam()`. This is the usual default for single-cell data. It picks its own number of clusters. bHIVE depends on its random seed, so we run it three times. ```{r compare} truth <- sce$population_id timed <- function(expr) { start <- proc.time()[["elapsed"]] clusters <- expr list(clusters = clusters, seconds = proc.time()[["elapsed"]] - start) } runs <- list( "bHIVE, seed 1" = timed({ set.seed(1); clusterRows(X, param) }), "bHIVE, seed 2" = timed({ set.seed(2); clusterRows(X, param) }), "bHIVE, seed 3" = timed({ set.seed(3); clusterRows(X, param) }), "FlowSOM" = timed(FlowSOM::GetMetaclusters(FlowSOM::FlowSOM( X, compensate = FALSE, transform = FALSE, scale = FALSE, colsToUse = seq_len(ncol(X)), nClus = 14, seed = 1))), "k-means" = timed({ set.seed(1); clusterRows(X, KmeansParam(centers = 14)) }), "SNN graph + Louvain" = timed(clusterRows( X, SNNGraphParam(k = 15, cluster.fun = "louvain")))) comparison <- data.frame( method = names(runs), clusters = vapply(runs, function(r) length(unique(r$clusters)), integer(1)), ARI = vapply(runs, function(r) pairwiseRand(truth, r$clusters, mode = "index"), numeric(1)), seconds = vapply(runs, `[[`, numeric(1), "seconds"), row.names = NULL) knitr::kable(comparison, digits = c(0, 0, 3, 1)) ``` ```{r compare-summary, include = FALSE} ari <- setNames(comparison$ARI, comparison$method) n_cl <- setNames(comparison$clusters, comparison$method) bhive_ari <- c(ari[grepl("^bHIVE", names(ari))], ari_42) ``` The graph-based method reaches an ARI of `r fmt(ari[["SNN graph + Louvain"]])` with only `r n_cl[["SNN graph + Louvain"]]` clusters. ARI is dominated by the large populations, so merging a few small ones costs little. FlowSOM reaches `r fmt(ari[["FlowSOM"]])` and k-means `r fmt(ari[["k-means"]])`. Across the three seeds here and the run above, bHIVE ranges from `r fmt(min(bhive_ari))` to `r fmt(max(bhive_ari))`. A single bHIVE run can match the graph method or fall well below FlowSOM. That seed dependence is the main practical difference. If the goal is one stable partition of a dataset, FlowSOM or a graph-based method is the safer choice. With bHIVE, run several seeds and compare them, for example with `bluster::pairwiseRand()`. The `ConvergentSelector` module keeps the antibodies that recur across runs. FlowSOM also builds prototypes, its map nodes, and can map new samples onto them. So prototypes alone are not what sets bHIVE apart. The differences are how bHIVE learns them, by selection, mutation and suppression rather than on a fixed grid, and that the same model also learns from labels. The next section shows that. # Label transfer to new cells Label transfer trains a classifier on annotated cells and applies it to new ones. We train bHIVE on half of the gated cells and predict the other half. For classification, bHIVE seeds every class with at least `minPerClass` antibodies (default 5), so rare populations get their own antibodies. ```{r transfer} set.seed(2) train <- sample(ncol(sce), ncol(sce) / 2) test <- setdiff(seq_len(ncol(sce)), train) set.seed(3) ref <- bHIVE(sce[, train], y = "population_id", assay.type = "exprs", nAntibodies = 150, maxIter = 20, verbose = FALSE) model <- metadata(ref)$bHIVE$model predicted <- model$predict(X[test, ]) ``` As a baseline we use the nearest class centroid. Each population is one mean vector, and each test cell takes the label of the closest mean. ```{r baseline} centroids <- apply(X[train, ], 2, function(v) tapply(v, truth[train], mean)) dist2 <- vapply(seq_len(nrow(centroids)), function(j) colSums((t(X[test, ]) - centroids[j, ])^2), numeric(length(test))) baseline <- rownames(centroids)[max.col(-dist2)] score <- function(p) { hit <- p == truth[test] c(accuracy = mean(hit), balanced_accuracy = mean(tapply(hit, truth[test], mean))) } knitr::kable(rbind(bHIVE = score(predicted), "nearest centroid" = score(baseline)), digits = 3) ``` Balanced accuracy averages the recall of each population, so rare populations count as much as common ones. ```{r recall, fig.width = 8, fig.height = 4.5} recall <- function(p) tapply(p == truth[test], truth[test], mean) per_pop <- data.frame( population = rep(levels(truth), 2), method = rep(c("bHIVE", "nearest centroid"), each = nlevels(truth)), recall = c(recall(predicted), recall(baseline))) per_pop$population <- factor(per_pop$population, levels = names(sort(table(truth)))) ggplot(per_pop, aes(population, recall, fill = method)) + geom_col(position = "dodge") + coord_flip() + scale_fill_viridis_d(end = 0.8) + labs(x = NULL, y = "Recall on held-out cells") + theme_minimal() ``` ```{r transfer-summary, include = FALSE} acc <- rbind(bHIVE = score(predicted), centroid = score(baseline)) r_bhive <- recall(predicted) r_base <- recall(baseline) worst <- names(which.min(r_bhive)) ``` The nearest-centroid baseline beats bHIVE here. Its accuracy is `r fmt(acc["centroid", "accuracy"])` against `r fmt(acc["bHIVE", "accuracy"])`, and its balanced accuracy `r fmt(acc["centroid", "balanced_accuracy"])` against `r fmt(acc["bHIVE", "balanced_accuracy"])`. It matches or beats bHIVE on `r sum(r_base >= r_bhive)` of `r length(r_bhive)` populations. bHIVE's per-class seeding keeps every population represented, and its lowest recall is `r fmt(min(r_bhive))` for `r worst`. The gated populations in this dataset are compact and well separated in marker space, so one mean per population is enough. bHIVE spends several antibodies on each population. That helps when a population is spread out or has more than one mode, which is not the case here. For label transfer on data like these, a simple centroid classifier is the better choice. The fitted model can label any cells measured on the same panel. The cells the original analysis left unassigned are one example. ```{r unassigned} unassigned <- which(rowData(se)$population_id == "unassigned") set.seed(4) new_cells <- asinh(assay(se, "exprs")[sample(unassigned, 2000), type_markers] / 5) head(sort(table(model$predict(new_cells)), decreasing = TRUE)) ``` These cells were not gated by hand, so these labels have no reference to check them against. Treat them as hypotheses to check against the markers. # The R6 engine and immune modules `bHIVE()` is a thin wrapper over the R6 class `AINet`. Using `AINet` directly gives a reusable model object and access to the immune modules. Here `coverageBoost` adds antibodies where the data are covered worst, before the 14 clusters are formed. ```{r ainet} net <- AINet$new(nAntibodies = 40, maxIter = 15, targetK = 14, coverageBoost = TRUE, affinityParams = list(alpha = "auto"), verbose = FALSE) set.seed(5) net$fit(X) net pairwiseRand(truth, net$result$assignments, mode = "index") ``` `AINet` also accepts modules for somatic hypermutation (`SHMEngine`), idiotypic regulation (`IdiotypicNetwork`), germinal center selection (`GerminalCenter`), microenvironment-dependent mutation (`Microenvironment`), V(D)J initialization (`VDJLibrary`), two-signal activation (`ActivationGate`), memory (`MemoryPool`) and isotype switching (`ClassSwitcher`). The second vignette demonstrates each one: ```{r modules-vignette, eval = FALSE} vignette("immune-modules", package = "bHIVE") ``` # Other tools in the package - `honeycombHIVE()` stacks bHIVE layers. Each layer clusters the prototypes of the layer below, which gives a hierarchy of populations. - `swarmbHIVE()` runs a grid search over any bHIVE setting and scores each run with an internal clustering index or a classification metric. - `bHIVEmodel` and `honeycombHIVEmodel` are custom models for `r CRANpkg("caret")`, for cross-validated tuning of bHIVE classifiers. ```{r other-tools} set.seed(6) layers <- honeycombHIVE(X, task = "clustering", layers = 2, nAntibodies = 40, maxIter = 10, affinityParams = list(alpha = "auto"), verbose = FALSE) vapply(layers, function(l) nrow(l$prototypes), integer(1)) set.seed(7) tuned <- swarmbHIVE(X, task = "clustering", grid = data.frame(targetK = c(10, 14, 20)), metric = "davies_bouldin", maxIter = 10, verbose = FALSE) tuned$results ``` # Session information ```{r session} sessionInfo() ```