Contents

1 Introduction

The bluster package provides a few diagnostics for quantitatively examining the cluster output. These diagnostics We will demonstrate on another dataset from the scRNAseq package, clustered with graph-based methods via the clusterRows() generic as described in the previous vignette.

library(scRNAseq)
sce <- GrunPancreasData()

# Quality control to remove bad cells.
library(scuttle)
qcstats <- perCellQCMetrics(sce)
qcfilter <- quickPerCellQC(qcstats, sub.fields="altexps_ERCC_percent")
sce <- sce[,!qcfilter$discard]

# Normalization by library size.
sce <- logNormCounts(sce)

# Feature selection.
library(scran)
dec <- modelGeneVar(sce)
hvgs <- getTopHVGs(dec, n=1000)

# Dimensionality reduction.
set.seed(1000)
library(scater)
sce <- runPCA(sce, ncomponents=20, subset_row=hvgs)

# Clustering.
library(bluster)
mat <- reducedDim(sce)
clust.info <- clusterRows(mat, NNGraphParam(), full=TRUE)
clusters <- clust.info$clusters
table(clusters)
## clusters
##   1   2   3   4   5   6   7   8   9  10  11  12 
## 285 171 161  59 174  49  70 137  69  65  28  23

2 Computing the silhouette width

The silhouette width is a standard metric to quantify the separation between clusters generated by any procedure. A cell with a large positive width is closer to other cells from the same cluster compared to cells from different clusters. On the other hand, low or negative widths indicate that cells from different clusters are not well separated.

The exact silhouette calculation is rather computationally intensive so bluster implements an approximation instead. This is provided in the approxSilhouette() function, which returns the width for each cell and its closest (non-self) cluster. Clusters consisting of cells with lower widths may warrant some more care during interpretation.

sil <- approxSilhouette(mat, clusters)
sil
## DataFrame with 1291 rows and 3 columns
##             cluster    other     width
##            <factor> <factor> <numeric>
## D2ex_1            1        7  0.375510
## D2ex_2            1        7  0.370596
## D2ex_3            1        7  0.409271
## D2ex_4            5        3  0.257069
## D2ex_5            5        3  0.249883
## ...             ...      ...       ...
## D17TGFB_91        6        2  0.126607
## D17TGFB_92        1        5  0.403106
## D17TGFB_93        1        7  0.384365
## D17TGFB_94        6        2  0.243933
## D17TGFB_95        6        2  0.204860
boxplot(split(sil$width, clusters))

The function also returns the identity of the closest “other” cluster for each cell. This can be helpful to identify which clusters are easily confused to each other, based on how many of one cluster’s cells are closer to the other cluster.

best.choice <- ifelse(sil$width > 0, clusters, sil$other)
table(Assigned=clusters, Closest=best.choice)
##         Closest
## Assigned   1   2   3   4   5   6   7   8   9  10  11  12
##       1  275   0   7   0   3   0   0   0   0   0   0   0
##       2    0 171   0   0   0   0   0   0   0   0   0   0
##       3    0   0 150   0  11   0   0   0   0   0   0   0
##       4    0   0   0  59   0   0   0   0   0   0   0   0
##       5    0   0  21   0 153   0   0   0   0   0   0   0
##       6    0  26   0   0   0  23   0   0   0   0   0   0
##       7   14   0   8  10   2   0  18   0   7   2   8   1
##       8    0   0   0   1   0   0   0 136   0   0   0   0
##       9    0   0   0   0   0   0   0  10  58   0   0   1
##       10   0   0   0   0   0   0   0   0   0  65   0   0
##       11   0   0   0   0   0   0   0   0   0   0  28   0
##       12   0   0   0   0   0   0   0   0   0   0   0  23

3 Computing the neighborhood purity

Another diagnostic uses the percentage of neighbors for each cell that belong to the same cluster. Well-separated clusters should exhibit high percentages (i.e., “purities”) as cells from different clusters do not mix. Low purities are symptomatic of overclustering where cluster boundaries become more ambiguous.

The neighborPurity() function computes the purity of the neighborhood for each cell. Clusters with systematically low purities may warrant some more care during interpretation. By default, we perform some weighting so that large clusters do not have large purities simply because there are few cells assigned to other clusters in the dataset.

