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.

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 bluster, so it also runs through bluster::clusterRows() and bluster::clusterSweep().

bHIVE does not replace graph-based clustering or 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

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

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

## class: SingleCellExperiment 
## dim: 32 5000 
## metadata(0):
## assays(1): exprs
## rownames(32): CD45RA CD133 ... CD64 CD41
## rowData names(0):
## colnames(5000): cell1 cell2 ... cell4999 cell5000
## colData names(2): patient_id population_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

sort(table(sce$population_id), decreasing = TRUE)

## 
##            CD4_T_cells              Monocytes            CD8_T_cells 
##                   1275                   1022                    936 
##         Mature_B_cells            Pre_B_cells         CD16-_NK_cells 
##                    780                    302                    187 
## CD34+CD38+CD123-_HSPCs         CD16+_NK_cells                   pDCs 
##                    166                    109                     71 
##              Basophils       CD34+CD38lo_HSCs            Pro_B_cells 
##                     59                     34                     25 
##         Plasma_B_cells CD34+CD38+CD123+_HSPCs 
##                     20                     14

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.

set.seed(42)
sce <- bHIVE(sce, assay.type = "exprs", targetK = 14,
             affinityParams = list(alpha = "auto"),
             maxIter = 20, verbose = FALSE)
table(sce$bHIVE)

## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14 
##  332   57 2211  197  242  391  134   83  296  123  143  214  183  394

ari_42 <- pairwiseRand(sce$population_id, sce$bHIVE, mode = "index")
ari_42

## [1] 0.5428569

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 0.54. 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.

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

Reading the rows shows what a run separates and what it merges. A population on more than one row was split. This run splits Mature_B_cells, Monocytes. A population with no row of its own was merged into its neighbors. Here that is CD16+_NK_cells, CD34+CD38+CD123+_HSPCs, CD34+CD38lo_HSCs, CD8_T_cells, Plasma_B_cells, Pro_B_cells. 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.

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 bluster. It works anywhere a BlusterParam does.

param <- bHIVEParam(targetK = 14, affinityParams = list(alpha = "auto"),
                    maxIter = 20)
param

## class: bHIVEParam
## nAntibodies: [default]
## beta: [default]
## epsilon: [default]
## maxIter: 20
## k: [default]
## affinityFunc: [default]
## distFunc: [default]
## initMethod: [default]
## scale: [default]
## targetK: 14
## epsilonQuantile: [default]
## affinityParams(1): alpha
## extra(0):

set.seed(42)
clusters <- clusterRows(X, param)
table(clusters)

## clusters
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14 
##  332   57 2211  197  242  391  134   83  296  123  143  214  183  394

The rows of X are the cells of sce, so the labels go straight into colData().

colLabels(sce) <- clusters
nlevels(colLabels(sce))

## [1] 14

bluster::clusterSweep() varies any setting of the parameter object. Here it tries three cluster counts.

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

##  targetK.8 targetK.14 targetK.20 
##  0.8913295  0.9113102  0.7927102

Comparison with FlowSOM and bluster methods

We compare bHIVE with three methods used for cytometry and single-cell data:

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.

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

method	clusters	ARI	seconds
bHIVE, seed 1	14	0.805	1.7
bHIVE, seed 2	14	0.869	1.7
bHIVE, seed 3	14	0.967	1.7
FlowSOM	14	0.897	5.1
k-means	14	0.583	0.0
SNN graph + Louvain	8	0.971	0.6

The graph-based method reaches an ARI of 0.97 with only 8 clusters. ARI is dominated by the large populations, so merging a few small ones costs little. FlowSOM reaches 0.90 and k-means 0.58. Across the three seeds here and the run above, bHIVE ranges from 0.54 to 0.97. 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.

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.

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)

	accuracy	balanced_accuracy
bHIVE	0.952	0.884
nearest centroid	0.973	0.944

Balanced accuracy averages the recall of each population, so rare populations count as much as common ones.

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

The nearest-centroid baseline beats bHIVE here. Its accuracy is 0.97 against 0.95, and its balanced accuracy 0.94 against 0.88. It matches or beats bHIVE on 12 of 14 populations. bHIVE’s per-class seeding keeps every population represented, and its lowest recall is 0.29 for CD34+CD38+CD123-_HSPCs.

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.

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

## 
##    CD4_T_cells      Monocytes    CD8_T_cells Mature_B_cells    Pre_B_cells 
##            514            458            338            258            119 
## CD16-_NK_cells 
##             95

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.

net <- AINet$new(nAntibodies = 40, maxIter = 15, targetK = 14,
                 coverageBoost = TRUE,
                 affinityParams = list(alpha = "auto"), verbose = FALSE)
