OSCA单细胞数据分析笔记9—Clustering

对应原版教程第10章http://bioconductor.org/books/release/OSCA/overview.html

“物以类聚，人以群分” 分群步骤即将基因表达（降维结果）相似的细胞归为同一个群体，往往对应一种特定的细胞类型或者细胞轨迹状态。从这一步开始，就可以开始叙述我们的生物学故事了~

笔记要点

• 1、clustering是一个显微镜
• 2、基于图聚类的分群
• 3、其它分群算法（k均值与层次聚类）
• 4、分群结果评价

一、clustering是一个显微镜

• 细胞分群结果是通过基因表达相似度的计算过程，人为贴上的标签；即本质是一个数学计算问题。
• 如果把分群比作一个显微镜，那么我们可以根据不同的放大倍数(resolution分辨率)，得到不同的结果。脱离于生物学背景知识，来谈论哪个分群结果是“最佳”的问题，是没有意义。（例如是想观察细胞的primary cell type还是specific cell sub-type。）
• 分群是我们分析scRNA-seq的一个工具，是真正开始结合生物学背景知识的开始。我们可以灵活采用不同的算法、分辨率获得我们“满意”的分群结果。

二、基于图聚类的分群

2.1 算法简介

可简单分为3步

• （1）计算所有细胞两两间的距离(欧几里得距离)，确定每个细胞的Top K nearest neighbors(KNN)；
• （2）根据上述关系，计算细胞(节点)两两间的相关性(边)(节点之间的边的权重)；
• （3）根据保留的KNN网络，划分出内部互相连接关系远高于内外部互相连接关系的cluster。如上分别对应3个问题：

1. 选择多少个最近邻居；
2. 如何度量细胞相关性；
3. 采用什么划分cluster的算法。

2.2 scran包分群实操

• 示例数据

sce.pbmc   #来源参考原教程

• 参数设置
  为每个细胞确定10个最邻近细胞；基于highest average rank of the shared neighbors，计算两两细胞间的关联性；使用igraph包提供的Walktrap算法进行cluster的划分

library(scran)
g <- buildSNNGraph(sce.pbmc, k=10, use.dimred = 'PCA')
clust <- igraph::cluster_walktrap(g)$membership
table(clust)
#clust
# 1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
#205 508 541  56 374 125  46 432 302 867  47 155 166  61  84  16

• 结合t-SNE可视化cluster分布

library(scater)
colLabels(sce.pbmc) <- factor(clust)
plotReducedDim(sce.pbmc, "TSNE", colour_by="label")

值得注意的是：t-SNE二维转换与cluster计算是基于分别PCA降维结果的独立计算的过程。经过上图的可视化呈现，可以起到相互验证的作用。但利用igraph包的系列可视化方式就是基于KNN产生的细胞间关系，进行排布的，如下代码。

set.seed(11000)
reducedDim(sce.pbmc, "force") <- igraph::layout_with_fr(g)
plotReducedDim(sce.pbmc, colour_by="label", dimred="force")

2.3 参数调整对cluster结果的影响

（1） Top n nearest neighbour的选择，即buildSNNGraph()的k=参数设置（*）

• 一般来说，k越大，簇内互相关联程度越紧密，即cluster越大，分得的cluster数越少；k越小，分得的cluster数就越多。

g.5 <- buildSNNGraph(sce.pbmc, k=5, use.dimred = 'PCA')
clust.5 <- igraph::cluster_walktrap(g.5)$membership
table(clust.5)
#  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22 
#523 302 125  45 172 573 249 439 293  95 772 142  38  18  62  38  30  16  15   9  16  13

g.50 <- buildSNNGraph(sce.pbmc, k=50, use.dimred = 'PCA')
clust.50 <- igraph::cluster_walktrap(g.50)$membership
table(clust.50)
#  1   2   3   4   5   6   7   8   9  10 
#869 514 194 478 539 944 138 175  89  45

（2）其它评价细胞间关联性（权重weight）的方法

• Based on the number of nearest neighbors that are shared between two cells

g.num <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="number")

• Based on the Jaccard index of the two sets of neighbors.

g.jaccard <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="jaccard")

（3）其它划分cluster的算法

• igraph包提供的各种方法