pure <- neighborPurity(mat, clusters)
pure
## DataFrame with 1291 rows and 2 columns
##               purity  maximum
##            <numeric> <factor>
## D2ex_1      1.000000        1
## D2ex_2      1.000000        1
## D2ex_3      1.000000        1
## D2ex_4      0.951603        5
## D2ex_5      1.000000        5
## ...              ...      ...
## D17TGFB_91  0.959277        6
## D17TGFB_92  1.000000        1
## D17TGFB_93  1.000000        1
## D17TGFB_94  0.986538        6
## D17TGFB_95  0.898817        6
boxplot(split(pure$purity, clusters))

The function also returns the identity of the other cluster with the highest percentage. This can again be useful to identify the relationships between clusters based on which pairs have the greatest intermingling in their neighborhoods.

table(Assigned=clusters, Max=pure$maximum)
##         Max
## Assigned   1   2   3   4   5   6   7   8   9  10  11  12
##       1  279   0   0   0   0   0   6   0   0   0   0   0
##       2    0 170   0   0   0   0   0   0   0   0   0   1
##       3    0   0 158   0   3   0   0   0   0   0   0   0
##       4    0   0   0  59   0   0   0   0   0   0   0   0
##       5    1   0   6   0 165   0   2   0   0   0   0   0
##       6    0  18   0   0   0  31   0   0   0   0   0   0
##       7    0   0   0   1   0   0  67   0   1   0   1   0
##       8    0   0   0   0   0   0   0 136   1   0   0   0
##       9    0   0   0   0   0   0   0   5  64   0   0   0
##       10   0   0   0   0   0   0   0   0   0  65   0   0
##       11   0   0   0   0   0   0   0   0   0   0  28   0
##       12   0   0   0   0   0   0   0   0   0   0   0  23

4 Computing per-cluster RMSD

The root-mean-squared deviation (RMSD) for each cluster represents the dispersion of cells within each cluster. A large RMSD value indicates that a cluster has high internal heterogeneity, making it a good candidate for further subclustering.

rmsd <- clusterRMSD(mat, clusters)
barplot(rmsd)

Alternatively, we can compute the within-cluster sum of squares (WCSS), a metric commonly seen in \(k\)-means clustering. One could pick a “sensible” choice for \(k\) by computing the WCSS for a range of values and picking the value the WCSS begins to plateau - see Section 8 for more details.

clusterRMSD(mat, clusters, sum=TRUE)
##         1         2         3         4         5         6         7         8 
## 60368.794  9537.387 26566.200 19009.475 30701.306  8244.534 31829.182 14682.089 
##         9        10        11        12 
## 10652.096  6094.771  2644.843  1956.403

5 Computing graph modularity

For graph-based methods, we can compute the cluster modularity within clusters and between pairs of clusters. Specifically, we examine the ratio of observed to expected edge weights for each pair of clusters (closely related to the modularity score used in many cluster_* functions from igraph). We would usually expect to see high observed weights between cells in the same cluster with minimal weights between clusters, indicating that the clusters are well-separated. Large off-diagonal entries indicate that the corresponding pair of clusters are closely related.

g <- clust.info$objects$graph
ratio <- pairwiseModularity(g, clusters, as.ratio=TRUE)

library(pheatmap)
pheatmap(log10(ratio+1), cluster_cols=FALSE, cluster_rows=FALSE,
    col=rev(heat.colors(100)))

This may be better visualized with a force-directed layout:

cluster.gr <- igraph::graph_from_adjacency_matrix(log2(ratio+1), 
    mode="upper", weighted=TRUE, diag=FALSE)

# Increasing the weight to increase the visibility of the lines.
set.seed(1100101)
plot(cluster.gr, edge.width=igraph::E(cluster.gr)$weight*5,
    layout=igraph::layout_with_lgl)

We can also tune the resolution of the clustering post hoc with the mergeCommunities() function. This will iteratively merge the most closely related pair of clusters together until the desired number of clusters is reached. For example, if we wanted to whittle down the number of clusters to 10, we could do:

merged <- mergeCommunities(g, clusters, number=10)
table(merged)
## merged
##   1   2   3   5   6   7   9  10  11  12 
## 285 171 161 174  49 129 206  65  28  23

6 Comparing two clusterings

To compare two clusterings, the pairwiseRand() function computes the adjusted Rand index (ARI). High ARIs indicate that the two clusterings are similar with respect to how they partition the observations, and an ARI of 1 means that the clusterings are identical.

hclusters <- clusterRows(mat, HclustParam(cut.dynamic=TRUE))
pairwiseRand(clusters, hclusters, mode="index")
## [1] 0.4512108

