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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## ----modules-vignette, eval = FALSE-------------------------------------------
# vignette("immune-modules", package = "bHIVE")

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