Contents

1 Introduction

The bluster package implements a variety of diagnostics to quantify the behavior of a clustering method. These diagnostics are most helpful when they can be computed for each cluster, allowing us to prioritize certain clusters for more careful analysis. Poorly separated clusters can be prioritized for manual inspection to determine the relevance of the differences to its neighboring clusters, while heterogeneous clusters may be subject to further clustering to identify internal structure. We can also quantify the effect of different algorithms or parameter values, providing some insight into how our dataset responds to changes in the clustering method.

Let’s 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()

# Performing all of the preprocessing to get to the PCA,
# but skipping t-SNE, UMAP and the internal clustering.
library(scrapper)
sce <- analyze.se(sce,
    more.rna.qc.args=c(altexp.proportions="ERCC"),
    more.tsne.args=NULL,
    more.umap.args=NULL,
    more.cluster.graph.args=NULL
)$x

# Clustering with bluster.
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  13  14  15  16  17 
## 193 101 165 125  46 128 125  44  45  63  71  22  15  25  42  57  24

2 Computing the silhouette width

The silhouette width is an established metric for evaluating cluster separation. For each cell, we compute the average distance to all cells in the same cluster. We also compute the average distance to all cells in another cluster, taking the minimum of the averages across all other clusters. The silhouette width for each cell is defined as the difference between these two values divided by their maximum. Cells with large positive silhouette widths are closer to other cells in the same cluster than to cells in different clusters. Thus, clusters with large positive silhouette widths are well-separated from other clusters.

The silhouette width is a natural diagnostic for hierarchical clustering where the distance matrix is already available. For larger datasets, we instead use an approximate approach that uses the root of the average squared distances rather than the average distance itself. The approximation avoids the time-consuming calculation of pairwise distances that would otherwise make this metric impractical. 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            3        7   0.157041
## D2ex_2            3        7   0.161503
## D2ex_3            3        7   0.236920
## D2ex_4            4        1  -0.205589
## D2ex_5            4        1  -0.113506
## ...             ...      ...        ...
## D17TGFB_91        5       2   0.0944543
## D17TGFB_92        7       3   0.2468357
## D17TGFB_93        7       3  -0.0749748
## D17TGFB_94        5       11  0.2150985
## D17TGFB_95        5       2   0.1722279
boxplot(split(sil$width, clusters), xlab="cluster", ylab="silhouette")

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  13  14  15  16  17
##       1  178   0   0   0   0   0   0   0   0   0   0   0   0   0  15   0   0
##       2    0  97   0   0   0   0   0   0   0   0   4   0   0   0   0   0   0
##       3    0   0 143   0   0   0  22   0   0   0   0   0   0   0   0   0   0
##       4   66   0   3  17   0   0   3   4   0   0   0   0   7   3  21   1   0
##       5    0  20   0   0  20   0   0   0   0   0   6   0   0   0   0   0   0
##       6    0   1   0   0   2 116   0   1   0   0   0   0   0   0   0   8   0
##       7    0   0  10   0   0   0 114   0   0   0   0   0   0   0   1   0   0
##       8    0   0   0   0   0   0   0  33   0   0   0   0   9   2   0   0   0
##       9    0   0   6   0   0   0   2   0  14   0   0   0   1  13   0   2   7
##       10   0   0   0   0   0   0   0   0   0  63   0   0   0   0   0   0   0
##       11   0   0   0   0   0   0   0   0   0   0  71   0   0   0   0   0   0
##       12   0   0   0   0   0   0   0   0   0   0   0  22   0   0   0   0   0
##       13   0   0   0   0   0   0   0   0   0   0   0   0  15   0   0   0   0
##       14   0   0   0   0   0   0   0   0   0   0   0   0   0  25   0   0   0
##       15   0   0   0   0   0   0   0   0   0   0   0   0   0   0  42   0   0
##       16   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  57   0
##       17   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  24

