Assorted clustering diagnostics
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 24Computing 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.
## 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
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 24A 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 @ref(clustering-parameter-sweeps)).
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.
## 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
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.
## 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 24The 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.
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.
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 @ref(clustering-parameter-sweeps) for more details.
## 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.608Computing 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.
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
## 3 4 6 7 8 11 12 15 16 17
## 165 125 128 125 59 218 85 235 82 69Comparing 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.4385518pairwiseRand() 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 @ref(linking-clusters)).
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.
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.0000000The 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.
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.
## [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"## 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 infomapWe 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.
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 cae79ab UNW- 50 98 --
## + attr: name (v/c), weight (e/n)
## + edges from cae79ab (vertex names):
## [1] walktrap.3 --louvain.1 walktrap.4 --louvain.1 walktrap.7 --louvain.1
## [4] walktrap.9 --louvain.1 walktrap.4 --louvain.2 walktrap.1 --louvain.3
## [7] walktrap.4 --louvain.3 walktrap.15--louvain.3 walktrap.4 --louvain.4
## [10] walktrap.8 --louvain.4 walktrap.13--louvain.4 walktrap.6 --louvain.5
## [13] walktrap.10--louvain.6 walktrap.12--louvain.6 walktrap.14--louvain.7
## [16] walktrap.16--louvain.7 walktrap.2 --louvain.8 walktrap.5 --louvain.8
## [19] walktrap.6 --louvain.8 walktrap.11--louvain.8 walktrap.9 --louvain.9
## [22] walktrap.17--louvain.9 walktrap.3 --louvain.10 walktrap.4 --louvain.10
## + ... omitted several edgesThe 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:
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.756863
## 2 2 8 0.459091
## 3 3 1 0.792271
## 4 4 2 0.680000
## 5 5 8 0.209091
## ... ... ... ...
## 13 13 4 0.208333
## 14 14 7 0.304878
## 15 15 3 0.164706
## 16 16 7 0.695122
## 17 17 9 0.705882Comparing 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
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 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] 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 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.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] htmltools_0.5.9 S4Arrays_1.13.0 BiocBaseUtils_1.15.1
## [19] dynamicTreeCut_1.63-1 AnnotationHub_4.3.2 curl_7.1.0
## [22] BiocNeighbors_2.7.3 Rhdf5lib_2.1.0 SparseArray_1.13.2
## [25] rhdf5_2.57.12 sass_0.4.10 alabaster.base_1.13.2
## [28] bslib_0.12.0 alabaster.sce_1.13.0 httr2_1.3.0
## [31] cachem_1.1.0 buildtools_1.0.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 maketools_1.3.2 Biostrings_2.81.6
## [97] knitr_1.51 gridExtra_2.3.1 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] Rcpp_1.1.2 GenomeInfoDb_1.49.1 dbplyr_2.6.0
## [115] png_0.1-9 kohonen_3.0.13 XML_3.99-0.24
## [118] parallel_4.6.1 blob_1.3.0 ClusterR_1.3.6
## [121] AnnotationFilter_1.37.0 bitops_1.1-0 viridisLite_0.4.3
## [124] alabaster.se_1.13.0 scales_1.4.0 crayon_1.5.3
## [127] rlang_1.3.0 KEGGREST_1.53.6