1 Overview

Some of the materials originate in the Hemberg group course material with some of the text copied with a few edits. Also see the OSCA book’s “Basic” and “Advanced” chapters on clustering. In particular, please read the overview with regard to the comments on the “correctness” of any given clustering result.

Once we have normalized the data and removed confounders we can carry out analyses that are relevant to the biological questions at hand. The exact nature of the analysis depends on the data set. One of the most promising applications of scRNA-seq is de novo discovery and annotation of cell-types based on transcription profiles. This requires the identification of groups of cells based on the similarities of the transcriptomes without any prior knowledge of the label a.k.a. unsupervised clustering. To avoid the challenges caused by the noise and high dimensionality of the scRNA-seq data, clustering is performed after feature selection and dimensionality reduction. For data that has not required batch correction this would usually be based on the PCA output. As our data has required batch correction we will use the “corrected” reducedDims data.

We will focus here on graph-based clustering, however, it is also possible to apply hierarchical clustering and k-means clustering on smaller data sets - see the OSCA book for details. Graph-base clustering is a more recent development and better suited for scRNA-seq, especially large data sets.

2 Load the data

library(scater)
library(scran)
library(bluster)
library(cluster)
library(igraph)
library(pheatmap)
library(patchwork)
library(tidyverse)

We will use the data set generated in the previous session. This contains 7 samples from the Caron data set. For the purposes of these materials, in the interests of time, each sample has been downsampled to only contain 500 cells.

sce <- readRDS("R_objects/Caron_batch_corrected.500.rds")

table(sce$SampleName)

## 
## ETV6-RUNX1_1 ETV6-RUNX1_2 ETV6-RUNX1_3 ETV6-RUNX1_4      PBMMC_1      PBMMC_2 
##          500          500          500          500          500          500 
##      PBMMC_3 
##          500

3 Graph-based clustering overview

Graph-based clustering entails building a nearest-neighbour (NN) graph using cells as nodes and their similarity as edges, then identifying ‘communities’ of cells within the network. A graph-based clustering method has three key parameters:

How many neighbors are considered when constructing the graph
What scheme is used to weight the edges
Which community detection algorithm is used to define the clusters

3.1 Connecting nodes (cells) based on nearest neighbours

Two types of NN graph may be used: “K nearest-neighbour” (KNN) and “shared nearest-neighbour” (SNN). In a KNN graph, two nodes (cells), say A and B, are connected by an edge if the distance between them is amongst the k smallest distances from A to other cells. In an SNN graph A and B are connected if the distance is amongst the k samllest distances from A to other cells and also among the k smallest distance from B to other cells.

In the figure above, if k is 5, then A and B would be connected in a KNN graph as B is one of the 5 closest cells to A, however, they would not be connected in an SNN graph as B has 5 other cells that are closer to it than A.

The value of k can be roughly interpreted as the anticipated size of the smallest subpopulation” (see scran’s buildSNNGraph() manual).

The plot below shows the same data set as a network built using three different numbers of neighbours: 5, 15 and 25 (from here).

3.2 Weighting the edges

The edges between nodes (cells) can be weighted based on the similarity of the cells; edges connecting cells that are more closely related will have a higher weight. The three common methods for this weighting are (see the bluster package documentation for the makeSNNGraph function):

rank - the weight is based on the highest rank of the shared nearest neighbours
number - the weight is based the number of nearest neighbours in common between the two cells
jaccard - the Jaccard index of the two cells’ sets of nearest neighbours.

3.3 Grouping nodes (cells) into clusters

Clusters are identified using an algorithm that interprets the connections of the graph to find groups of highly interconnected cells. A variety of different algorithms are available to do this, in these materials we will focus on three methods: walktrap, louvain and leiden. See the OSCA book for details of others available in scran.

3.4 Modularity

Several methods to detect clusters (‘communities’) in networks rely on a metric called “modularity”. For a given partition of cells into clusters, modularity measures how separated clusters are from each other, based on the difference between the observed and expected weight of edges between nodes. For the whole graph, the closer to 1 the better.

3.5 Pros and Cons of graph based clustering

Pros:
- fast and memory efficient (avoids the need to construct a distance matrix for all pairs of cells)
- no assumptions on the shape of the clusters or the distribution of cells within each cluster
- no need to specify a number of clusters to identify (but the size of the neighbourhood used affects the size of clusters)
Cons:
- loss of information beyond neighboring cells, which can affect community detection in regions with many cells.