clust.louvain <- igraph::cluster_louvain(g)$membership
clust.infomap <- igraph::cluster_infomap(g)$membership
clust.fast <- igraph::cluster_fast_greedy(g)$membership
clust.labprop <- igraph::cluster_label_prop(g)$membership
clust.eigen <- igraph::cluster_leading_eigen(g)$membership

Pipelines involving scran default to rank-based weights followed by Walktrap clustering. In contrast, Seurat uses Jaccard-based weights followed by Louvain clustering.

2.4 评价cluter结果

• 主要目的即是为了不同cluster间的分离度是否足够明显；
• 在笔记最后，会介绍一些通用的评价方法，这里介绍一种专门适用KNN图聚类的评价方法；
• 在之前已经明白了计算细胞间关联性（weight）是图聚类算法的重要步骤（第二步）。设A=某一cluster间任意两两细胞的实际关系（weight）总和；B=某一cluster间两两细胞的预期关系（来自于所有cluster两两细胞间的关系随机分布）。若A-B的差值越大，则说明这个cluster的内联度越显著。
• 为了全面评价同一cluster间的内联度与不同cluster两两间的分离度，使用bluster包的pairwiseModularity()函数进行如上的计算。唯一的区别是用比值ratio代替了差值。

ratio <- pairwiseModularity(g, clust, as.ratio=TRUE)
dim(ratio)

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

如下图，理想情况是对角线的结果最明显，上三角的结果越小越好。

• 可以看到部分cluster间的关联程度还是存在的，我们可以通过igraph进行可视化，进一步查看。

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(11001010)
plot(cluster.gr, edge.width=igraph::E(cluster.gr)$weight*5,
     layout=igraph::layout_with_lgl)

三、其它分群算法

3.1 k均值法

（1）算法简介

• 首先确定想要分成k个cluster；然后在所有细胞中随机抽取k个细胞作为cluster中心点；将中心点以外所有细胞分配到相距最近的中心点，从而形成第一次分群；计算上一次分群的新的中心点，将新的中心点以外所有细胞重新分配到相距最近的中心点；如此往复，直至分群结果稳定下来。
• 根据上述定义，虽然k-means计算速度很快，但存在一些不适合scRNA-seq的地方。例如需要提前指定k的数目，需要多次尝试出比较合适的值；倾向于形成的cluster呈围绕中心点的圆形分布，可能不符合细胞类型的真实分布

（2）kmeans()函数

set.seed(100)
clust.kmeans <- kmeans(reducedDim(sce.pbmc, "PCA"), centers=10)
table(clust.kmeans$cluster)
#  1   2   3   4   5   6   7   8   9  10 
#548  46 408 270 539 199 148 783 163 881

#结合t-SNE进行可视化
colLabels(sce.pbmc) <- factor(clust.kmeans$cluster)
plotReducedDim(sce.pbmc, "TSNE", colour_by="label")

• 如上图 k均值的分群结果在t-SNE分布还是较为理想的。

（3）评价k-means分群结果

• 首先计算每个cluster内所有点到中心点的距离平方和(WCSS, within-cluster sum of squars)
• 然后计算每个cluster内细胞据中心点的平均距离，最后开根号的结果（RMSD，root-mean-squared-deciation）作为该cluster的评价指标。
• 若cluster的RMSD越小，表明该群的分布越紧凑，即效果越好

ncells <- tabulate(clust.kmeans$cluster) #统计1-10群细胞数量
tab <- data.frame(wcss=clust.kmeans$withinss, ncells=ncells)
tab$rms <- sqrt(tab$wcss/tab$ncells)
tab
#       wcss      ncells   rms
# 1  20626.970    548  6.135182
# 2   6530.806     46 11.915286
# 3   9899.650    408  4.925835
# 4   6429.640    270  4.879906
# 5  12413.267    539  4.798977
# 6  13467.078    199  8.226406
# 7   4718.144    148  5.646180
# 8  27010.471    783  5.873341
# 9   7703.558    163  6.874670
# 10  9469.524    881  3.278507

• 可以顺便可视化下基于cluster中心点的层次聚类关系（关于层析聚类算法，之后会说一下）

cent.tree <- hclust(dist(clust.kmeans$centers), "ward.D2")
plot(cent.tree)

