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

对应原版教程第10章http://bioconductor.org/books/release/OSCA/overview.html
“物以类聚，人以群分” 分群步骤即将基因表达（降维结果）相似的细胞归为同一个群体，往往对应一种特定的细胞类型或者细胞轨迹状态。从一步开始，就可以开始叙述我们的生物学故事了~

源网，侵删~

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

• 细胞分群结果是通过基因表达相似度的计算过程，人为贴上的标签；即本质是一个数学计算问题。
• 如果把分群比作一个显微镜，那么我们可以根据不同的放大倍数(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

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

library(scater) colLabels(sce.pbmc)

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

set.seed(11000) reducedDim(sce.pbmc, "force")

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

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

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

g.5

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

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

g.num

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

g.jaccard

（3）其它划分cluster的算法

• igraph包提供的各种方法

clust.louvain

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

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

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

cluster.gr

3.1 k均值法

（1）算法简介

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

（2）kmeans()函数

set.seed(100) clust.kmeans

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

  

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

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

ncells

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

cent.tree

image.png

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

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

library(cluster) set.seed(110010101) gaps

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

set.seed(0101010) kgraph.clusters

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

• “砍树”分群：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

4 分群结果评价

4.1 cluter间的分离度

(1) 轮廓系数 silhouette width

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

sil

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

# Performing the calculations on the PC coordinates, like before. sil.approx 0, clust, sil.data$other)) sil.data$cluster

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

  

(2) clutering purity

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

pure.pbmc

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

  

4.2 不同clustering 结果的比较

• pheatmap热图形式

tab

• clustree树分布形式

library(clustree) combined

• Rand index指标
  这个指标常用于检测分类模型的预测结果。例如我有2个苹果，2个香蕉，2个芒果；根据模型对这6个水果的分类，使用Rand index指标表示预测结果与真实结果的相似性；
  简单来说，首先A=6个水果所有两两组合的可能性，即(6*5)/(2*1)=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了~