|
嘿，大家好我又来了！回忆一下，在前两期的文章中我已经给大家介绍了常见的聚类如层次聚类、k-means等，以及为单细胞专门开发的SC3聚类。今天带来本系列的最后一个聚类方法——谱聚类（Spectral Clustering）的介绍。听这个名字大家可能会觉得很陌生，但是如果你有过单细胞的分析经历，可能会在无意中使用过它。没错，它就是大名鼎鼎的R包Seurat自带的默认聚类方法。简单来说，Seurat中的谱聚类基于共享最近邻图（SNN）和模块化优化的聚类算法来识别细胞簇。它首先计算k最近邻(k-nearest neighbors)并构造SNN图， 然后优化模块化功能以确定具体类群。听起来好像有点难懂，没办法，高大上的方法背后一般都是深奥的理论作为支撑，并且这个算法最早居然在2013年发表在物理学的期刊The European Physical Journal B上。我会尽量简化其中所涉及到的数学相关知识，更多地传达方法背后生物学意义。第一个难点是谱聚类是一种基于图论方法，如果你傻傻地分不清图论的“图”和图像的“图”，建议先去学习一下基本的离散数学。谱聚类是从图论中演化出来的算法，它的主要思想是把所有的细胞看为高维空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低；距离较近的两个点之间的边权重值较高。然后通过对所有数据点组成的图进行“切图”，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。谱聚类中的图采用的是无向有权图，也就是说，每个细胞（图中的点）之间的连接是没有方向性的。但是每条边上是有权重信息的。类比一下地图，两地之间的路是没有方向（双向的），但是地点与地点之间路的远近距离是不同的。对于一个无向有权图G，一般用点的集合V和边的集合E来描述，即为G(V,E)。其中V即为我们数据集里面所有的点(v1,v2,...vn)的集合。对于V中的任意两个点，定义权重wij为点vi和点vj之间的权重。由于是无向图，所以wij=wji如果两点之间没有边连接，那么wij=0。OK，如果你看到这都能懂，那么接下来也没啥问题了。第二个难点如何把我们的单细胞的数据构建成为无向有权图，需要明确的是，我们手中最初只有细胞的表达量矩阵，这个每个细胞中各个基因的表达情况，属于“节点信息”。那么哪些细胞之间需要有边相连，边的权重又该如何计算呢？这里引入了邻接矩阵(W)的概念，我们首先定义一种计算细胞之间距离度量的方法，来计算每条边之间的权重wij，又称相似性sij，然后获得整个的邻接矩阵W。构建邻接矩阵W的常用方法有三类：ϵ-邻近法，K邻近法和全连接法。Seurat中选用的是k-近邻法，那我们也就把k-近邻法作为例子讲解。首先利用KNN算法遍历所有的样本点，取每个细胞最近的k个细胞作为近邻，然后设置只有与此细胞距离最近的k个点之间的wij>0，其他全设置为0。但是这种方法会造成重构之后的邻接矩阵W非对称，我们后面的算法需要对称邻接矩阵。为了解决这种问题，一般采取下面两种方法之一(1)只要第一个点在另一个点的k近邻中即可保留sij（2）必须两个点互为k近邻才可保留sij。由此我们就能用计算出的权重矩阵构建出如下所示的图。 
然后是度矩阵(D)的概念：对于图中的任意一个点vi，它的度di定义为和它相连的所有边的权重之和，即
从每个点度的定义，我们可以得到一个n x n的度矩阵D,它是一个对角矩阵，只有主对角线有值，对应第i行的第i个点的度的值
好了，最后一个重要的矩阵：拉普拉斯矩阵（L）登场。拉普拉斯矩阵定义很简单：L = D -W , W是上面所说的邻接矩阵，D是度矩阵。但是拉普拉斯矩阵有着很重要的性质。终于到了“切图”的环节，实际就是将图切成相互没有连接的k个子图。目标就是让子图内的点权重和高，子图间的点权重和低。由于这里的求解过程十分复杂全是公式，具体的过程就不讲了，实际就是一个最优化问题。可以理解为通过D，L这两个矩阵所给的信息，在图的部分边进行切分，让完整的图变成许多子图，每一个子图就构成了一个类群，如下图所示
当然，如果你只想在Seurat中简单的用一下，直接用FindClusters()这个函数就好了。就是简单的一句话’pbmc <- FindClusters(pbmc, resolution = int)’。最常调整参数有resolution和pc.use，pc.use决定哪些维度用来进行聚类，如果之前做过PCA，会选一般前10-20维度。resolution参数决定类群的“粒度”，官网建议将此参数设置在0.4-1.2之间通常可为包含大约3K个cell的数据集返回良好的聚类结果。对于较大的数据集，最佳resolution值通常会增加。值得一提的是，FindClusters()函数是不需要输入希望划分类群数量的，由算法自动划分。调整resolution会得到不同的类群的结果，越高细胞类群数量越多，想获得更多的类群，可以尝试高于1.0的值。最后加上一句umap或者t-SNE的绘图命令就能可视化聚类结果了。下图是我用Seurat自带的pbmc数据的聚类结果，可以发现效果还是非常不错的。 
好啦，以上就是谱聚类的所有内容了。如果你不太懂算法原理，可以不用过于纠结，能在R语言中灵活运用才是解决生物信息问题的关键。
|