4 Implementation

The implementation of clustering in R is carried out using functions from a number of different packages, in particular the bluster and igraph packages. scran provides a handy “wrapper” function clusterCells that allows us use a variety of different algorithms with one simple command.

By default clusterCells just returns a vector containing the cluster number for each cell. We can also retrieve the intermediate statistics (varying according to the algorithm used) and the SNN graph by specifying the bluster argument full = TRUE. If you are only interested in retrieving the clusters, this isn’t necessary but in this first instance we will retrieve the graph and visualise it. The default algorithm for clusterCells is Walktrap with k is set to 10 by default. The default edge weighting is “rank”.

clustering1 <- clusterCells(sce, use.dimred="corrected", full=TRUE)

This has defined 24 clusters with varying numbers of cells:

table(clustering1$clusters)

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 133 154 243  66 194  74  14 151 313  64  34 143 204 498 844  31  80  24  56  56 
##  21  22  23  24 
##  67  17  24  16

The number of cells in the data set is large and plotting all the cells would take too long, so we randomly choose 1000 nodes (cells) in the network before plotting the resulting smaller network. Adding sample data to the graph and plotting the results are done using the igraph package. Cells can be color-coded by sample type:

# extract the graph
snn.gr <- clustering1$objects$graph

# Add Sample group to vertices (nodes, ie cells)
V(snn.gr)$SampleGroup <- as.character(colData(sce)$SampleGroup)

# pick 1000 nodes randomly
set.seed(1423)
selectedNodes <- sample(3500, 1000)

# subset graph for these 1000 randomly chosen nodes
snn.gr.subset <- subgraph(snn.gr, selectedNodes)

# set colors for clusters
grps <-  V(snn.gr.subset)$SampleGroup
cols <- c("dodgerblue", "lightyellow")[as.numeric(factor(grps))]
names(cols) <- grps

# plot graph
plot.igraph(snn.gr.subset,
  layout = layout_with_fr(snn.gr.subset),
  vertex.size = 3, 
  vertex.label = NA,
  vertex.color = cols,
  frame.color = cols,
  main = "default parameters"
)

# add legend
legend('bottomright',
       legend=unique(names(cols)),
       pch=21,
       pt.bg=unique(cols),
       pt.cex=1, cex=.6, bty="n", ncol=1)

More commonly we will visualise the clusters by superimposing them on a t-SNE or UMAP plot. We can store the clusters in the sce object colData.

sce$Clusters1 <- clustering1$clusters
plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by="Clusters1",
               text_by = "Clusters1")

5 The Walktrap method

The walktrap method relies on short random walks (a few steps) through the network. These walks tend to be ‘trapped’ in highly-connected regions of the network. Node similarity is measured based on these walks. Nodes are first each assigned their own community. Pairwise distances are computed and the two closest communities are grouped. These steps are repeated a given number of times to produce a dendrogram. Hierarchical clustering is then applied to the distance matrix. The best partition is that with the highest modularity. The original article describing the algorithm is Pons P, Latapy M (2006) Computing communities in large networks using random walks. J Graph Algorithms Appl 10(2):191–218

Walktrap is the default algorithm for clusterCells, k is set to 10 by default and the default edge weighting is “rank”. To explicitly request a specific algorithm and to set the k to a different number of nearest neighbours, we use a SNNGraphParam object from the bluster package (which is the package clusterCells is using under the hood).

Let’s set the k to 15 but keep the other parameters the same. This time we will just return the clusters:

sce$walktrap15 <- clusterCells(sce, 
                           use.dimred = "corrected", 
                           BLUSPARAM = SNNGraphParam(k = 15, 
                                                     cluster.fun = "walktrap"))

This time we have defined 16 clustering. As a general rule, increasing k will tend to decrease the number of clusters (not always, but generally).

table(sce$walktrap15)

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 438 260 124  84 322  51 154 125  38 960 142 501  26  57 202  16

We can visualise the assignment of cells from different samples to the clusters using a heatmap. This gives us an overview of how well each cluster is represented across the samples and the replicates. Several clusters (2, 8, 9 and 16) are present in the PBMMC samples, but absent from the ETV6_RUNX1 samples for instance.

w15_table <- log(table(sce$walktrap15, sce$SampleName)+1)
pheatmap(w15_table, cluster_rows = TRUE, cluster_cols = FALSE)