（4）关于中心点(centroid)数量k的选择

• A：关于k的最佳取值，我觉得针对scRNA-seq还是要多尝试不同的分群结果。教程也介绍了一种可以用来选择最佳k值的思路：计算每次分群结果的gap统计量，一般越高越好。可使用cluster包的相关函数，相关代码如下，具体原理可参考原教程。

library(cluster)
set.seed(110010101)
gaps <- clusGap(reducedDim(sce.pbmc, "PCA"), kmeans, K.max=20) #耗时
best.k <- maxSE(gaps$Tab[,"gap"], gaps$Tab[,"SE.sim"])
best.k
# 8

• B：需要提前k值的k均值法可以设定一个较大的值，达到图聚类法难以达到的分辨率。虽然进行生物水平的可解释性不高，但可实现从所有细胞中，抽取k个有代表性表达情况的细胞的目的，用于某些特定的分析场景。
• 例如 clusterRows {bluster}提供一种联合图聚类与k-均值聚类的方法，可明显的优势是相对于单纯图聚类大大提高了分析速度。
• 简单举例来说：首先使用k均值法，获得所有细胞的k个代表性细胞（一般取较大的值，如1000等）。然后使用图聚类法对这1000个中心点细胞矩形聚类。最后把归于相应中心点的其余细胞分到中心点所在的图聚类结果里。
• 相关代码如下

set.seed(0101010)
kgraph.clusters <- clusterRows(reducedDim(sce.pbmc, "PCA"),
    TwoStepParam(
        first=KmeansParam(centers=1000),  #1000个中心点
        second=NNGraphParam(k=5)  #5个最近邻居
    )
)
table(kgraph.clusters)
#  1   2   3   4   5   6   7   8   9  10  11  12 
#191 854 506 541 541 892  46 120  29 132  47  86

3.2 层次聚类法

（1）算法简介

• 第一步将每个细胞(n)看做一个单独的cluster，然后计算所有cluter两两间的“距离”，将相距最近的两个cluster归为一个cluster；第二步再次计算（n-1）个cluster两两间的“距离”，将相距最近的两个cluster归为一个cluster；如此往复循坏，直至归为最后一个最终的cluster。不同分群的结果即取决于距离阈值的切分（如上图所示，阈值=4）
• 该过程中最关键步骤是如何度量cluster间的距离，尤其是对于从第二步开始，不同cluster的细胞数是不一致的。相关算法也有很多，例如

complete linkage aims to merge clusters with the smallest maximum distance between their elements, while Ward’s method aims to minimize the increase in within-cluster variance.

• 层次聚类最直接的结果就是得到如上图所示的系统发生树(dendrogram)，从某种角度代表了细胞从上至上或者从上至下的关系。但另一方面，层次聚类法往往不适于动辄成千上万的细胞的计算，相对来说适合一些小的数据集；或者某一特定细胞群；或者结合k-均值的结果。

（2）hclust()函数实操

sce.416b
#含有185个细胞的sce子集。获取方式见原教程

dist.416b <- dist(reducedDim(sce.416b, "PCA"))
tree.416b <- hclust(dist.416b, "ward.D2") #Ward’s method

# Making a prettier dendrogram.
library(dendextend)
tree.416b$labels <- seq_along(tree.416b$labels)
dend <- as.dendrogram(tree.416b, hang=0.1)

combined.fac <- paste0(sce.416b$block, ".", 
                       sub(" .*", "", sce.416b$phenotype))
labels_colors(dend) <- c(
  `20160113.wild`="blue",
  `20160113.induced`="red",
  `20160325.wild`="dodgerblue",
  `20160325.induced`="salmon"
)[combined.fac][order.dendrogram(dend)]

plot(dend)

• “砍树”分群：cutree()函数，有两个参数，可分别设置。k=表示想分成多少个cluster；h=表示基于什么距离阈值(上图的纵坐标)分隔；

cutree(dend,k=4)
cutree(dend,h=200)

• “智能砍树”：dynamicTreeCut包的cutreeDynamic()函数，可自动寻找一个最佳的阈值，进行分群

library(dynamicTreeCut)

# minClusterSize needs to be turned down for small datasets.
# deepSplit controls the resolution of the partitioning.
clust.416b <- cutreeDynamic(tree.416b, distM=as.matrix(dist.416b),
                            minClusterSize=10, deepSplit=1)