A useful aspect of the silhouette width is that it naturally captures both underclustering and overclustering. Cells in heterogeneous clusters will have a large average distance to other cells in the same cluster; all else being equal, this will decrease their widths compared to cells in homogeneous clusters. Conversely, cells in overclustered datasets will have low average distances to cells in adjacent clusters, again resulting in low widths. One can exploit this to obtain a “sensible” initial choice for the number of clusters by maximizing the average silhouette width, (see comments in Section 8).

3 Computing the neighborhood purity

Another metric to assess cluster separation is the degree to which cells from multiple clusters intermingle in expression space. The “clustering purity” is defined for each cell as the proportion of neighboring cells that are assigned to the same cluster, after some weighting to adjust for differences in the number of cells between clusters. Well-separated clusters should exhibit little intermingling and thus high purity values for all member cells. 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      0.824468        3
## D2ex_2      0.872093        3
## D2ex_3      0.886728        3
## D2ex_4      0.467312        1
## D2ex_5      1.000000        4
## ...              ...      ...
## D17TGFB_91  0.611599        5
## D17TGFB_92  0.923481        7
## D17TGFB_93  0.599359        7
## D17TGFB_94  0.918395        5
## D17TGFB_95  0.624503        5
boxplot(split(pure$purity, clusters), xlab="cluster", ylab="purity")

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  13  14  15  16  17
##       1  187   0   0   1   0   0   0   0   0   0   0   0   0   0   5   0   0
##       2    0  98   0   0   0   0   0   0   0   0   3   0   0   0   0   0   0
##       3    0   0 160   0   0   0   5   0   0   0   0   0   0   0   0   0   0
##       4   10   0   0 111   0   0   0   0   0   0   0   0   1   0   3   0   0
##       5    0  17   0   0  25   0   0   0   0   0   4   0   0   0   0   0   0
##       6    0   0   0   0   0 120   0   0   0   1   0   0   0   1   0   6   0
##       7    0   0   4   0   0   0 121   0   0   0   0   0   0   0   0   0   0
##       8    0   0   0   0   0   0   0  44   0   0   0   0   0   0   0   0   0
##       9    0   0   0   0   0   0   0   0  44   0   0   0   0   0   0   0   1
##       10   0   0   0   0   0   0   0   0   0  63   0   0   0   0   0   0   0
##       11   0  25   0   0   0   0   0   0   0   0  46   0   0   0   0   0   0
##       12   0   0   0   0   0   0   0   0   0   0   0  22   0   0   0   0   0
##       13   0   0   0   0   0   0   0   0   0   0   0   0  15   0   0   0   0
##       14   0   0   0   0   0   0   0   0   0   0   0   1   0  24   0   0   0
##       15   0   0   0   0   0   0   0   0   0   0   0   0   0   0  42   0   0
##       16   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  57   0
##       17   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  24

The main difference between the cluster purity and silhouette width is that the former ignores the intra-cluster variance. This provides a simpler interpretation of cluster separation; a low silhouette width may still occur in well-separated clusters if the internal heterogeneity is high, while no such complication exists for the cluster purity. Comparing these two metrics can give us an indication of which clusters are heterogeneous, well-separated, or both.

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. Of course, this is subject to a number of caveats. Clusters generated from cells with low library sizes will naturally have larger RMSDs due to the greater impact of sequencing noise. Immediately adjacent clusters may also have high RMSDs if there are many cells on the boundaries between clusters, inflating the sum of squares from the center.

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

Alternatively, we can compute the within-cluster sum of squares (WCSS). This is a natural diagnostic for \(k\)-means clustering given that the algorithm aims to find a clustering that minimizes the WCSS. 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 
## 51886.793  7803.356 45233.285 64782.391 10212.559 23472.442 33808.943 22838.157 
##         9        10        11        12        13        14        15        16 
## 29106.698  9554.349  6210.724  2943.964  3597.425  6238.233  8256.786  8154.253 
##        17 
##  3072.608

5 Computing graph modularity

