Flexible clustering for Bioconductor
Introduction
The bluster
package provides a flexible and extensible framework for clustering in
Bioconductor packages/workflows. At its core is the
clusterRows() generic that controls dispatch to different
clustering algorithms. We will demonstrate on some single-cell RNA
sequencing data from the scRNAseq
package; our aim is to cluster cells into cell populations based on
their PC coordinates.
library(scRNAseq)
sce <- ZeiselBrainData()
# Trusting the authors' quality control, and going straight to normalization.
library(scuttle)
sce <- logNormCounts(sce)
# Feature selection based on highly variable genes.
library(scran)
dec <- modelGeneVar(sce)
hvgs <- getTopHVGs(dec, n=1000)
# Dimensionality reduction for work (PCA) and pleasure (t-SNE).
set.seed(1000)
library(scater)
sce <- runPCA(sce, ncomponents=20, subset_row=hvgs)
sce <- runUMAP(sce, dimred="PCA")
mat <- reducedDim(sce, "PCA")
dim(mat)## [1] 3005 20Based on distance matrices
Hierarchical clustering
Our first algorithm is good old hierarchical clustering, as
implemented using hclust() from the stats package.
This automatically sets the cut height to half the dendrogram
height.
library(bluster)
hclust.out <- clusterRows(mat, HclustParam())
plotUMAP(sce, colour_by=I(hclust.out))
Advanced users can achieve greater control of the procedure by
passing more parameters to the HclustParam() constructor.
Here, we use Ward’s criterion for the agglomeration with a dynamic tree
cut from the dynamicTreeCut
package.
## class: HclustParam
## metric: [default]
## method: ward.D2
## cut.fun: cutreeDynamic
## cut.params(0):
## dist.fun: stats::dist
## clust.fun: stats::hclust
Affinity propagation
Another option is to use affinity propagation, as implemented using the apcluster package. Here, messages are passed between observations to decide on a set of exemplars, each of which form the center of a cluster.
This is not particularly fast as it involves the calculation of a square similarity matrix between all pairs of observations. So, we’ll speed it up by taking analyzing a subset of the data:
set.seed(1000)
sub <- sce[,sample(ncol(sce), 200)]
ap.out <- clusterRows(reducedDim(sub), AffinityParam())
plotUMAP(sub, colour_by=I(ap.out))
The q parameter is probably the most important and
determines the resolution of the clustering. This can be set to any
value below 1, with smaller (possibly negative) values corresponding to
coarser clusters.
set.seed(1000)
ap.out <- clusterRows(reducedDim(sub), AffinityParam(q=-2))
plotUMAP(sub, colour_by=I(ap.out))
With a fixed number of clusters
k-means clustering
A classic algorithm is k-means clustering, as implemented
using the kmeans() function. This requires us to pass in
the number of clusters, either as a number:
set.seed(100)
kmeans.out <- clusterRows(mat, KmeansParam(10))
plotUMAP(sce, colour_by=I(kmeans.out))
Or as a function of the number of observations, which is useful for vector quantization purposes:
## class: KmeansParam
## centers: variable
## iter.max: 10
## nstart: 1
## algorithm: Hartigan-Wong
A variant of this approach is mini-batch k-means, as implemented in the mbkmeans package. This uses mini-batching to approximate the full k-means algorithm for greater speed.
set.seed(100)
mbkmeans.out <- clusterRows(mat, MbkmeansParam(20))
plotUMAP(sce, colour_by=I(mbkmeans.out))
Self-organizing maps
We can also use self-organizing maps (SOMs) from the kohonen package. This allocates observations to nodes of a simple neural network; each node is then treated as a cluster.
The key feature of SOMs is that they apply additional topological constraints on the relative positions of the nodes. This allows us to naturally determine the relationships between clusters based on the proximity of the nodes.
set.seed(1000)
som.out <- clusterRows(mat, SomParam(100), full=TRUE)
par(mfrow=c(1,2))
plot(som.out$objects$som, "counts")
grid <- som.out$objects$som$grid$pts
text(grid[,1], grid[,2], seq_len(nrow(grid)))
Graph-based clustering
We can build shared or direct nearest neighbor graphs and perform
community detection with igraph. Here,
the number of neighbors k controls the resolution of the
clusters.
set.seed(101) # just in case there are ties.
graph.out <- clusterRows(mat, NNGraphParam(k=10))
plotUMAP(sce, colour_by=I(graph.out))
It is again straightforward to tune the procedure by passing more arguments such as the community detection algorithm to use.
## class: SNNGraphParam
## k: 20
## type: rank
## BNPARAM: KmknnParam
## num.threads: 1
## cluster.fun: louvain
## cluster.args(0):
Density-based clustering
We also implement a version of the DBSCAN algorithm for density-based clustering. This focuses on identifying masses of continuous density.
Unlike the other methods, it will consider certain points to be noise and discard them from the output.
## Mode FALSE TRUE
## logical 1789 1216The resolution of the clustering can be modified by tinkering with
the core.prop. Smaller values will generally result in
smaller clusters as well as more noise points.
## Mode FALSE TRUE
## logical 565 2440Two-phase clustering
We also provide a wrapper for a hybrid “two-step” approach for handling larger datasets. Here, a fast agglomeration is performed with k-means to compact the data, followed by a slower graph-based clustering step to obtain interpretable meta-clusters. (This dataset is not, by and large, big enough for this approach to work particularly well.)
set.seed(100) # for the k-means
two.out <- clusterRows(mat, TwoStepParam())
plotUMAP(sce, colour_by=I(two.out))
Each step is itself parametrized by BlusterParam
objects, so it is possible to tune them individually:
## class: TwoStepParam
## first:
## class: KmeansParam
## centers: variable
## iter.max: 10
## nstart: 1
## algorithm: Hartigan-Wong
## second:
## class: SNNGraphParam
## k: 5
## type: rank
## BNPARAM: KmknnParam
## num.threads: 1
## cluster.fun: walktrap
## cluster.args(0):set.seed(100) # for the k-means
two.out <- clusterRows(mat, TwoStepParam())
plotUMAP(sce, colour_by=I(two.out))
Obtaining full clustering statistics
Sometimes the vector of cluster assignments is not enough. We can
obtain more information about the clustering procedure by setting
full=TRUE in clusterRows(). For example, we
could obtain the actual graph generated by
NNGraphParam():
## IGRAPH 7a7b07b U-W- 3005 87028 --
## + attr: weight (e/n)
## + edges from 7a7b07b:
## [1] 1-- 2 1-- 3 2-- 3 1-- 4 2-- 4 3-- 4 4-- 5 1-- 5 2-- 5 3-- 5
## [11] 1-- 6 2-- 6 3-- 6 4-- 6 5-- 6 1-- 7 5-- 7 6-- 7 2-- 7 4-- 7
## [21] 3-- 7 3-- 8 1-- 8 2-- 8 4-- 8 6-- 8 7-- 8 5-- 8 1-- 9 2-- 9
## [31] 4-- 9 5-- 9 7-- 9 3-- 9 6-- 9 8-- 9 1--10 5--10 9--10 3--10
## [41] 4--10 3--11 9--11 10--11 1--11 5--11 2--13 5--13 6--13 7--13
## [51] 12--13 8--13 12--14 13--14 12--15 13--15 14--15 2--15 5--15 6--15
## [61] 7--15 12--16 13--16 14--16 15--16 14--17 12--17 15--17 16--17 13--17
## [71] 12--18 14--18 15--18 16--18 13--18 17--18 19--20 12--20 14--20 15--20
## + ... omitted several edgesThe assignment vector is still reported in the clusters
entry:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 166 216 256 543 218 89 583 92 123 186 91 42 94 64 58 58 60 35 31Further comments
clusterRows() enables users or developers to easily
switch between clustering algorithms by changing a single argument.
Indeed, by passing the BlusterParam object across
functions, we can ensure that the same algorithm is used through a
workflow. It is also helpful for package functions as it provides
diverse functionality without compromising a clean function signature.
However, the true power of this framework lies in its extensibility.
Anyone can write a clusterRows() method for a new
clustering algorithm with an associated BlusterParam
subclass, and that procedure is immediately compatible with any workflow
or function that was already using clusterRows().
Session information
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.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.26.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=C
## [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] bluster_1.17.0 scater_1.35.0
## [3] ggplot2_3.5.1 scran_1.35.0
## [5] scuttle_1.17.0 scRNAseq_2.20.0
## [7] SingleCellExperiment_1.29.1 SummarizedExperiment_1.37.0
## [9] Biobase_2.67.0 GenomicRanges_1.59.1
## [11] GenomeInfoDb_1.43.2 IRanges_2.41.2
## [13] S4Vectors_0.45.2 BiocGenerics_0.53.3
## [15] generics_0.1.3 MatrixGenerics_1.19.0
## [17] matrixStats_1.4.1 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] sys_3.4.3 jsonlite_1.8.9 magrittr_2.0.3
## [4] ggbeeswarm_0.7.2 GenomicFeatures_1.59.1 gypsum_1.3.0
## [7] farver_2.1.2 rmarkdown_2.29 BiocIO_1.17.1
## [10] zlibbioc_1.52.0 vctrs_0.6.5 memoise_2.0.1
## [13] Rsamtools_2.23.1 RCurl_1.98-1.16 benchmarkme_1.0.8
## [16] htmltools_0.5.8.1 S4Arrays_1.7.1 dynamicTreeCut_1.63-1
## [19] AnnotationHub_3.15.0 curl_6.0.1 BiocNeighbors_2.1.2
## [22] Rhdf5lib_1.29.0 SparseArray_1.7.2 rhdf5_2.51.1
## [25] sass_0.4.9 alabaster.base_1.7.2 bslib_0.8.0
## [28] alabaster.sce_1.7.0 httr2_1.0.7 cachem_1.1.0
## [31] buildtools_1.0.0 GenomicAlignments_1.43.0 igraph_2.1.2
## [34] iterators_1.0.14 lifecycle_1.0.4 pkgconfig_2.0.3
## [37] rsvd_1.0.5 Matrix_1.7-1 R6_2.5.1
## [40] fastmap_1.2.0 GenomeInfoDbData_1.2.13 digest_0.6.37
## [43] colorspace_2.1-1 AnnotationDbi_1.69.0 dqrng_0.4.1
## [46] irlba_2.3.5.1 ExperimentHub_2.15.0 RSQLite_2.3.9
## [49] beachmat_2.23.4 labeling_0.4.3 filelock_1.0.3
## [52] fansi_1.0.6 httr_1.4.7 abind_1.4-8
## [55] compiler_4.4.2 doParallel_1.0.17 withr_3.0.2
## [58] bit64_4.5.2 BiocParallel_1.41.0 viridis_0.6.5
## [61] DBI_1.2.3 apcluster_1.4.13 HDF5Array_1.35.2
## [64] alabaster.ranges_1.7.0 alabaster.schemas_1.7.0 rappdirs_0.3.3
## [67] DelayedArray_0.33.3 rjson_0.2.23 tools_4.4.2
## [70] vipor_0.4.7 mbkmeans_1.23.0 beeswarm_0.4.0
## [73] glue_1.8.0 restfulr_0.0.15 rhdf5filters_1.19.0
## [76] grid_4.4.2 cluster_2.1.7 gtable_0.3.6
## [79] ensembldb_2.31.0 BiocSingular_1.23.0 ScaledMatrix_1.15.0
## [82] metapod_1.15.0 utf8_1.2.4 XVector_0.47.0
## [85] foreach_1.5.2 ggrepel_0.9.6 BiocVersion_3.21.1
## [88] pillar_1.9.0 limma_3.63.2 benchmarkmeData_1.0.4
## [91] dplyr_1.1.4 BiocFileCache_2.15.0 lattice_0.22-6
## [94] gmp_0.7-5 FNN_1.1.4.1 rtracklayer_1.67.0
## [97] bit_4.5.0.1 tidyselect_1.2.1 locfit_1.5-9.10
## [100] maketools_1.3.1 Biostrings_2.75.2 knitr_1.49
## [103] gridExtra_2.3 ProtGenerics_1.39.0 edgeR_4.5.1
## [106] xfun_0.49 statmod_1.5.0 UCSC.utils_1.3.0
## [109] lazyeval_0.2.2 yaml_2.3.10 evaluate_1.0.1
## [112] codetools_0.2-20 tibble_3.2.1 alabaster.matrix_1.7.4
## [115] BiocManager_1.30.25 cli_3.6.3 uwot_0.2.2
## [118] munsell_0.5.1 jquerylib_0.1.4 Rcpp_1.0.13-1
## [121] dbplyr_2.5.0 png_0.1-8 kohonen_3.0.12
## [124] XML_3.99-0.17 parallel_4.4.2 blob_1.2.4
## [127] ClusterR_1.3.3 AnnotationFilter_1.31.0 bitops_1.0-9
## [130] viridisLite_0.4.2 alabaster.se_1.7.0 scales_1.3.0
## [133] crayon_1.5.3 rlang_1.1.4 KEGGREST_1.47.0