table(clust.416b)
# 1  2  3  4 
#78 69 24 14
labels_colors(dend) <- clust.416b[order.dendrogram(dend)]
plot(dend)

4 分群结果评价

4.1 cluter间的分离度

(1) 轮廓系数 silhouette width

• 方法简介
  对于某一个cluster的某一个细胞来说，设A=该细胞到所在cluster所有细胞的平均距离；B=该细胞到其它所有cluster里的所有细胞的平均值的最小值(即与该细胞相距最近的其它cluster)。所以对于该细胞的轮廓系数=(B-A)/max{A,B}。该值越接近1，表明分群越合理；负值则表示该细胞可能是个“叛徒”
• silhouette()函数计算

sil <- silhouette(clust.416b, dist = dist.416b)
plot(sil)

• 针对大的scRNA-seq数据集，推荐使用approxSilhouette()函数采用近似的方法计算所有细胞的轮廓系数。

# Performing the calculations on the PC coordinates, like before.
sil.approx <- approxSilhouette(reducedDim(sce.pbmc, "PCA"), clusters=clust)

sil.data <- as.data.frame(sil.approx)
sil.data$closest <- factor(ifelse(sil.data$width > 0, clust, sil.data$other))
sil.data$cluster <- factor(clust)

ggplot(sil.data, aes(x=cluster, y=width, colour=closest)) + 
  ggbeeswarm::geom_quasirandom(method="smiley")

• 如下图，横坐标为既定的分群结果，每个点代表一个细胞，点的颜色即根据所有cluster所属细胞的轮廓系数标注的分群评价。如果轮廓系数为负值，即归为相距最近的其它类。

(2) clutering purity

• 简单来说就是：在既定的分群结果里，对于一个特定cluster的每个细胞来说，其近邻细胞都为该cluster的比例；比例越接近1越好。
• 若一个cluster里的所有细胞的purity的中位数大于0.9，那么表明分群结果理想；
• neighborPurity函数

pure.pbmc <- neighborPurity(reducedDim(sce.pbmc, "PCA"), clusters=clust)

pure.data <- as.data.frame(pure.pbmc)
pure.data$maximum <- factor(pure.data$maximum)
pure.data$cluster <- factor(clust)

ggplot(pure.data, aes(x=cluster, y=purity, colour=maximum)) +
    ggbeeswarm::geom_quasirandom(method="smiley")

• 如下图，横坐标为既定的分群结果，每个点代表一个细胞，点的颜色即根据每个细胞的近邻细胞中所属cluster占比最多的所决定。

4.2 不同clustering 结果的比较

• pheatmap热图形式

tab <- table(Walktrap=clust, Louvain=clust.louvain)   #2.3
#行为Walktrap， 列为Louvain
tab <- tab/rowSums(tab) #Louvain结果相对于Walktrap的分布比例(按行看)
pheatmap(tab, color=viridis::viridis(100), cluster_cols=FALSE, cluster_rows=FALSE)

• clustree树分布形式

library(clustree)
combined <- cbind(k.50=clust.50, k.10=clust, k.5=clust.5)
clustree(combined, prefix="k.", edge_arrow=FALSE)

• Rand index指标
  这个指标常用于检测分类模型的预测结果。例如我有2个苹果，2个香蕉，2个芒果；根据模型对这6个水果的分类，使用Rand index指标表示预测结果与真实结果的相似性；
  简单来说，首先A=6个水果所有两两组合的可能性，即(6x5)/(2x1)=15：而B=所有组合的真实分布与预测分布要么均归为1类（苹果1与苹果2预测归为1类），要么均归为不同类(苹果1与香蕉1预测归为不同类)的次数。则Rand index=B/A，越接近1，越好。

pairwiseRand(clust, clust.5, mode="index")
# 0.7796

最后关于subclustering，即在clustering结果基础上，对某一个cluster进一步分群，是提高细胞亚群分辨率的一种方法。步骤算法基本同上，但可能也会遇到一些坑，详见原教程最后，这里就不展开介绍了。

---

以上学习了scRNA-seq分群涉及到的一些知识点。下一步就是找出每个cluster的marker gene了~

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注