Of course, a single number is not particularly useful, so clusterRand() also provides the capability to break down the ARI into its contributions from each cluster or cluster pair. Specifically, for each cluster or cluster pair in a “reference” clustering (here, clusters), we see whether it is preserved in the “alternative” clustering (here, hclusters). Large values on the diagonal indicate that the reference cluster is recapitulated; large values off the diagonal indicate that the separation between the corresponding pair of clusters is also maintained. Conversely, low diagonal values indicate that the corresponding cluster is fragmented in the alternative, and low off-diagonal values can be used as a diagnostic for loss of separation.

ratio <- pairwiseRand(clusters, hclusters, mode="ratio")

library(pheatmap)
pheatmap(ratio, cluster_cols=FALSE, cluster_rows=FALSE,
    col=viridis::viridis(100), breaks=seq(-1, 1, length=101))

Explicit mappings between two clusterings can be performed using linkClusters() (see Section 9). Alternatively, we can quantify the degree of “nesting” of one clustering within another with nestedClusters(); this can be useful for verifying that the higher-resolution clustering is indeed nested within its coarser counterpart.

7 Bootstrapping cluster stability

We can use bootstrapping to evaluate the effect of sampling noise on the stability of a clustering procedure. The bootstrapStability() function will return the ARI of the original clusters against those generated from bootstrap replicates, averaged across multiple bootstrap iterations. High values indicate that the clustering is robust to sample noise.

set.seed(1001010)
ari <-bootstrapStability(mat, clusters=clusters, 
    mode="index", BLUSPARAM=NNGraphParam())
ari
## [1] 0.6824535

Advanced users may also set mode="ratio" to obtain a more detailed breakdown of the effect of noise on each cluster (pair).

set.seed(1001010)
ratio <-bootstrapStability(mat, clusters=clusters,
    mode="ratio", BLUSPARAM=NNGraphParam())

library(pheatmap)
pheatmap(ratio, cluster_cols=FALSE, cluster_rows=FALSE,
    col=viridis::viridis(100), breaks=seq(-1, 1, length=101))

8 Clustering parameter sweeps

The clusterSweep() function provides a convenient way to test multiple combinations of parameter settings. Given a BlusterParam object and a set of values for each parameter, the function will repeat the clustering ith each combination of parameters. The example below uses graph-based clustering with a variety of k as well as different community detection algorithms.

combinations <- clusterSweep(mat, BLUSPARAM=SNNGraphParam(),
    k=c(5L, 10L, 15L, 20L), cluster.fun=c("walktrap", "louvain", "infomap"))

This yields a list containing all clusterings and the corresponding parameter combinations used to generate them. The function will attempt to generate some sensible name for each combination, though this may require some manual curation for large numbers of parameters.

colnames(combinations$clusters)
##  [1] "k.5_cluster.fun.walktrap"  "k.10_cluster.fun.walktrap"
##  [3] "k.15_cluster.fun.walktrap" "k.20_cluster.fun.walktrap"
##  [5] "k.5_cluster.fun.louvain"   "k.10_cluster.fun.louvain" 
##  [7] "k.15_cluster.fun.louvain"  "k.20_cluster.fun.louvain" 
##  [9] "k.5_cluster.fun.infomap"   "k.10_cluster.fun.infomap" 
## [11] "k.15_cluster.fun.infomap"  "k.20_cluster.fun.infomap"
combinations$parameters
## DataFrame with 12 rows and 2 columns
##                                   k cluster.fun
##                           <integer> <character>
## k.5_cluster.fun.walktrap          5    walktrap
## k.10_cluster.fun.walktrap        10    walktrap
## k.15_cluster.fun.walktrap        15    walktrap
## k.20_cluster.fun.walktrap        20    walktrap
## k.5_cluster.fun.louvain           5     louvain
## ...                             ...         ...
## k.20_cluster.fun.louvain         20     louvain
## k.5_cluster.fun.infomap           5     infomap
## k.10_cluster.fun.infomap         10     infomap
## k.15_cluster.fun.infomap         15     infomap
## k.20_cluster.fun.infomap         20     infomap

We can combine this with some of the metrics defined above to quantify cluster separation as a function of the clustering parameters. This allows us to quickly determine which parameters have a noticeable impact on the results. In principle, we could choose the clustering with the greatest separation for further analysis; however, this tends to be disappointing as it often favors overly broad clusters.

set.seed(10)
nclusters <- 3:25
kcombos <- clusterSweep(mat, BLUSPARAM=KmeansParam(centers=5), centers=nclusters)