Most clusters comprise cells from several replicates of the same sample type, cluster 10 appears to be predominantly cells from the ETV6-RUNX samples.

We can visualise this on the TSNE:

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by="walktrap15", 
               text_by = "walktrap15")

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by="walktrap15", 
               text_by = "walktrap15",
               other_fields = list("SampleGroup")) +
  facet_wrap(vars(SampleGroup))

The different clustering algorithms may have additional parameters, specific to the algorithm, that can be adjusted. With the walktrap algorithm we could also tweak the number of “steps” in each walk. The default is 4, but we could, for example, change this to 10 by adding the parameter cluster.args = list(steps = 10) to the SNNGraphParam object in the clusterCells command.

6 The Louvain method

With the Louvain method, nodes are also first assigned their own community. This hierarchical agglomerative method then progresses in two-step iterations:

nodes are re-assigned one at a time to the community for which they increase modularity the most, if at all.
a new, ‘aggregate’ network is built where nodes are the communities formed in the previous step.

These two steps are repeated until modularity stops increasing. The diagram below is copied from this article.

We now apply the Louvain approach, store its outcome in the SCE object and show cluster sizes.

sce$louvain15 <- clusterCells(sce, 
                           use.dimred = "corrected", 
                           BLUSPARAM = SNNGraphParam(k = 15, 
                                                     cluster.fun = "louvain"))

table(sce$louvain15)

## 
##   1   2   3   4   5   6   7   8   9  10  11 
## 424  58 838 463 212 157 496 236 140 214 262

The t-SNE plot shows cells color-coded by cluster membership:

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by = "louvain15", 
               text_by = "louvain15")

If we split by sample type we can see differences in the clusters between the sample groups:

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by="louvain15", 
               text_by = "louvain15",
               other_fields = list("SampleGroup")) +
   facet_wrap(vars(SampleGroup))

7 The Leiden method

The Leiden method improves on the Louvain method by guaranteeing that at each iteration clusters are connected and well-separated. The method includes an extra step in the iterations: after nodes are moved (step 1), the resulting partition is refined (step2) and only then the new aggregate network made, and refined (step 3). The diagram below is copied from this article.

For this exercise please run the clustering again, this time using the “leiden” method.

Set the k to 20 and add the results of the clustering to the sce object in a new column called “leiden20”.

How many clusters does this result in?

Visualize the clusters by plotting the t-SNE with the cells coloured according to your new clustering.

Hint

You will need to change the k parameter and the cluster.fun parameter in the SNNGraphParam object used in the clusterCells function.

Answer

First run the clustering with clusterCells:

sce$leiden20 <- clusterCells(sce, 
                           use.dimred = "corrected", 
                           BLUSPARAM = SNNGraphParam(k = 20, 
                                                     cluster.fun = "leiden"))

We can quickly look at the results by summarising using table.

table(sce$leiden20)

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15 
## 941 332  62 428 209 508   3   1 147 120 116 140 214 262  17

There are 15 clusters, although cluster 7 contains only 3 cells and cluster 8 contains only 1 cell.

The t-SNE plot shows cells color-coded by cluster membership:

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by = "leiden20", 
               text_by = "leiden20")

8 Assessing cluster behaviour

A variety of metrics are available to aid us in assessing the behaviour of a particular clustering method on our data. These can help us in assessing how well defined different clusters within a single clustering are in terms of the relatedness of cells within the cluster and the how well separated that cluster is from cells in other clusters, and to compare the results of different clustering methods or parameter values (e.g. different values for k).

We will consider “Silhouette width” and “Modularity”. Further details and other metrics are described in the “Advanced” section of the OSCA book.

8.1 Silhouette width

The silhouette width (so named after the look of the traditional graph for plotting the results) is a measure of how closely related cells within cluster are to one another versus how closely related cells in the cluster are to cells in other clusters. This allows us to assess cluster separation.

For each cell in the cluster we calculate the the average distance to all other cells in the cluster and the average distance to all cells not in the cluster. The cells silhouette width is the difference between these divided by the maximum of the two values. Cells with a large silhouette are strongly related to cells in the cluster, cells with a negative silhouette width are more closely related to other clusters.

We will use the approxSilhouette function from the bluster package. The resulting table gives us the silhouette width for each cell, the cluster it belongs to, and which other cluster it is most closely related to.

sil.approx <- approxSilhouette(reducedDim(sce, "corrected"),
                               clusters=sce$leiden20)
sil.approx