For graph-based methods, we can compute the cluster modularity within clusters and between pairs of clusters. The pairwiseModularity() function computes the ratio between the observed sum of edge weights and the expected sum for random edges, either for the same cluster (on the diagonal) or for two different clusters (off-diagonal). We use the ratio instead of the difference as the former is less sensitive to the number of cells in each cluster.

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

Well-defined clusters should have most of its edges occurring between cells in the same cluster, manifesting as a large modularity ratio on the diagonal. Conversely, large off-diagonal entries indicate that the corresponding pair of clusters are closely related as they are connected by many edges.

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

One useful approach is to use the ratio matrix to form another graph where the nodes are clusters rather than cells. Edges between nodes are weighted according to the ratio of observed to expected edge weights between cells in those clusters. We can then repeat our graph operations on this new cluster-level graph to explore the relationships between clusters. For example, we could obtain clusters of clusters, or we could simply create a new cluster-based layout for visualization:

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
##   3   4   6   7   8  11  12  15  16  17 
## 165 125 128 125  59 218  85 235  82  69

6 Comparing two clusterings

We can quantify the agreement between two clusterings by computing the Rand index with pairwiseRand(). This is defined as the proportion of pairs of cells that retain the same status (i.e., both cells in the same cluster, or each cell in different clusters) in both clusterings. In practice, we usually compute the adjusted Rand index (ARI) where we subtract the number of concordant pairs expected under random permutations of the clusterings; this accounts for differences in the size and number of clusters within and between clusterings. A larger ARI indicates that the clusters are preserved, up to a maximum value of 1 for identical clusterings. In and of itself, the magnitude of the ARI has little meaning, and it is best used to assess the relative similarities of different clusterings (e.g., “Walktrap is more similar to Louvain than either are to Infomap”). Nonetheless, if one must have a hard-and-fast rule, experience suggests that an ARI greater than 0.5 corresponds to “good” similarity between two clusterings.

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

pairwiseRand() can also 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.

adj.ratio <- pairwiseRand(clusters, hclusters, mode="ratio")
library(pheatmap)
pheatmap(adj.ratio, cluster_cols=FALSE, cluster_rows=FALSE,
    col=viridis::viridis(100), breaks=seq(-1, 1, length=101))