sil <- vapply(as.list(kcombos$clusters), function(x) mean(approxSilhouette(mat, x)$width), 0)
plot(nclusters, sil, xlab="Number of clusters", ylab="Average silhouette width")

pur <- vapply(as.list(kcombos$clusters), function(x) mean(neighborPurity(mat, x)$purity), 0)
plot(nclusters, pur, xlab="Number of clusters", ylab="Average purity")

wcss <- vapply(as.list(kcombos$clusters), function(x) sum(clusterRMSD(mat, x, sum=TRUE)), 0)
plot(nclusters, wcss, xlab="Number of clusters", ylab="Within-cluster sum of squares")

9 Linking clusters

If we have many clusterings, we can identify corresponding clusters with the linkClusters() function. This constructs a graph where edges are formed between pairs of clusters from different clusterings, based on the number of cells assigned to both clusters. Re-using some of the clusterings from our previous sweep, we might do:

linked <- linkClusters(
    list(
        walktrap=combinations$clusters$k.10_cluster.fun.walktrap,
        louvain=combinations$clusters$k.10_cluster.fun.louvain,
        infomap=combinations$clusters$k.10_cluster.fun.infomap
    )
)
linked
## IGRAPH bc17ed4 UNW- 45 88 -- 
## + attr: name (v/c), weight (e/n)
## + edges from bc17ed4 (vertex names):
##  [1] walktrap.1 --louvain.1  walktrap.7 --louvain.1  walktrap.1 --louvain.2 
##  [4] walktrap.3 --louvain.2  walktrap.5 --louvain.2  walktrap.3 --louvain.3 
##  [7] walktrap.7 --louvain.3  walktrap.11--louvain.3  walktrap.3 --louvain.4 
## [10] walktrap.4 --louvain.5  walktrap.7 --louvain.5  walktrap.6 --louvain.6 
## [13] walktrap.8 --louvain.6  walktrap.10--louvain.7  walktrap.2 --louvain.8 
## [16] walktrap.6 --louvain.8  walktrap.3 --louvain.9  walktrap.5 --louvain.9 
## [19] walktrap.8 --louvain.10 walktrap.9 --louvain.10 walktrap.12--louvain.11
## [22] walktrap.1 --louvain.12 walktrap.3 --louvain.12 walktrap.7 --louvain.12
## + ... omitted several edges

This can be used to visualize the relationships between clusters, or to identify metaclusters across clusterings with community detection algorithms:

meta <- igraph::cluster_walktrap(linked)
plot(linked, mark.groups=meta)

By default, the edge weights are computed by dividing the number of shared cells with the smaller of the total number of cells in either cluster. This favors strong edges between a large cluster in one clustering and smaller subcluster in another (finer) clustering. Alternative weighting schemes will favour a 1:1 mapping between clusterings, which can be easier to interpret.

10 Comparing multiple clusterings

The compareClusterings() function will return a symmetric matrix of the ARIs between pairs of different clusterings. This is helpful for visualizing the relationships between different clusterings, e.g., to see which parameters most contribute to differences between clusterings.

aris <- compareClusterings(combinations$clusters)
g <- igraph::graph.adjacency(aris, mode="undirected", weighted=TRUE)
meta2 <- igraph::cluster_walktrap(g)
plot(g, mark.groups=meta2)

We can also identify groups of clusterings, typically corresponding to parameter combinations that yield more-or-less similar results. This allows us to prune out combinations that are largely redundant prior to downstream analyses.

Session information