## DataFrame with 3500 rows and 3 columns
##                        cluster    other      width
##                       <factor> <factor>  <numeric>
## 1_CGACTTCGTCCAGTTA-1         1        2 -0.0131760
## 1_AGAATAGCATACGCTA-1         2        1  0.0161969
## 1_TGACTAGAGAACTCGG-1         2        1  0.2809436
## 1_CTTAACTGTTATGCGT-1         3        8  0.3387946
## 1_CCCAGTTTCAAGCCTA-1         1        2  0.2535488
## ...                        ...      ...        ...
## 11_CCCAGTTCACATCCGG-1       14        7  0.0653069
## 11_TTAGTTCGTTAAAGAC-1       13        9  0.5494922
## 11_GATGAAATCTGGGCCA-1       6         9  0.4776268
## 11_ACGGAGAGTTAAGATG-1       13        2  0.2663340
## 11_TGGCGCAAGCGTTGCC-1       13        2  0.5178875

We can view the results in as a beeswarm plot. We colour each cell according to either its current cluster or, if the cell has a negative silhouette width, the cluster that it is closest to.

plotSilBeeswarm <- function(silDat){
  silTab <- silDat %>% 
    as.data.frame() %>% 
    mutate(closestCluster = ifelse(width > 0, cluster, other) %>% factor())
  
  plt <- silTab %>% 
      ggplot(aes(x=cluster, y=width, colour=closestCluster)) +
        ggbeeswarm::geom_quasirandom(method="smiley", alpha=0.6) +
        theme_bw()
  
  plt <- scater:::.resolve_plot_colours(plt, silTab$closestCluster, "closestCluster")
  plt
}

p1 <- plotSilBeeswarm(sil.approx)
p2 <- plotReducedDim(sce, 
                     dimred = "TSNE_corrected", 
                     colour_by="leiden20", 
                     text_by = "leiden20")
p1 + p2

We could also look at the correspondence between different clusters by plotting these numbers on a grid showing for each cluster number of cells in that cluster that are closer to another cluster, colouring each tile by the proportion of the total cells in the cluster that it contains. Ideally we would like to see a strong diagonal band and only a few off-diagonal tiles containing small number of cells.

plotSilGrid <- function(silDat){
  silDat %>% 
    as.data.frame() %>% 
    mutate(closestCluster = ifelse(width > 0, cluster, other) %>% factor()) %>% 
    count(cluster, closestCluster,  name="olap") %>% 
    group_by(cluster) %>% 
    mutate(total  = sum(olap)) %>% 
    mutate(proportion = olap / total) %>% 
    mutate(proportion = ifelse(cluster == closestCluster, proportion, -proportion)) %>% 
    ggplot(aes(x = cluster, y = closestCluster)) +
      geom_tile(aes(fill = proportion)) +
      geom_text(aes(label = olap), size=5) +
      scale_fill_gradientn(colors = c("#fc8d59", "#ffffbf", "#91cf60"),
                            limits = c(-1, 1)) +
      geom_vline(xintercept=seq(0.5, 30.5, by=1)) +
      geom_hline(yintercept=seq(0.5, 30.5, by=1), colour="lightgrey", linetype=2) +
      guides(fill = "none") +
      theme(
          aspect.ratio = 1,
          panel.background = element_blank())
}
plotSilGrid(sil.approx)

From these two plots we can see that clusters 7, 8 and 12 appear to have a good degree of separation, however, clusters 7 and 8 only contains few cells, whilst there are many cells in other clusters that appear closer to them than they are to their assigned cluster. Perhaps clusters 7 and 8 needs to be merged with cluster 1.

Let’s do the same plots with the walktrap clusters generated with k=15.

sil.approx <- approxSilhouette(reducedDim(sce, "corrected"),
                               clusters=sce$walktrap15)

wp1 <- plotSilBeeswarm(sil.approx)

wp2 <- plotReducedDim(sce, 
                     dimred = "TSNE_corrected", 
                     colour_by="walktrap15", 
                     text_by = "walktrap15")

wp3 <- plotSilGrid(sil.approx)

wp1 + wp2 + wp3

This clustering appears to have generated a set of clusters with slightly better separatedness than the Leiden method with a k of 20.

And again with the louvain clusters:

sil.approx <- approxSilhouette(reducedDim(sce, "corrected"),
                               clusters=sce$louvain15)

lp1 <- plotSilBeeswarm(sil.approx)

