Clustering real data: an OSCA workflow with clusterLibreville
Overview
This vignette follows the basic workflow of Orchestrating
Single-Cell Analysis with Bioconductor (OSCA) on a real dataset, and
clusters it with clusterLibreville: quality control,
normalisation, feature selection, PCA and UMAP as in OSCA, then
deterministic Leiden clustering, with a look at what determinism does
and does not buy.
Each step below names its package: quality control and normalisation from scuttle, pre-clustering for size factors, variance modelling and PCA from scran, UMAP and plots from scater, graphs from bluster.
Data: mouse brain cells (Zeisel et al. 2015)
The Zeisel brain data (3005 cells), one of the datasets in scRNAseq, comes with the authors’ cell type labels, so the clusters can be checked against them.
sce <- scRNAseq::ZeiselBrainData()
sce
#> class: SingleCellExperiment
#> dim: 20006 3005
#> metadata(0):
#> assays(1): counts
#> rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
#> rowData names(1): featureType
#> colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12
#> 1772058148_F03
#> colData names(9): tissue group # ... level1class level2class
#> reducedDimNames(0):
#> mainExpName: gene
#> altExpNames(2): repeat ERCC
table(sce$level1class)
#>
#> astrocytes_ependymal endothelial-mural interneurons
#> 224 235 290
#> microglia oligodendrocytes pyramidal CA1
#> 98 820 939
#> pyramidal SS
#> 399Preprocessing, as in OSCA
Quality control on library size, detected genes, mitochondrial and spike-in proportions:
is.mito <- grepl("^mt-", rownames(sce))
stats <- scuttle::perCellQCMetrics(sce, subsets = list(Mt = is.mito))
qc <- scuttle::quickPerCellQC(
stats, percent_subsets = c("altexps_ERCC_percent", "subsets_Mt_percent"))
sce <- sce[, !qc$discard]
ncol(sce)
#> [1] 2816Normalisation by deconvolution, variance modelling and the top 10% of
highly variable genes, then PCA and UMAP.
scran::quickCluster() is only a coarse pre-clustering, so
that scran::computeSumFactors() pools similar cells when
estimating size factors; it is not the clustering this vignette is
about:
set.seed(1000)
clusters <- scran::quickCluster(sce)
sce <- scran::computeSumFactors(sce, cluster = clusters)
sce <- scuttle::logNormCounts(sce)
dec <- scran::modelGeneVar(sce)
hvg <- scran::getTopHVGs(dec, prop = 0.1)
set.seed(1001)
sce <- scran::fixedPCA(sce, rank = 50, subset.row = hvg)
set.seed(1002)
sce <- scater::runUMAP(sce, dimred = "PCA")Clustering
On a SingleCellExperiment, leiden() builds
the shared nearest neighbour graph of a reduced dimension (as
bluster::makeSNNGraph() does) and returns the labels as a
factor, community 1 the largest:
colLabels(sce) <- leiden(sce, use.dimred = "PCA", k = 10, seed = 0)
table(colLabels(sce))
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14
#> 559 444 392 285 264 199 167 129 110 84 72 48 35 28
scater::plotReducedDim(sce, "UMAP", colour_by = "label", text_by = "label")
The clusters against the published cell types:
table(cluster = colLabels(sce), type = sce$level1class)
#> type
#> cluster astrocytes_ependymal endothelial-mural interneurons microglia
#> 1 9 0 0 0
#> 2 0 0 4 0
#> 3 5 2 3 2
#> 4 3 4 274 0
#> 5 0 0 2 0
#> 6 2 0 7 0
#> 7 10 151 0 3
#> 8 126 1 0 2
#> 9 0 1 0 0
#> 10 23 1 0 0
#> 11 1 0 0 70
#> 12 0 0 0 0
#> 13 0 0 0 0
#> 14 0 0 0 1
#> type
#> cluster oligodendrocytes pyramidal CA1 pyramidal SS
#> 1 550 0 0
#> 2 0 436 4
#> 3 1 22 357
#> 4 0 4 0
#> 5 0 261 1
#> 6 110 46 34
#> 7 2 0 1
#> 8 0 0 0
#> 9 0 109 0
#> 10 0 60 0
#> 11 1 0 0
#> 12 48 0 0
#> 13 35 0 0
#> 14 27 0 0
bluster::pairwiseRand(colLabels(sce), sce$level1class, mode = "index")
#> [1] 0.5525948At resolution 1 Leiden finds more clusters than there are published
types. Interneurons, astrocytes, endothelial-mural cells and microglia
each map to one cluster, while CA1 pyramidal cells and oligodendrocytes
are split into several subpopulations (the authors’ finer
level2class labels split them too), and one smaller cluster
mixes types. The adjusted Rand index counts every split against the
clustering, so it is low here although most clusters are pure; the
resolution scan below shows how it changes.
Determinism, and what it does not fix
The same seed gives the same clusters at any number of threads, so a result is reproducible on a laptop and on a many-core server:
Determinism makes each seed reproducible; it does not make the partition seed-independent. Modularity has many near-equal optima, so different seeds can still disagree on some cells:
runs <- lapply(0:4, function(s) leiden(sce, seed = s))
agree <- outer(seq_along(runs), seq_along(runs), Vectorize(function(a, b)
bluster::pairwiseRand(runs[[a]], runs[[b]], mode = "index")))
round(agree, 3)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.000 0.928 0.965 1.000 0.951
#> [2,] 0.928 1.000 0.913 0.928 0.906
#> [3,] 0.965 0.913 1.000 0.965 0.915
#> [4,] 1.000 0.928 0.965 1.000 0.951
#> [5,] 0.951 0.906 0.915 0.951 1.000Because every run is reproducible and cheap, many seeds (or resolutions) can be run and summarised, e.g. into a consensus or a per-cell stability score.
Did it converge?
details = TRUE returns how the run ended: the iterations
it took, whether the last one changed nothing (a fixed point), why it
stopped, and the modularity. By default Leiden runs until an iteration
changes nothing; with a fixed n_iterations it may stop
before:
str(leiden(sce, seed = 0, details = TRUE)[-1])
#> List of 4
#> $ iterations: int 5
#> $ converged : logi TRUE
#> $ stop : chr "converged"
#> $ modularity: num 0.78
str(leiden(sce, seed = 0, n_iterations = 1, details = TRUE)[-1])
#> List of 4
#> $ iterations: int 1
#> $ converged : logi FALSE
#> $ stop : chr "n_iterations"
#> $ modularity: num 0.777Resolution
Higher resolutions give more, smaller clusters:
res <- c(0.25, 0.5, 1, 2)
by_res <- lapply(res, function(r) leiden(sce, resolution = r, seed = 0))
data.frame(resolution = res,
clusters = vapply(by_res, nlevels, integer(1)),
ARI = vapply(by_res, function(l) bluster::pairwiseRand(
l, sce$level1class, mode = "index"), numeric(1)))
#> resolution clusters ARI
#> 1 0.25 5 0.7772167
#> 2 0.50 9 0.7150961
#> 3 1.00 14 0.5525948
#> 4 2.00 23 0.2889107plots <- lapply(c(1, 4), function(i) {
sce$res <- by_res[[i]]
r <- res[i]
scater::plotReducedDim(sce, "UMAP", colour_by = "res") +
ggplot2::ggtitle(paste("resolution", r))
})
patchwork::wrap_plots(plots)
Using clusterLibreville in bluster
cluster_libreville() has igraph’s community-function
interface, so it plugs into bluster’s graph-based clustering, for
example through bluster::clusterRows(); it gives the same
labels as leiden(sce) above:
Compared with the usual Bioconductor choices
scran::clusterCells() (OSCA’s entry point, now
bluster::clusterRows()) clusters the same SNN graph with
walktrap by default. bluster’s cluster.fun = "leiden"
passes no arguments to igraph::cluster_leiden(), so
igraph’s defaults apply: the CPM objective at resolution 1 (not
modularity) and 2 iterations. Each choice on the same graph:
# the SNN graph that leiden(sce) builds, clustered by each method
g <- bluster::makeSNNGraph(reducedDim(sce, "PCA"), k = 10, type = "rank")
choices <- list(
"walktrap (OSCA default)" = function(g) igraph::cluster_walktrap(g),
"igraph leiden, defaults (CPM)" = function(g) igraph::cluster_leiden(g),
"igraph leiden, modularity, 2 iterations" = function(g)
igraph::cluster_leiden(g, objective_function = "modularity"),
"clusterLibreville, 2 iterations" = function(g)
cluster_libreville(g, seed = 0, n_iterations = 2),
"igraph leiden, modularity, converged" = function(g)
igraph::cluster_leiden(g, objective_function = "modularity",
n_iterations = -1),
"clusterLibreville, converged" = function(g)
cluster_libreville(g, seed = 0))
run <- lapply(choices, function(f) {
# clustering only, on the SNN graph built above; best of 3
secs <- vapply(1:3, function(i) {
set.seed(0)
system.time(comm <<- f(g))[["elapsed"]]
}, numeric(1))
list(labels = factor(igraph::membership(comm)), seconds = min(secs))
})
labels <- lapply(run, `[[`, "labels")
ari <- function(a, b) bluster::pairwiseRand(a, b, mode = "index")
ref <- labels[["clusterLibreville, converged"]]
data.frame(clusters = vapply(labels, nlevels, integer(1)),
ARI_types = vapply(labels, ari, numeric(1), b = sce$level1class),
ARI_libreville = vapply(labels, ari, numeric(1), b = ref),
modularity = vapply(labels, function(l) modularity(g, l),
numeric(1)),
seconds = vapply(run, `[[`, numeric(1), "seconds"))
#> clusters ARI_types ARI_libreville
#> walktrap (OSCA default) 16 0.6320194 0.7939179
#> igraph leiden, defaults (CPM) 58 0.1482093 0.3082168
#> igraph leiden, modularity, 2 iterations 14 0.5519030 0.9989124
#> clusterLibreville, 2 iterations 14 0.5482494 0.9912963
#> igraph leiden, modularity, converged 14 0.5525948 1.0000000
#> clusterLibreville, converged 14 0.5525948 1.0000000
#> modularity seconds
#> walktrap (OSCA default) 0.7563095 1.011
#> igraph leiden, defaults (CPM) 0.6581287 0.027
#> igraph leiden, modularity, 2 iterations 0.7795432 0.029
#> clusterLibreville, 2 iterations 0.7795028 0.035
#> igraph leiden, modularity, converged 0.7796365 0.047
#> clusterLibreville, converged 0.7796365 0.050With the modularity objective, igraph’s Leiden and
clusterLibreville find essentially the same partition, at 2
iterations (igraph’s default) and run to convergence
(n_iterations = -1; igraph’s R package converges, but
python-igraph up to 0.11.9 runs a single iteration for -1).
Leiden through bluster with its default arguments optimises CPM
instead and splits the data into many small clusters: pass
objective_function = "modularity" (or a CPM resolution
suited to the graph’s weights) when using it. Walktrap is a different
algorithm, with its own trade-offs.
Does performance matter?
Not on a dataset this size. The seconds column above is
the clustering step alone, on the machine that built this vignette;
igraph runs on one thread and clusterLibreville on all
cores. At 2 iterations the two take about the same time, and run to
convergence igraph is the faster on 2816 cells:
clusterLibreville pays fixed costs (a thread pool, graph
conversion and validation) that a small graph doesn’t amortise.
The balance flips with size. On the 68,579 PBMCs of TENxPBMCData
("pbmc68k", the same preprocessing and the same kind of SNN
graph), the clustering step took:
| method | threads | seconds |
|---|---|---|
| igraph Leiden, 2 iterations | 1 | 49 |
| igraph Leiden, converged | 1 | 73 |
clusterLibreville, converged |
1 | 42 |
clusterLibreville, converged |
8 | 16 |
All runs converge to the same partition (10 clusters, adjusted Rand
index 1.000), and clusterLibreville’s labels are identical
at 1, 2, 4 and 8 threads. These times come from an 8-core laptop with
other jobs running, so the multi-threaded ones are pessimistic; they
were measured outside this vignette, which would otherwise take minutes
to build. The crossover lies somewhere between the two dataset sizes. So
clusterLibreville pays off in speed once datasets reach
tens of thousands of cells; its other advantage, reproducibility at any
number of threads, holds at every size.
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] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#> [3] Biobase_2.73.2 GenomicRanges_1.65.4
#> [5] Seqinfo_1.3.2 IRanges_2.47.5
#> [7] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [9] generics_0.1.4 MatrixGenerics_1.25.0
#> [11] matrixStats_1.5.0 clusterLibreville_0.99.0
#> [13] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_2.0.0
#> [4] magrittr_2.0.5 ggbeeswarm_0.7.3 GenomicFeatures_1.65.0
#> [7] gypsum_1.9.0 farver_2.1.2 rmarkdown_2.32
#> [10] BiocIO_1.23.3 vctrs_0.7.3 memoise_2.0.1
#> [13] Rsamtools_2.29.0 RCurl_1.98-1.20 htmltools_0.5.9
#> [16] S4Arrays_1.13.2 BiocBaseUtils_1.15.1 AnnotationHub_4.3.2
#> [19] curl_8.1.0 BiocNeighbors_2.7.3 Rhdf5lib_2.1.0
#> [22] SparseArray_1.13.4 rhdf5_2.57.20 sass_0.4.10
#> [25] alabaster.base_1.13.4 bslib_0.12.0 alabaster.sce_1.13.1
#> [28] httr2_1.3.0 cachem_1.1.0 buildtools_1.0.0
#> [31] GenomicAlignments_1.49.2 igraph_2.3.4 lifecycle_1.0.5
#> [34] pkgconfig_2.0.3 rsvd_1.0.6 Matrix_1.7-6
#> [37] R6_2.6.1 fastmap_1.2.0 digest_0.6.39
#> [40] patchwork_1.3.2 AnnotationDbi_1.75.2 scRNAseq_2.27.0
#> [43] scater_1.41.2 RSpectra_0.16-2 dqrng_0.4.1
#> [46] irlba_2.4.1 ExperimentHub_3.3.2 RSQLite_3.53.3
#> [49] beachmat_2.29.3 labeling_0.4.3 filelock_1.0.3
#> [52] httr_1.4.9 abind_1.4-8 compiler_4.6.1
#> [55] withr_3.0.3 bit64_4.8.6 S7_0.2.2
#> [58] BiocParallel_1.47.0 viridis_0.6.5 DBI_1.3.0
#> [61] HDF5Array_1.41.3 alabaster.ranges_1.13.1 alabaster.schemas_1.13.0
#> [64] rappdirs_0.3.4 DelayedArray_0.39.8 rjson_0.2.23
#> [67] bluster_1.23.1 tools_4.6.1 vipor_0.4.7
#> [70] otel_0.2.0 beeswarm_0.4.0 glue_1.8.1
#> [73] h5mread_1.5.3 restfulr_0.0.17 rhdf5filters_1.25.4
#> [76] grid_4.6.1 cluster_2.1.8.3 gtable_0.3.6
#> [79] ensembldb_2.37.3 BiocSingular_1.29.2 ScaledMatrix_1.21.0
#> [82] metapod_1.21.0 XVector_0.53.0 ggrepel_0.9.8
#> [85] BiocVersion_3.24.0 pillar_1.11.1 limma_3.99.2
#> [88] dplyr_1.2.1 BiocFileCache_3.3.0 lattice_0.23-1
#> [91] FNN_1.1.4.1 rtracklayer_1.73.0 bit_4.6.0
#> [94] tidyselect_1.2.1 locfit_1.5-9.12 maketools_1.3.2
#> [97] Biostrings_2.81.9 scuttle_1.23.2 knitr_1.52
#> [100] gridExtra_2.3.1 ProtGenerics_1.45.0 edgeR_4.99.10
#> [103] xfun_0.62 statmod_1.5.2 UCSC.utils_1.9.0
#> [106] lazyeval_0.2.3 yaml_2.3.12 evaluate_1.0.5
#> [109] codetools_0.2-20 cigarillo_1.3.1 tibble_3.3.1
#> [112] alabaster.matrix_1.13.4 BiocManager_1.30.27 cli_3.6.6
#> [115] uwot_0.2.5 jquerylib_0.1.4 Rcpp_1.1.2
#> [118] GenomeInfoDb_1.49.1 dbplyr_2.6.0 png_0.1-9
#> [121] XML_3.99-0.25 parallel_4.6.1 ggplot2_4.0.3
#> [124] blob_1.3.0 AnnotationFilter_1.37.0 scran_1.41.1
#> [127] sparseMatrixStats_1.25.0 bitops_1.1-0 viridisLite_0.4.3
#> [130] alabaster.se_1.13.1 scales_1.4.0 crayon_1.5.3
#> [133] rlang_1.3.0 cowplot_1.2.0 KEGGREST_1.53.6