set.seed(5)
net$fit(X)
net

## <AINet>
##   Repertoire: 52 antibodies x 32 features
##   Task: clustering
##   Iterations: 15

pairwiseRand(truth, net$result$assignments, mode = "index")

## [1] 0.604857

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:

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 caret, for cross-validated tuning of bHIVE classifiers.

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

## [1] 28 14

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

##   targetK metric_value
## 1      10     1.570817
## 2      14     1.665611
## 3      20     1.851999

Session information

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggplot2_4.0.3               bluster_1.23.1             
##  [3] HDCytoData_1.33.1           flowCore_2.25.1            
##  [5] ExperimentHub_3.3.2         AnnotationHub_4.3.2        
##  [7] BiocFileCache_3.3.0         dbplyr_2.6.0               
##  [9] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
## [11] Biobase_2.73.2              GenomicRanges_1.65.4       
## [13] Seqinfo_1.3.2               IRanges_2.47.5             
## [15] S4Vectors_0.51.10           BiocGenerics_0.59.12       
## [17] generics_0.1.4              MatrixGenerics_1.25.0      
## [19] matrixStats_1.5.0           bHIVE_0.99.8               
## [21] BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.3.0                   httr2_1.3.0                
##  [3] rlang_1.3.0                 magrittr_2.0.5             
##  [5] otel_0.2.0                  compiler_4.6.1             
##  [7] RSQLite_3.53.3              png_0.1-9                  
##  [9] vctrs_0.7.3                 pkgconfig_2.0.3            
## [11] crayon_1.5.3                fastmap_1.2.0              
## [13] backports_1.5.1             XVector_0.53.0             
## [15] labeling_0.4.3              rmarkdown_2.32             
## [17] purrr_1.2.2                 bit_4.6.0                  
## [19] xfun_0.61                   cachem_1.1.0               
## [21] jsonlite_2.0.0              blob_1.3.0                 
## [23] DelayedArray_0.39.8         tweenr_2.0.3               
## [25] BiocParallel_1.47.0         broom_1.0.13               
## [27] clusterCrit_1.3.0           parallel_4.6.1             
## [29] cluster_2.1.8.3             R6_2.6.1                   
## [31] bslib_0.12.0                RColorBrewer_1.1-3         
## [33] car_3.1-5                   jquerylib_0.1.4            
## [35] Rcpp_1.1.2                  knitr_1.52                 
## [37] BiocBaseUtils_1.15.1        Matrix_1.7-6               
## [39] igraph_2.3.4                tidyselect_1.2.1           
## [41] abind_1.4-8                 yaml_2.3.12                
## [43] codetools_0.2-20            curl_8.1.0                 
## [45] lattice_0.23-1              tibble_3.3.1               
## [47] withr_3.0.3                 KEGGREST_1.53.6            
## [49] S7_0.2.2                    evaluate_1.0.5             
## [51] Rtsne_0.17                  polyclip_1.10-7            
## [53] ConsensusClusterPlus_1.77.0 Biostrings_2.81.9          
## [55] ggpubr_1.0.0                pillar_1.11.1              
## [57] BiocManager_1.30.27         filelock_1.0.3             
## [59] carData_3.0-6               colorRamps_2.3.4           
## [61] BiocVersion_3.24.0          scales_1.4.0               
## [63] glue_1.8.1                  maketools_1.3.2            
## [65] tools_4.6.1                 ggnewscale_0.5.2           
## [67] BiocNeighbors_2.7.3         sys_3.4.3                  
## [69] ggsignif_0.6.4              buildtools_1.0.0           
## [71] XML_3.99-0.25               grid_4.6.1                 
## [73] tidyr_1.3.2                 RProtoBufLib_2.25.0        
## [75] AnnotationDbi_1.75.2        ggforce_0.5.0              
## [77] Formula_1.2-6               cli_3.6.6                  
## [79] rappdirs_0.3.4              cytolib_2.25.0             
## [81] S4Arrays_1.13.2             viridisLite_0.4.3          
## [83] dplyr_1.2.1                 gtable_0.3.6               
## [85] rstatix_1.1.0               sass_0.4.10                
## [87] digest_0.6.39               FlowSOM_2.21.0             
## [89] SparseArray_1.13.4          farver_2.1.2               
## [91] memoise_2.0.1               htmltools_0.5.9            
## [93] lifecycle_1.0.5             httr_1.4.9                 
## [95] MASS_7.3-66                 bit64_4.8.6

Immune-inspired prototypes for cytometry and single-cell data with bHIVE