lp2 <- plotReducedDim(sce, 
                     dimred = "TSNE_corrected", 
                     colour_by="louvain15", 
                     text_by = "louvain15")

lp3 <- plotSilGrid(sil.approx)

lp1 + lp2 + lp3

There seems to be a greater degree of overlap between these clusters, perhaps more resolution would improve this clustering - we might consider reducing k.

8.2 Modularity to assess clusters quality

As mentioned earlier, the modularity metric is used in evaluating the separatedness of clusters. Some of the clustering algorithms, e.g. Louvain, seek to optimise this for the entire NN graph as part of their cluster detection. Modularity is a ratio between the observed weights of the edges within a cluster versus the expected weights if the edges were randomly distributed between all nodes. Rather than calculating a single modularity value for the whole graph, we can instead calculate a pair-wise modularity value between each pair of clusters using the pairwiseModularity function from the bluster package. For this we need to have the graph from the clustering, so we will rerun the walktrap clustering with k=15 to obtain this. We can plot the resulting ratios on a heatmap. We would expect the highest modularity values to be on the diagonal.

walktrap15 <- clusterCells(sce, 
                           use.dimred = "corrected", 
                           BLUSPARAM = SNNGraphParam(k = 15, 
                                                     cluster.fun = "walktrap"),
                           full = TRUE)
g <- walktrap15$objects$graph
ratio <- pairwiseModularity(g, walktrap15$clusters, as.ratio=TRUE)

hm1 <- pheatmap(log2(ratio+1),
         cluster_rows=FALSE, 
         cluster_cols=FALSE,
         color=colorRampPalette(c("white", "blue"))(100))

We can compare this to the silhouette width grid

wp4 <- ggplotify::as.ggplot(hm1)
wp2 + wp3 + wp4

Largely, this reflects what we saw from the silhouette widths, but also reveals some additional inter-connectedness between other clusters. We can also visualise this as network graph where nodes are clusters and the edge weights are the modularity.

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

set.seed(11001010)
plot(cluster.gr, 
     edge.width=igraph::E(cluster.gr)$weight*5,
     layout=igraph::layout_with_lgl)

8.3 Comparing two sets of clusters

We can assess the concordance between different clustering methods to get a better idea of how they eachtreat the data, e.g. does one cluster from one method equate to just one cluster in the other or is it a combination of different clusters. This may be revealing about the underlying biology. We will use the Jaccard index as measure of concordance between clusters. A value of 1 represents perfect concordance between clusters (i.e. they contain exactly the same cells).

jacc.mat <- linkClustersMatrix(sce$louvain15, sce$walktrap15)
rownames(jacc.mat) <- paste("Louvain", rownames(jacc.mat))
colnames(jacc.mat) <- paste("Walktrap", colnames(jacc.mat))
pheatmap(jacc.mat, color=viridis::viridis(100), cluster_cols=FALSE, cluster_rows=FALSE)

We can see that Louvain clusters 2, 6, 7, and 9 are equivalent to walktrap clusters 14, 7, 12, and 11 respectively. The remaining Louvain clusters are combinations of cells from various walktrap clusters. We may want to look at marker genes for these clusters to assess what these two different views are telling us about the biology.

9 Cluster sweep

As we have seen, there are a number of different parameters we can change to alter the final clustering result - primarily the k used to build the NN graph, the edge weighting method and the clustering algorithm. There is no one gold standard that will fit all data, so, in most cases, it is necessary to assess a number of different clusterings to obtain one that provides a view of the data that suits our biological interpretations. The clusterSweep function allows us to apply a range of different parameters in one go and obtain the clustering for each.

For example, suppose we wish to assess the effect of different values of k on the walktrap clustering. We can parallelize this process to make it faster.

out <- clusterSweep(reducedDim(sce, "corrected"),
                    BLUSPARAM = NNGraphParam(),
                    k = as.integer(c(5, 10, 15, 20, 25)),
                    cluster.fun = "walktrap",
                    BPPARAM=BiocParallel::MulticoreParam(7))

The resulting object is a list containing a DataFrame with the clusters for each combination of the clustering parameters and a corresponding DataFrame showing the parameters used to generate each of these:

out$clusters[,1:4]

