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

生信技能树jimmy

发布于 2021-07-02 18:25:28

1.8K1

发布于 2021-07-02 18:25:28

文章被收录于专栏：单细胞天地单细胞天地

对应原版教程第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了~

往期回顾

NC单细胞文章复现(三)：复杂热图

scPhere——用地球仪来展示降维结果

2021第一期生信入门微信群答疑精选200题

开机，写bug

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

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2021-06-02，如有侵权请联系 cloudcommunity@tencent.com 删除

聚类算法

编程算法

linux

本文分享自单细胞天地微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

聚类算法

编程算法

linux

登录后参与评论

0 条评论

热度

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

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

一、clustering是一个显微镜

二、基于图聚类的分群

2.1 算法简介

2.2 scran包分群实操

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

2.4 评价cluter结果

三、其它分群算法

3.1 k均值法

（1）算法简介

（2）`kmeans()`函数

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

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

3.2 层次聚类法

（1）算法简介

（2）`hclust()`函数实操

4 分群结果评价

4.1 cluter间的分离度

(1) 轮廓系数 silhouette width

(2) clutering purity

4.2 不同clustering 结果的比较

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

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

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

一、clustering是一个显微镜

二、基于图聚类的分群

2.1 算法简介

2.2 scran包分群实操

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

2.4 评价cluter结果

三、其它分群算法

3.1 k均值法

（1）算法简介

（2）kmeans()函数

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

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

3.2 层次聚类法

（1）算法简介

（2）hclust()函数实操

4 分群结果评价

4.1 cluter间的分离度

(1) 轮廓系数 silhouette width

(2) clutering purity

4.2 不同clustering 结果的比较

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

（2）`kmeans()`函数

（2）`hclust()`函数实操