sessionInfo()
## R Under development (unstable) (2024-10-21 r87258)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              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: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12             bluster_1.17.0             
##  [3] scater_1.35.0               ggplot2_3.5.1              
##  [5] scran_1.35.0                scuttle_1.17.0             
##  [7] scRNAseq_2.19.1             SingleCellExperiment_1.29.0
##  [9] SummarizedExperiment_1.37.0 Biobase_2.67.0             
## [11] GenomicRanges_1.59.0        GenomeInfoDb_1.43.0        
## [13] IRanges_2.41.0              S4Vectors_0.45.0           
## [15] BiocGenerics_0.53.0         MatrixGenerics_1.19.0      
## [17] matrixStats_1.4.1           BiocStyle_2.35.0           
## 
## loaded via a namespace (and not attached):
##   [1] BiocIO_1.17.0            bitops_1.0-9             filelock_1.0.3          
##   [4] tibble_3.2.1             XML_3.99-0.17            lifecycle_1.0.4         
##   [7] httr2_1.0.5              edgeR_4.5.0              doParallel_1.0.17       
##  [10] lattice_0.22-6           ensembldb_2.31.0         alabaster.base_1.7.0    
##  [13] magrittr_2.0.3           limma_3.63.0             sass_0.4.9              
##  [16] rmarkdown_2.28           jquerylib_0.1.4          yaml_2.3.10             
##  [19] metapod_1.15.0           RColorBrewer_1.1-3       cowplot_1.1.3           
##  [22] DBI_1.2.3                abind_1.4-8              zlibbioc_1.53.0         
##  [25] AnnotationFilter_1.31.0  RCurl_1.98-1.16          rappdirs_0.3.3          
##  [28] GenomeInfoDbData_1.2.13  ggrepel_0.9.6            irlba_2.3.5.1           
##  [31] kohonen_3.0.12           alabaster.sce_1.7.0      dqrng_0.4.1             
##  [34] codetools_0.2-20         DelayedArray_0.33.0      tidyselect_1.2.1        
##  [37] UCSC.utils_1.3.0         farver_2.1.2             gmp_0.7-5               
##  [40] ScaledMatrix_1.15.0      viridis_0.6.5            BiocFileCache_2.15.0    
##  [43] dynamicTreeCut_1.63-1    GenomicAlignments_1.43.0 jsonlite_1.8.9          
##  [46] BiocNeighbors_2.1.0      iterators_1.0.14         foreach_1.5.2           
##  [49] tools_4.5.0              Rcpp_1.0.13              glue_1.8.0              
##  [52] gridExtra_2.3            SparseArray_1.7.0        xfun_0.48               
##  [55] dplyr_1.1.4              HDF5Array_1.35.0         gypsum_1.3.0            
##  [58] withr_3.0.2              BiocManager_1.30.25      fastmap_1.2.0           
##  [61] rhdf5filters_1.19.0      fansi_1.0.6              digest_0.6.37           
##  [64] rsvd_1.0.5               R6_2.5.1                 colorspace_2.1-1        
##  [67] RSQLite_2.3.7            utf8_1.2.4               generics_0.1.3          
##  [70] rtracklayer_1.67.0       FNN_1.1.4.1              httr_1.4.7              
##  [73] S4Arrays_1.7.0           uwot_0.2.2               pkgconfig_2.0.3         
##  [76] gtable_0.3.6             blob_1.2.4               XVector_0.47.0          
##  [79] htmltools_0.5.8.1        bookdown_0.41            ProtGenerics_1.39.0     
##  [82] scales_1.3.0             alabaster.matrix_1.7.0   ClusterR_1.3.3          
##  [85] png_0.1-8                knitr_1.48               rjson_0.2.23            
##  [88] curl_5.2.3               cachem_1.1.0             rhdf5_2.51.0            
##  [91] BiocVersion_3.21.1       parallel_4.5.0           vipor_0.4.7             
##  [94] AnnotationDbi_1.69.0     restfulr_0.0.15          pillar_1.9.0            
##  [97] grid_4.5.0               alabaster.schemas_1.7.0  vctrs_0.6.5             
## [100] BiocSingular_1.23.0      dbplyr_2.5.0             beachmat_2.23.0         
## [103] cluster_2.1.6            beeswarm_0.4.0           evaluate_1.0.1          
## [106] tinytex_0.53             GenomicFeatures_1.59.0   magick_2.8.5            
## [109] cli_3.6.3                locfit_1.5-9.10          compiler_4.5.0          
## [112] Rsamtools_2.23.0         rlang_1.1.4              crayon_1.5.3            
## [115] mbkmeans_1.23.0          labeling_0.4.3           ggbeeswarm_0.7.2        
## [118] alabaster.se_1.7.0       viridisLite_0.4.2        BiocParallel_1.41.0     
## [121] munsell_0.5.1            Biostrings_2.75.0        lazyeval_0.2.2          
## [124] Matrix_1.7-1             ExperimentHub_2.15.0     benchmarkme_1.0.8       
## [127] bit64_4.5.2              Rhdf5lib_1.29.0          KEGGREST_1.47.0         
## [130] statmod_1.5.0            alabaster.ranges_1.7.0   highr_0.11              
## [133] apcluster_1.4.13         AnnotationHub_3.15.0     igraph_2.1.1            
## [136] memoise_2.0.1            bslib_0.8.0              benchmarkmeData_1.0.4   
## [139] bit_4.5.0