## DataFrame with 3500 rows and 4 columns
##      k.5_cluster.fun.walktrap k.10_cluster.fun.walktrap
##                      <factor>                  <factor>
## 1                          19                        9 
## 2                          19                        9 
## 3                          19                        9 
## 4                          27                        20
## 5                          8                         15
## ...                       ...                       ...
## 3496                       9                         5 
## 3497                       15                        1 
## 3498                       20                        14
## 3499                       3                         16
## 3500                       15                        1 
##      k.15_cluster.fun.walktrap k.20_cluster.fun.walktrap
##                       <factor>                  <factor>
## 1                           10                        7 
## 2                           5                         7 
## 3                           5                         7 
## 4                           14                        10
## 5                           10                        9 
## ...                        ...                       ...
## 3496                        2                          3
## 3497                        8                          1
## 3498                        12                         8
## 3499                        9                          5
## 3500                        8                          1

out$parameters

## DataFrame with 5 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.25_cluster.fun.walktrap        25    walktrap

We can then combine this cluster sweep with the metrics for assessing cluster behaviour in order to get a overview of the effects of these parameter changes that may enable us to make some decisions as to which clustering or clusterings we may wish to investigate further.

Here we will just look at the mean silhouette width and the number of clusters.

df <- as.data.frame(out$parameters)

# get the number of clusters
df$num.clusters <- apply(out$clusters, 2, max)

# get the mean silhouette width
getMeanSil <- function(cluster) {
    sil <- approxSilhouette(reducedDim(sce, "corrected"), cluster)
    mean(sil$width)
}
df$silhouette <- map_dbl(as.list(out$clusters), getMeanSil)

nclPlot <- ggplot(df, aes(x = k, y = num.clusters)) + 
                  geom_line(lwd=2)
silPlot <- ggplot(df, aes(x = k, y = silhouette)) + 
                  geom_line(lwd=2)
nclPlot + silPlot

Based on our previous analysis and knowledge of the biology we may feel that 12 clusters represents a good number clusters, and we can see here that k = 25, provides this. Also k = 25 gives us a better silhouette score than lower values of k. On the other hand, perhaps k = 15 provides greater resolution of cell types (more clusters) with only a slight decrease in the silhouette score.

Earlier we looked at the Jaccard index as a means of comparing two different clusterings. We could apply the same method here:

jacc.mat <- linkClustersMatrix(out$clusters$k.15_cluster.fun.walktrap, 
                               out$clusters$k.25_cluster.fun.walktrap)
rownames(jacc.mat) <- paste("Walktrap_15", rownames(jacc.mat))
colnames(jacc.mat) <- paste("Walktrap_25", colnames(jacc.mat))
pheatmap(jacc.mat, 
         color = viridis::viridis(100), 
         cluster_cols = FALSE, 
         cluster_rows = FALSE)

The OSCA book provides some additional methods for comparing different clusterings that can be combined with the cluster sweep results to assess cluster behaviour under different parameters.

In this section, we have just done a sweep changing the k, but it is also possible to combine this with multiple clustering algorithms and multiple edge weightings.

For this exercise, you will rerun clusterSweepwith additional parameters.This time:

test the walktrap, louvain and leiden methods
test k set to 10, 15, 20 and 25

This will test 12 different clusterings.

Hint

This time, as well as setting the k parameter to a vector that contains 10, 15, 20 and 25, you will also need to provide the cluster.fun parameter with a vector containing the names of the three different clustering methods.

Answer

out <- clusterSweep(reducedDim(sce, "corrected"),
                    BLUSPARAM = NNGraphParam(),
                    k = as.integer(c(10, 15, 20, 25)),
                    cluster.fun = c("walktrap", "louvain", "leiden"),
                    BPPARAM=BiocParallel::MulticoreParam(7))

Now plot the number of clusters generated by each clustering and the mean silhouette width. You will need to adjust the plotting parameters to plot a different coloured line for each clustering algorithm (walktap, louvain and leiden).

df <- as.data.frame(out$parameters)

# get the number of clusters
df$num.clusters <- apply(out$clusters, 2, max)

# get the mean silhouette width
getMeanSil <- function(cluster) {
    sil <- approxSilhouette(reducedDim(sce, "corrected"), cluster)
    mean(sil$width)
}
df$silhouette <- map_dbl(as.list(out$clusters), getMeanSil)

nclPlot <- FIXME

silPlot <- FIXME

nclPlot + silPlot

Hint

You will need to add the term aes(colour = cluster.fun) into the geom_line function. This will instruct ggplot2 to plot a different coloured line for each clustering function.

Answer

df <- as.data.frame(out$parameters)

