对应原版教程第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个问题：
- 选择多少个最近邻居；
- 如何度量细胞相关性；
- 采用什么划分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")
代码语言：javascript复制值得注意的是：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
scrandefault to rank-based weights followed by Walktrap clustering. In contrast,Seuratuses 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()函数
代码语言：javascript复制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()函数实操
代码语言：javascript复制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了~
往期回顾
NC单细胞文章复现(三)：复杂热图
scPhere——用地球仪来展示降维结果
2021第一期生信入门微信群答疑精选200题
开机，写bug
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
- 生信爆款入门-2021第4期
- 数据挖掘（GEO,TCGA,单细胞）2021第4期
- 明码标价之共享96线程384G内存服务器
