---
title: Assorted clustering diagnostics
author: 
- name: Aaron Lun
  email: infinite.monkeys.with.keyboards@gmail.com
output: BiocStyle::html_document
package: bluster
vignette: >
  %\VignetteIndexEntry{2. Clustering diagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}    
---

```{r, echo=FALSE}
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
library(BiocStyle)
```

# Introduction 

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

```{r}
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)
```

# 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 `r Biocpkg("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.

```{r}
sil <- approxSilhouette(mat, clusters)
sil
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.

```{r}
best.choice <- ifelse(sil$width > 0, clusters, sil$other)
table(Assigned=clusters, Closest=best.choice)
```

# 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.

```{r}
pure <- neighborPurity(mat, clusters)
pure
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.

```{r}
table(Assigned=clusters, Max=pure$maximum)
```

# 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.

```{r}
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 \@ref(clustering-parameter-sweeps) for more details.

```{r}
clusterRMSD(mat, clusters, sum=TRUE)
```

# 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 `r CRANpkg("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.

```{r}
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: 

```{r}
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:

```{r}
merged <- mergeCommunities(g, clusters, number=10)
table(merged)
```

# 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.

```{r}
hclusters <- clusterRows(mat, HclustParam(cut.dynamic=TRUE))
pairwiseRand(clusters, hclusters, mode="index")
```

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.

```{r}
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 \@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

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.

```{r}
set.seed(1001010)
ari <-bootstrapStability(mat, clusters=clusters, 
    mode="index", BLUSPARAM=NNGraphParam())
ari
```

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

```{r}
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))
```

# 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.

```{r}
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.

```{r}
colnames(combinations$clusters)
combinations$parameters
```

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.

```{r}
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")
```

# 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:

```{r}
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
```

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

```{r}
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.

# 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.

```{r}
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 {-}

```{r}
sessionInfo()
```