# get the number of clusters
df$num.clusters <- apply(out$clusters, 2, max)

# get the mean silhouette width
getMeanSil <- function(cluster) {
    sil <- approxSilhouette(reducedDim(sce, "corrected"), cluster)
    mean(sil$width)
}
df$silhouette <- map_dbl(as.list(out$clusters), getMeanSil)

nclPlot <- ggplot(df, aes(x = k, y = num.clusters)) + 
                  geom_line(aes(colour=cluster.fun), lwd=2)
silPlot <- ggplot(df, aes(x = k, y = silhouette)) + 
                  geom_line(aes(colour=cluster.fun), lwd=2)
nclPlot + silPlot

Bonus: if you have time, select 1 or 2 of the clusterings, calculate the silhouette widths and examine them further using the beeswarm plot and the silhouette grid.

Hint

We have not yet added the clusters to our sce object, so when running approxSilhouette you will need to use the clustering directly from the clusters object in the out list generated by clusterSweep.

Answer

For example, we could take a look at the leiden with k = 25, which results in 12 clusters and has the highest mean silhouette width.

sil.approx <- approxSilhouette(reducedDim(sce, "corrected"), 
                               clusters=out$clusters$k.25_cluster.fun.leiden)

plotSilBeeswarm(sil.approx)

plotSilGrid(sil.approx)

Finally, we can add all (or a subset) of the clusterings from clusterSweep to our SCE object.

colData(sce) <- cbind(colData(sce), DataFrame(out$clusters))

10 Finalise clustering selection

When you have come to a decision about which clustering to use it is convenient to add it to colData column called “label” using the colLabels function. This means downstream code does not need to be changed should you later decide to switch to a different clustering, you’d simply need to change the contents of the “label” column. This also makes the code easily re-usable for different analyses.

For now we will use the Leiden k=25 clustering.

colLabels(sce) <- sce$k.25_cluster.fun.leiden

11 Expression of known marker genes

If we expect our clusters to represent known cell types for which there are well established marker genes, we can now start to investigate the clusters by plotting in parallel the expression of these genes. This can also help us in assessing if our clustering has satisfactorily partitioned our cells.

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               colour_by = "label", 
               text_by = "label") +
  ggtitle("Leiden k=25 clusters")

Having identified clusters, we now display the level of expression of cell type marker genes to quickly match clusters with cell types. For each marker we will plot its expression on a t-SNE, and show distribution across each cluster on a violin plot.

We will be using gene symbols to identify the marker genes, so we will switch the rownames in the SCE object to be gene symbols. We use the scater function uniquifyFeatureNames to do this as there are a few duplicated gene symbols.

rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ID, rowData(sce)$Symbol)

11.1 B-cells markers

Let’s look at the expression of the B cell marker CD79A.

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               by_exprs_values = "logcounts",
               colour_by = "CD79A",
               text_by = "label")

plotExpression(sce, 
               exprs_values = "logcounts",
               x = "label", 
               colour_by = "label",
               features=c("CD79A"))

A large number of the clusters seem to represent B cells. This is what we expect from a cancer of B cells

11.2 Monocyte markers

We have a very clear cluster of monocytes

plotReducedDim(sce, 
               dimred = "TSNE_corrected",
               by_exprs_values = "logcounts",
               colour_by = "LYZ",
               text_by = "label")

plotExpression(sce, 
               exprs_values = "logcounts",
               x = "label", 
               colour_by = "label",
               features=c("LYZ"))

11.3 Save data

Write SCE object to file.

saveRDS(sce, file="results/Caron_clustering_material.rds")

