3.5 聚类
刘小泽写于2020.6.30 + 7.2
1 前言
聚类是一种无监督的学习过程，在scRNA分析中，根据表达谱相似性对不同细胞进行分组。
在后面根据marker基因对各个cluster进行注释后，cluster就可以代表不同的细胞类型或状态。因此，聚类是将枯燥的数据与感兴趣生物学意义联系起来的一个桥梁。
不过这里需要明确一下，聚类分的群和细胞类型还不是一回事
分群是根据经验计算的结果
细胞类型有真正的生物学意义
至于这种问题：“我应该设置多少cluster数量？”，则是个伪命题。因为我们可以根据需要，使用不同方法，得到不同数量的clusters，每个cluster都是高维空间数据的一部分。
不过倒是可以问：“我们的clusters与细胞类型之间的拟合程度如何？”。不过这个问题很难回答，因为取决于自己的分析目标。有的人对于能得到主要的细胞类型就很满意，而有的人还想再进一步，得到细胞亚群；还有的，更想基于细胞亚群，达到不同细胞状态的分辨率（如代谢活性、胁迫反应等）。另外，两种截然不同的分群结果，可能都具有生物学意义，只是关注点不同而已。
对于分群聚类的操作，更像是显微镜的操作，可以通过调整参数（使用不同的镜头组合）来获得不同的分辨率；另外还可以探索多种聚类方法，看到数据的不同面。
数据准备之PBMC
具体的操作在上一章也提到过
# 下载
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
"cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
untar(raw.path, exdir=file.path(tempdir(), "pbmc4k"))
# 读取
suppressMessages(library(DropletUtils))
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
> dim(sce.pbmc)
[1] 33694 737280
# 基因注释之ID整合
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
# 基因注释之添加位置信息
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID,
column="SEQNAME", keytype="GENEID")
# 空液滴检测
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
# 质控
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]
# 归一化（去卷积）
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
# 利用数据分布来表示技术噪音，并挑选HVGs
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)
# 降维（三种方法）
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc, subset.row=top.pbmc, technical=dec.pbmc)
set.seed(100000)
sce.pbmc <- runTSNE(sce.pbmc, dimred="PCA")
set.seed(1000000)
sce.pbmc <- runUMAP(sce.pbmc, dimred="PCA")
# 看下结果
sce.pbmc
## class: SingleCellExperiment
## dim: 33694 3922
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
## colData names(3): Sample Barcode sizeFactor
## reducedDimNames(3): PCA TSNE UMAP
## altExpNames(0):补充知识点之聚类与分类
分类：类别是已知的，通过对已知分类的数据进行训练和学习，找到这些不同类的特征，再对未分类的数据进行分类。属于有监督学习。
聚类：事先不知道数据会分为几类，通过聚类分析将数据聚合成几个群体。聚类不需要对数据进行训练和学习。属于无监督学习。
2 基于图形的聚类 | Graph-based clustering
2.1 是什么？
这种聚类方法在Seurat中使用最流行。最基础的想法是：我们首先构建一个图，其中每个节点都是一个细胞，它与高维空间中最邻近的细胞相连。连线基于细胞之间的相似性计算权重，权重越高，表示细胞间关系更密切。
如果一群细胞之间的权重高于另一群细胞，那么这一群细胞就被当做一个群体 “communiity”
这个方法最大的优点就是：可缩放性。只需要根据 k-nearest neighbor（KNN）进行搜索，自己去寻找去拓展。相比于k-means、 Gaussian mixture models等方法，不需要预先对cluster形状以及cluster中细胞分布做一个估计。
这种方法使得每个细胞都被强制连接到一定数量的相邻细胞，这减少了仅由一两个离群细胞组成的无意义群体的风险
不过这个方法也有缺点：
缺点一：
它最后只保留了邻近细胞之间的联系，除此以外没有保存。
缺点二：
另外正是由于它主要关注相邻的权重，导致一些噪音也可能聚集在一起，形成一个community，从而影响聚类的分辨率，可能导致数据异质性的过度解读
2.2 怎么做？
在使用时，需要考虑几个问题：
构建这个图需要设置几个邻居？
用什么方法确定边的权重？
用什么检测community的方法来定义cluster？
下面的例子中，就定义了使用一个细胞与10个邻居的关系，进而构建shared nearest neighbor graph（SNNG），如果两个细胞的邻居中有一个是共享的，则通过一条边连接，之后再根据highest average rank of the shared neighbors定义边的权重(Xu and Su 2015)
首先使用scran包的buildSNNGraph()函数，并且使用PCA的前几个PCs；之后利用 igraph包中的Walktrap方法去鉴定clusters
之后可以把cluster信息放回到SingleCellExperiment对象中，保存成一个因子
从上面的计算方法可以看出，有一个很重要的参数是k ，含义是：the number of nearest neighbors used to construct the graph。如果k设置越大，得到的图之间联通程度越高，cluster也越大。因此这个参数也是可以不断尝试的
这个结果可以利用force-directed的布局出图，它与t-SNE、UMAP的降维结果都是紧密相关的
2.3 其他的参数
关于如何计算权重，在buildSNNGraph()函数中可以设置：
type="number"：根据近邻数量计算type="jaccard"：根据Jaccard index of the two sets of neighbors计算
上面得到的各种结果，都是来自igraph包的graph对象，因此可以用于igraph包的各种community检测方法：
最后就可以比较两种不同的聚类方法差异
可以看出，Infomap比Walktrap得到了更多的cluster（比较精细），而fast-greedy结果相对较少（比较粗糙）
可以通过cut_at()指定分群的数量
scran默认使用rank-based的权重计算方法，Seurat使用Jaccard-based的方法（后面辅助Louvain聚类）
2.4 分群的评价
当使用graph-based方法时，模块化（modularity）是一个评价cluster的指标。modularity值越大，表示cluster内部的细胞与细胞之间连线（edge）数最多，内部分散效果越好，从而避免了不同cluster之间邻近的细胞形成连线。
之前通过cluster_walktrap的方法得到了18个cluster，下面就看看这个方法做的如何：会用到一个函数clusterModularity() ，并且最好设置一个参数as.ratio=TRUE ，表示两两cluster之间比较，观察到的权重与预期权重之比，这个比值（ratio）越大，表示正在比较的两个cluster之间的细胞间连线数越多，联系越密切。
下面看一个基于ratio的热图
从图中可以看到，深色区域（ratio较大的值）基本都集中在对角线，也就是同一个cluter之中。因此，同一个cluster中会存在最多的细胞间连线数量，而不是与其他的cluster，这不也正是设置cluster的目的么？一个cluster中的细胞的相关性就是应该高于其他cluster中的细胞
对角线说明了大部分的cluster划分合理，当然也不排除有一些对角线以外的第二高ratio值，表示那个cluster与其他的cluster之间也有细胞间连线。
下面看一个基于ratio的点线图
就像这种，和之前使用的graph不同，之前的点（node）是细胞，这里的node是cluster，看像不像分化轨迹做的”拟时序分析“？其实就是把细胞聚类后，再聚一次，找到它们的前后相互关系
这个图是这么做的：
3 k-means聚类
3.1 是什么
顾名思义，就是通过不断迭代找出k 个中心(centroid)，然后将各个细胞分配至最近的中心点。看似晦涩，来看看k-means 聚类解读
这个方法算法简单、易于实现，具有速度优势。但是存在几个严重的缺点，可能会导致分配的cluster不清不楚：
k值需要提前定义。假如真实有10种细胞类型，但我现在只定义k=9，那么有两个细胞类型就会硬生生地”合二为一“。但其他方法，例如上面的基于图形的，即使分辨率设置的很低（意思就是看的很模糊），也是能看出这两种不是同一种类型
根据它的计算方法，需要不断重复多运行几次，保证得到的cluster结果是稳定的
既然k-means是找k个中心，那么有中心就会有区域半径，就会限定数据的范围。但实际上，很多分组并没有常规的数据大小和结构，也就是说，一组数据不会这么”乖乖地“落在我们画好的范围中
看了上面的优缺点后，k-means依然成为了最好用的样本数据压缩方法之一。另外它可以不受细胞密度的影响，保证了即使是数量最多的细胞类型也不会主导下游分析
3.2 怎么做？
基础包中就有函数kmeans() ，我们这里基于PCA的结果，指定目标中心数为10，另外也是需要设置随机种子
当然，可以在一开始就设置k大一点（多找一些”核心人物“），避免发生细胞亚群重叠的情况。
3.3 分群的评价
之前了解了k-means的计算方法（计算每个小弟与各个中心大佬的距离，离哪个近就跟哪个），那么就可以用within-cluster sum of squares (WCSS) 方法对每个亚群进行评价。顾名思义，这个方法就是寻找最近的距离。
另外，基于WCSS又可以计算root-mean-squared deviation (RMSD)，反映的是偏离平均位置的程度。当然，它的计算也很简单，就是用每个亚群的WCSS值除以该亚群的细胞数量然后再开方。
RMSD值低的好。如果RMSD值低，就表示和其他亚群的混杂程度更小，更像是一个纯粹的亚群（例如rms值较小的8、10群，它们当中混杂的其他亚群数量很少）。反观rms值大的群，如16、19群，都混杂了其他亚群的细胞（尤其是19群，乍一看上去这个群颜色很单一，但放大看，就会发现其中除了蓝色还有灰色）
再结合上一张t-SNE的图，看到rms指标与t-SNE图上反映的cluster大小基本没有相关性（例如第3群比第12群细胞数量少，同样第12群rms大；第3群比第1群细胞数量少，却是第3群rms大）。这里也侧面反映了之前介绍t-SNE时说的：t-SNE就是一个可视化的工具，不要试图量化t-SNE结果上的cluster。上面的点不能说明什么问题，t -SNE的算法会让稠密的点变膨胀，让原本稀疏的点变紧凑，不可以用cluster的大小来衡量亚群的异质性。另外它并没有义务帮我们保留cluster的相对位置，因此我们不能根据点的位置远近来确定cluster之间的相似程度，每运行一次点的位置都是变动的。
那么可能会想，如果不通过t-SNE的结果来判断cluster之间的关系，那怎么看呢？
可以通过层次聚类
同样看到，19群格格不入，与之前t-SNE的可视化结果，以及rms结果都相吻合
3.4 k-means的更高级用法
之前提到过，k-means可以将邻近的细胞靠一个”中心大佬“来体现，用一个中心(centroid)来表示其余多个点，达到数据压缩的目的。因此，它可以用到其他更复杂的聚类方法中，作为铺垫先压缩一下数据。
例如scran包中有一个函数clusterSNNGraph() ，就是整合了k-means和graph-based
先利用k-means得到一些关键的中心点（centroid）
之后将中心点进行graph-based聚类
结果再把每个中心点（centroid）周围的”小弟们“接过来，放在中心点周围
这种方法相比于单纯的graph-based聚类，速度大大提升。（回顾一下之前说的graph-based方法）因为不用去为每个细胞找邻居，从而构建一个图。另外之前说graph-based这种方法使得每个细胞都被强制连接到一定数量的相邻细胞，但k-means的”核心大佬“加入减轻了这个问题。
注意到，上面的代码使用了一个参数kmeans.centers， 指的是k-means获得cluster的数量，它的设置取决于是要处理速度还是要数据还原度。设置越大越能表示细胞真实分布，但也为第二步聚类造成了计算量激增。
4 层次聚类
4.1 是什么
这是一种历史悠久的聚类方法，最后会给出样本的聚类树（dendrogram）。思路是这样的：贪婪地将样本聚集成簇，然后将这些簇聚集成更大的簇，依此类推，直到所有样本都属于一个簇为止。
层次聚类包含的算法也有很多种，差别在于如何层层聚集。例如complete linkage方法（hclust的默认方法）是保证合并簇之间样本距离最小， Ward’s方法是保证一个簇内部在聚合的同时方差增加的最少，也就是将每次合并时的信息损失最小化。
在实际使用中，层次聚类运行很慢，因此只适用于小型scRNA数据。因为需要计算一个细胞间的距离矩阵，如果动辄上千细胞，那么这个计算量会很大。
4.2 怎么做
之前用的PBMC数据集中细胞数太多，于是切换到sce.416b这个数据集，毕竟一二百个细胞还是可以算的
首先利用前几个PCs计算细胞间的距离
然后使用Ward‘s方法进行聚类
准备画复杂一点的聚类树
图中可以看到，处理前后确实形成了红、蓝分界
另外更推荐使用Ward’s方法，它在合并时受到的cluster差异影响最小（毕竟人家的出发点就是每次合并时的信息损失最小化，所以如果差异很大，那我就不合并呗）。
为了得到更清楚地分群结果，可以对树进行”修建“
最简单的是
cutree()，例如cutree(hc, 4)手动切分数量复杂一点但更好用的是
cutreeDynamic()，来自dynamicTreeCut这个包
一般这个结果与t-SNE的结果是吻合的
不过这里注意到，t-SNE中cluster2拆分开了，值得怀疑：是真的聚类出了问题，还是t-SNE本身的显示问题，毕竟我们知道t-SNE有时真的会把一个亚群拆开画。于是继续下面👇的探索评价
4.3 评价
可以借助”轮廓图“（silhouette width）检查分群的质量。会对每个细胞都计算一个silhouette width值，如果一个细胞的width值为正并且越大，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，分群效果越好；如果width为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。
图中可以反映的问题：
这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。
这里看到cluster2都是正值并且横坐标width数值还是最大的，因此说明了分群没有问题，之前t-SNE的结果的确是它本身的问题
如果发现width值较小，表示分群结果一般，还有可能是分群多度，本来属于一个群的，又被拆分成小群
利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果。因为即使图上看似合理的分群，可能实际上也不会得到更多的生物信息
5 评价分群的稳定性
什么叫做分群的稳定性？
就是分群之前操作的小小波动，也不会影响到最后的分群结果；高稳定性的分群结果其实也方便了别人进行重复。scran利用自助法（bootstrapping）进行评价
在统计学中，自助法（Bootstrap Method，Bootstrapping，或自助抽样法）是一种从给定训练集中有放回的均匀抽样，也就是说，每当选中一个样本，它等可能地被再次选中并被再次添加到训练集中。 自助法由Bradley Efron于1979年在《Annals of Statistics》上发表。
细胞有放回的抽样，得到一个”自助重复“数据集，然后用他进行聚类，看看是否能得到相同的亚群
结果返回一个coassignment的矩阵，表示随机从clusterX与clusterY中挑一个细胞，这两个细胞在经历”自助法“判断后，依然落在同一个cluster的概率。我们希望落在对角线上的概率更高（因为对角线表示自己和自己比较，说明多次重复后细胞还是属于这个cluster）
需要注意的是
这个自助法（Bootstrapping）是一个通用的评价cluster稳定性的方法，适用于各种聚类算法
高稳定性的亚群不一定是高质量的，因为如果一个亚群质量一直很差，它也很稳定
6 数据集取小再聚类 | Subclustering
首先通过常规操作得到分群结果，找到某些感兴趣的亚群，然后在某一个cluster中再一次进行特征基因挑选、聚类操作，从而增加分析的精确度。
举个例子：
首先利用整个PBMC数据集，根据T细胞的marker基因筛出可能是T细胞的亚群，然后推测其中某个亚群可能是记忆T细胞
现在推测cluster6可能是记忆T细胞
得到亚群cluster6的CD4+和CD8+ 的亚亚群
这种数据集取小再聚类的方法，可以简化亚群的解读。只需要在某一个亚群的背景下继续挖掘，就可能会获得重要信息，一层层递进。比如上面的思路：一堆细胞=》分18群=》鉴定T细胞的亚群=》取出某个T细胞亚群再去鉴定记忆T细胞
最后更新于
这有帮助吗？