# Alternatively, we can skip the adjustment to get values in [0, 1]:
ratio <- pairwiseRand(clusters, hclusters, mode="ratio", adjusted=FALSE)
library(pheatmap)
pheatmap(ratio, cluster_cols=FALSE, cluster_rows=FALSE,
    col=viridis::viridis(100), breaks=seq(0, 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

A desirable property of a given clustering is that it is stable to perturbations to the input data. Stable clusters are logistically convenient as small changes to upstream processing will not change the conclusions. Greater stability also increases the likelihood that those conclusions can be reproduced in an independent replicate study. We use bootstrapping to evaluate the stability of a clustering algorithm on a given dataset. That is, cells are sampled with replacement to create a “bootstrap replicate” dataset, and clustering is repeated on this replicate to see if the same clusters can be reproduced.

set.seed(1001010)
ari <-bootstrapStability(
    mat,
    clusters=clusters, 
    mode="ratio",
    BLUSPARAM=NNGraphParam()
)
ari
##            1         2         3          4           5         6           7
## 1  0.7000951 1.0000000 1.0000000 -0.5194724  1.00000000 1.0000000 1.000000000
## 2         NA 0.9663014 1.0000000  1.0000000 -0.06161268 1.0000000 1.000000000
## 3         NA        NA 0.8418279  0.8903865  1.00000000 1.0000000 0.006845124
## 4         NA        NA        NA  0.4616849  1.00000000 1.0000000 0.649037481
## 5         NA        NA        NA         NA  0.58358961 0.9088248 1.000000000
## 6         NA        NA        NA         NA          NA 0.8944236 1.000000000
## 7         NA        NA        NA         NA          NA        NA 0.906927260
## 8         NA        NA        NA         NA          NA        NA          NA
## 9         NA        NA        NA         NA          NA        NA          NA
## 10        NA        NA        NA         NA          NA        NA          NA
## 11        NA        NA        NA         NA          NA        NA          NA
## 12        NA        NA        NA         NA          NA        NA          NA
## 13        NA        NA        NA         NA          NA        NA          NA
## 14        NA        NA        NA         NA          NA        NA          NA
## 15        NA        NA        NA         NA          NA        NA          NA
## 16        NA        NA        NA         NA          NA        NA          NA
## 17        NA        NA        NA         NA          NA        NA          NA
##            8         9 10         11 12        13        14         15
## 1  1.0000000 0.9861999  1  1.0000000  1 1.0000000 1.0000000 -1.9271731
## 2  1.0000000 1.0000000  1  0.3346691  1 1.0000000 1.0000000  1.0000000
## 3  1.0000000 0.7224516  1  1.0000000  1 1.0000000 1.0000000  1.0000000
## 4  0.8734594 0.1149303  1  1.0000000  1 1.0000000 1.0000000  0.4718479
## 5  1.0000000 1.0000000  1 -2.7422963  1 1.0000000 1.0000000  1.0000000
## 6  1.0000000 1.0000000  1  1.0000000  1 1.0000000 0.9784837  1.0000000
## 7  1.0000000 0.9784937  1  1.0000000  1 1.0000000 1.0000000  1.0000000
## 8  0.8425927 0.9250305  1  1.0000000  1 0.1722527 1.0000000  1.0000000
## 9         NA 0.5364578  1  1.0000000  1 1.0000000 1.0000000  1.0000000
## 10        NA        NA  1  1.0000000  1 1.0000000 1.0000000  1.0000000
## 11        NA        NA NA  1.0000000  1 1.0000000 1.0000000  1.0000000
## 12        NA        NA NA         NA  1 1.0000000 1.0000000  1.0000000
## 13        NA        NA NA         NA NA 1.0000000 1.0000000  1.0000000
## 14        NA        NA NA         NA NA        NA 1.0000000  1.0000000
## 15        NA        NA NA         NA NA        NA        NA  1.0000000
## 16        NA        NA NA         NA NA        NA        NA         NA
## 17        NA        NA NA         NA NA        NA        NA         NA
##            16        17
## 1   1.0000000 1.0000000
## 2   1.0000000 1.0000000
## 3   1.0000000 1.0000000
## 4   1.0000000 1.0000000
## 5   1.0000000 1.0000000
## 6   0.2325991 1.0000000
## 7   1.0000000 1.0000000
## 8   1.0000000 1.0000000
## 9   1.0000000 0.9114592
## 10  1.0000000 1.0000000
## 11  1.0000000 1.0000000
## 12  1.0000000 1.0000000
## 13  1.0000000 1.0000000
## 14 -0.7236293 1.0000000
## 15  1.0000000 1.0000000
## 16  1.0000000 1.0000000
## 17         NA 1.0000000

The function returns a matrix of ARI-derived ratios for every pair of original clusters in clusters, averaged across bootstrap iterations. High ratios indicate that the clustering in the bootstrap replicates are highly consistent with that of the original dataset. More specifically, high ratios on the diagonal indicate that cells in the same original cluster are still together in the bootstrap replicates, while high ratios off the diagonal indicate that cells in the corresponding cluster pair are still separated.

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

Bootstrapping is a general approach for evaluating cluster stability that is compatible with any clustering algorithm. The ARI-derived ratio between cluster pairs is also more informative than a single stability measure for all/each cluster as the former considers the relationships between clusters, e.g., unstable separation between \(X\) and \(Y\) does not penalize the stability of separation between \(X\) and another cluster \(Z\). However, one should take these metrics with a grain of salt as bootstrapping only considers the effect of sampling noise, ignoring other factors that affect reproducibility in an independent study (e.g., batch effects, donor variation). In addition, it is possible for a poor separation to be highly stable, so a highly stable cluster may not necessarily represent some distinct subpopulation.

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. We could then use linkClusters(), clustree or similar functions to visualize the relationships between different clusterings.

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. We might then make some decisions on which clustering(s) to use for further analyses. For example, we might choose a few clusterings that span several different resolutions, so as to obtain a greater diversity of summaries of the data; conversely, we might be able to save some time and effort by ignoring redundant clusterings with similar values for our metrics.

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

We could even use the sweep to automatically choose the “best” clustering by optimizing one or more of these metrics. The simplest strategy is to maximize the silhouette width, though one can imagine more complex scores involving combinations of metrics. This approach is valid but any automatic choice should be treated as a suggestion rather than a rule. The clustering at the optimal value of a metric may not be the most scientifically informative clustering, given that well-separated clusters typically correspond to cell types that are already known. Conversely, poorly separated clusters will often be observed in non-trivial analyses of scRNA-seq data where the aim is to characterize closely related subtypes or states.

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 71127a2 UNW- 48 92 -- 
## + attr: name (v/c), weight (e/n)
## + edges from 71127a2 (vertex names):
##  [1] walktrap.3 --louvain.1 walktrap.7 --louvain.1 walktrap.4 --louvain.2
##  [4] walktrap.1 --louvain.3 walktrap.4 --louvain.3 walktrap.15--louvain.3
##  [7] walktrap.3 --louvain.4 walktrap.4 --louvain.4 walktrap.8 --louvain.4
## [10] walktrap.9 --louvain.4 walktrap.13--louvain.4 walktrap.17--louvain.4
## [13] walktrap.6 --louvain.5 walktrap.14--louvain.5 walktrap.16--louvain.5
## [16] walktrap.10--louvain.6 walktrap.12--louvain.6 walktrap.2 --louvain.7
## [19] walktrap.5 --louvain.7 walktrap.6 --louvain.7 walktrap.11--louvain.7
## [22] walktrap.4 --louvain.8 walktrap.7 --louvain.8 walktrap.3 --infomap.1
## + ... omitted several edges

The output is a graph where edges are formed between related clusters in the different clusterings. Each edge is weighted according to the Jaccard indices of the two clusters, i.e., the ratio of the number of shared cells over the total number of unique cells. Larger indices represent a stronger correspondence between the two clusters. 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)

We can identify the best corresponding clusters based on the Jaccard index. The magnitude of the index can be used as a measure of strength for the correspondence between those two clusters. A low index for a cluster indicates that no counterpart exists in the other clustering.

link.mat <- linkClustersMatrix(
    combinations$clusters$k.10_cluster.fun.walktrap,
    combinations$clusters$k.10_cluster.fun.louvain
)
best <- max.col(link.mat, ties.method="first")
DataFrame(
    Cluster=rownames(link.mat), 
    Corresponding=colnames(link.mat)[best], 
    Index=link.mat[cbind(seq_len(nrow(link.mat)), best)]
)
## DataFrame with 17 rows and 3 columns
##         Cluster Corresponding     Index
##     <character>   <character> <numeric>
## 1             1             3  0.750973
## 2             2             7  0.459091
## 3             3             1  0.964497
## 4             4             2  0.672000
## 5             5             7  0.209091
## ...         ...           ...       ...
## 13           13             4  0.102740
## 14           14             5  0.120192
## 15           15             3  0.163424
## 16           16             5  0.274038
## 17           17             4  0.164384

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_from_adjacency_matrix(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 version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 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.13             scater_1.41.2              
##  [3] ggplot2_4.0.3               scuttle_1.23.1             
##  [5] bluster_1.23.1              scrapper_1.7.8             
##  [7] scRNAseq_2.27.0             SingleCellExperiment_1.35.2
##  [9] SummarizedExperiment_1.43.0 Biobase_2.73.2             
## [11] GenomicRanges_1.65.1        Seqinfo_1.3.0              
## [13] IRanges_2.47.2              S4Vectors_0.51.7           
## [15] BiocGenerics_0.59.12        generics_0.1.4             
## [17] MatrixGenerics_1.25.0       matrixStats_1.5.0          
## [19] BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3       jsonlite_2.0.0           magrittr_2.0.5          
##   [4] magick_2.9.1             ggbeeswarm_0.7.3         GenomicFeatures_1.65.0  
##   [7] gypsum_1.9.0             farver_2.1.2             rmarkdown_2.31          
##  [10] BiocIO_1.23.3            vctrs_0.7.3              memoise_2.0.1           
##  [13] Rsamtools_2.29.0         RCurl_1.98-1.20          benchmarkme_1.0.8       
##  [16] tinytex_0.60             htmltools_0.5.9          S4Arrays_1.13.0         
##  [19] BiocBaseUtils_1.15.1     dynamicTreeCut_1.63-1    AnnotationHub_4.3.2     
##  [22] curl_8.0.0               BiocNeighbors_2.7.3      Rhdf5lib_2.1.0          
##  [25] SparseArray_1.13.2       rhdf5_2.57.12            sass_0.4.10             
##  [28] alabaster.base_1.13.2    bslib_0.12.0             alabaster.sce_1.13.0    
##  [31] httr2_1.3.0              cachem_1.1.0             GenomicAlignments_1.49.1
##  [34] igraph_2.3.3             iterators_1.0.14         lifecycle_1.0.5         
##  [37] pkgconfig_2.0.3          rsvd_1.0.5               Matrix_1.7-6            
##  [40] R6_2.6.1                 fastmap_1.2.0            digest_0.6.39           
##  [43] AnnotationDbi_1.75.2     irlba_2.3.7              ExperimentHub_3.3.2     
##  [46] RSQLite_3.53.3           beachmat_2.29.1          labeling_0.4.3          
##  [49] filelock_1.0.3           httr_1.4.8               abind_1.4-8             
##  [52] compiler_4.6.1           doParallel_1.0.17        bit64_4.8.4             
##  [55] withr_3.0.3              S7_0.2.2                 BiocParallel_1.47.0     
##  [58] viridis_0.6.5            DBI_1.3.0                apcluster_1.4.14        
##  [61] HDF5Array_1.41.3         alabaster.ranges_1.13.0  alabaster.schemas_1.13.0
##  [64] rappdirs_0.3.4           DelayedArray_0.39.6      rjson_0.2.23            
##  [67] tools_4.6.1              vipor_0.4.7              otel_0.2.0              
##  [70] mbkmeans_1.29.0          beeswarm_0.4.0           glue_1.8.1              
##  [73] h5mread_1.5.2            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.0      ScaledMatrix_1.21.0     
##  [82] XVector_0.53.0           foreach_1.5.2            ggrepel_0.9.8           
##  [85] BiocVersion_3.24.0       pillar_1.11.1            benchmarkmeData_2.0.0   
##  [88] dplyr_1.2.1              BiocFileCache_3.3.0      lattice_0.23-1          
##  [91] gmp_0.7-5.1              rtracklayer_1.73.0       bit_4.6.0               
##  [94] tidyselect_1.2.1         Biostrings_2.81.6        knitr_1.51              
##  [97] gridExtra_2.3.1          bookdown_0.47            ProtGenerics_1.45.0     
## [100] xfun_0.60                UCSC.utils_1.9.0         lazyeval_0.2.3          
## [103] yaml_2.3.12              evaluate_1.0.5           codetools_0.2-20        
## [106] cigarillo_1.3.1          tibble_3.3.1             alabaster.matrix_1.13.1 
## [109] BiocManager_1.30.27      cli_3.6.6                jquerylib_0.1.4         
## [112] dichromat_2.0-1          Rcpp_1.1.2               GenomeInfoDb_1.49.1     
## [115] dbplyr_2.6.0             png_0.1-9                kohonen_3.0.13          
## [118] XML_3.99-0.24            parallel_4.6.1           blob_1.3.0              
## [121] ClusterR_1.3.6           AnnotationFilter_1.37.0  bitops_1.1-0            
## [124] viridisLite_0.4.3        alabaster.se_1.13.0      scales_1.4.0            
## [127] crayon_1.5.3             rlang_1.3.0              cowplot_1.2.0           
## [130] KEGGREST_1.53.6