11.4 Session information

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DT_0.26                     forcats_0.5.2              
##  [3] stringr_1.5.0               dplyr_1.0.10               
##  [5] purrr_1.0.0                 readr_2.1.3                
##  [7] tidyr_1.2.1                 tibble_3.1.8               
##  [9] tidyverse_1.3.2             patchwork_1.1.2            
## [11] pheatmap_1.0.12             igraph_1.3.5               
## [13] cluster_2.1.4               bluster_1.6.0              
## [15] scran_1.24.1                scater_1.24.0              
## [17] ggplot2_3.4.0               scuttle_1.6.3              
## [19] SingleCellExperiment_1.18.1 SummarizedExperiment_1.28.0
## [21] Biobase_2.58.0              GenomicRanges_1.50.1       
## [23] GenomeInfoDb_1.34.3         IRanges_2.32.0             
## [25] S4Vectors_0.36.0            BiocGenerics_0.44.0        
## [27] MatrixGenerics_1.10.0       matrixStats_0.63.0         
## [29] knitr_1.40                 
## 
## loaded via a namespace (and not attached):
##  [1] googledrive_2.0.0         ggbeeswarm_0.6.0         
##  [3] colorspace_2.0-3          ellipsis_0.3.2           
##  [5] XVector_0.38.0            BiocNeighbors_1.14.0     
##  [7] fs_1.5.2                  rstudioapi_0.14          
##  [9] farver_2.1.1              ggrepel_0.9.2            
## [11] fansi_1.0.3               lubridate_1.8.0          
## [13] xml2_1.3.3                codetools_0.2-18         
## [15] sparseMatrixStats_1.8.0   cachem_1.0.6             
## [17] jsonlite_1.8.4            broom_1.0.1              
## [19] dbplyr_2.2.1              compiler_4.2.0           
## [21] httr_1.4.4                dqrng_0.3.0              
## [23] backports_1.4.1           assertthat_0.2.1         
## [25] Matrix_1.5-3              fastmap_1.1.0            
## [27] gargle_1.2.1              limma_3.52.4             
## [29] cli_3.5.0                 BiocSingular_1.12.0      
## [31] htmltools_0.5.4           tools_4.2.0              
## [33] rsvd_1.0.5                gtable_0.3.1             
## [35] glue_1.6.2                GenomeInfoDbData_1.2.9   
## [37] Rcpp_1.0.9                cellranger_1.1.0         
## [39] jquerylib_0.1.4           vctrs_0.5.1              
## [41] DelayedMatrixStats_1.18.2 xfun_0.33                
## [43] rvest_1.0.3               beachmat_2.12.0          
## [45] lifecycle_1.0.3           irlba_2.3.5.1            
## [47] statmod_1.4.37            googlesheets4_1.0.1      
## [49] edgeR_3.38.4              zlibbioc_1.44.0          
## [51] scales_1.2.1              hms_1.1.2                
## [53] parallel_4.2.0            RColorBrewer_1.1-3       
## [55] yaml_2.3.6                gridExtra_2.3            
## [57] yulab.utils_0.0.6         sass_0.4.4               
## [59] stringi_1.7.8             highr_0.9                
## [61] ScaledMatrix_1.4.1        BiocParallel_1.32.1      
## [63] rlang_1.0.6               pkgconfig_2.0.3          
## [65] bitops_1.0-7              evaluate_0.19            
## [67] lattice_0.20-45           labeling_0.4.2           
## [69] htmlwidgets_1.6.1         cowplot_1.1.1            
## [71] tidyselect_1.2.0          magrittr_2.0.3           
## [73] R6_2.5.1                  generics_0.1.3           
## [75] metapod_1.4.0             DelayedArray_0.24.0      
## [77] DBI_1.1.3                 pillar_1.8.1             
## [79] haven_2.5.1               withr_2.5.0              
## [81] RCurl_1.98-1.9            crayon_1.5.2             
## [83] modelr_0.1.9              utf8_1.2.2               
## [85] tzdb_0.3.0                rmarkdown_2.17           
## [87] viridis_0.6.2             locfit_1.5-9.7           
## [89] grid_4.2.0                readxl_1.4.1             
## [91] reprex_2.0.2              digest_0.6.31            
## [93] gridGraphics_0.5-1        munsell_0.5.0            
## [95] ggplotify_0.1.0           beeswarm_0.4.0           
## [97] viridisLite_0.4.1         vipor_0.4.5              
## [99] bslib_0.4.2

Introduction to single-cell RNA-seq analysis

Clustering

Jan 2023

1 Overview

2 Load the data

3 Graph-based clustering overview

3.1 Connecting nodes (cells) based on nearest neighbours

3.2 Weighting the edges

3.3 Grouping nodes (cells) into clusters

3.4 Modularity

3.5 Pros and Cons of graph based clustering

4 Implementation

5 The Walktrap method

6 The Louvain method

7 The Leiden method

8 Assessing cluster behaviour

8.1 Silhouette width

8.2 Modularity to assess clusters quality

8.3 Comparing two sets of clusters

9 Cluster sweep

10 Finalise clustering selection

11 Expression of known marker genes

11.1 B-cells markers

11.2 Monocyte markers

11.3 Save data

11.4 Session information