单细胞聚类方法调研

Review of Single Cell Clustering

Posted by Cheereus on February 22, 2020

单细胞聚类方法调研

背景略。

SNN-Cliq

首先构建 SNN 图。

首先使用欧氏距离(其他距离计算方法亦可)计算成对数据点之间的相似矩阵。

然后,对于每个数据点 $x_i$,根据相似矩阵列出 $k$ 最近邻居表,$x_i$ 本身作为其中的第一个。

构建 SNN 图:对于一对至少有一个共享 KNN 的数据点 $x_i,x_j$,指定一条边 $e(x_i,x_j)$,这条边的权值定义为 $k$ 与共同 KNN 的最高平均排名的差:

其中 $k$ 为最近邻表的大小,$rank(v,x_i)$ 表示节点 $v$ 在 $x_i$ 的最近邻表 $NN(x_i)$ 中的排名(排名以第 1 名为最近)。

共享最近邻的排名相比于其他普通相似度计算方法的优点是在高维空间中依然具有其意义。而且 SNN 图通常是稀疏的,有利于应用到大数据集上。

然后识别 SNN 图中的类。

最大团是一个完整(完全连接)的子图且不包括在更大的团中。尽管枚举出所有的最大团是一个 NP 难的问题,但每个节点所相关的最大团可以用启发式的方法有效地找到。

然而由于 SNN 图的一般稀疏性,团是罕见的。我们寻找伪团来取而代之,伪团足够密集但不一定完整(完全连接)。寻找与每个节点相关的伪团使用贪心算法来实现:

(1) 对于由节点 $v$ 构造的子图 $S$ (包括 $v$ 及其所有邻居节点和相关的边)中的每一个节点 $s$,计算一个局部的度量 $d$ 为从 $S$ 中其他节点入射到 $s$ 的边数。

(2) 在 $S$ 中选取具有最小 $d$ 值 $d_i$ 的节点 $s_i$,如果 $\frac{d_i}{\vert S\vert}<r$($\vert S\vert$ 为 $S$ 中的节点总数,$r$ 为预定义的阈值,$0<r\leq 1$),将其从 $S$ 中去除。

(3) 更新 $d$ 值,并重复步 (2) 直到没有节点会被去除,如果此时 $\vert S \vert \geq 3$ 即剩下的 $S$ 中至少有三个节点,则称此时的 $S$ 为伪集团。

(4) 找出所有的伪团后,要对结果进行去冗余,即删除所有完全包含在其他伪团中的结果。

参数 $r$ 定义了结果中各伪团的连接度,高的 $r$ 值表示各子图更为紧凑,经过测试,发现 $r=0.7$ 时有比较好的效果。

然后对于任意两个子图 $S_i$ 和 $S_j$,定义一个重叠指数 $O_{i,j}$:

如果重叠指数大于预定义的阈值 $0<m\leq 1$ 则将 $S_i$ 和 $S_j$ 合并。重复此操作直到没有子图能被合并为止,则此时剩下的子图即为我们所识别出来的类。

由于一个子图可能和多个子图重叠,从而有不同的合并顺序导致不同结果。为避免此情况,本方法对于有最大总大小 $\vert S_i\vert + \vert S_j \vert$ 的子图对给予高优先级。

上述过程所生成的类中依然会有重叠的部分,为了解决聚类过程中节点所属类别的唯一性问题,对于每个候选类别 $C$ 和目标节点 $v$ 定义一个接近度指标 $\text{Score}(C,v)$,由从 $C$ 中的节点入射到 $v$ 的边上的平均权重计算而来:

其中 $c_i$ 是 $C$ 中的一个点。然后将节点 $v$ 分入具有最大 $\text{Score}$ 值的类中,并将其从其他类中去除。这种做法会改变类的组成,且有可能会产生少于三个节点的类,这种较少数目的类被看做孤点,但在本文涉及的应用条件中未出现这种情况。

Corr

本方法核心思想在于新定义了一个体现成对细胞间相似性的差异相关系数(differentiability correlation)。

假设具有 $p$ 个基因表达量的 $n$ 个细胞数据,其第 $i$ 个细胞的基因表达情况表示为向量:

对于两个细胞 $i$ 和 $j$,定义两个基因集合 $V_{i_j}^+$ 和 $V_{i_j}^-$,其中 $V_{i_j}^+$ 表示细胞 $i$ 中表达水平高于除了细胞 $j$ 以外其他所有细胞的平均水平的基因, $V_{i_j}^-$ 则表示细胞 $i$ 中表达水平低于除了细胞 $j$ 以外其他所有细胞的平均水平的基因,即:

使用上述的 $V_{i_j}^+$ 和 $V_{i_j}^-$ 来给细胞 $i$ 中的基因 $k$ 定义差异状态,其规则为:

同理可以定义出细胞 $j$ 中基因 $k$ 的差异状态 $U_{j_{i} k}, k=1, \ldots, p$。于是可以定义细胞 $i$ 和细胞 $j$ 之间的不相似度:

其中

那么差异相关系数 $\left(1-S_{i j}\right)$ 则可以表示细胞 $i$ 和细胞 $j$ 之间相对于其他细胞的差异的相关性,并以此作为距离来进行层次聚类。

pcaReduce

用 $X_{n\times d}$ 表示基因表达矩阵,$n$ 为细胞数,$d$ 为基因数目,$x_i=\{x_{i1},\cdots,x_{id}\}$。再设 $Y_{n\times q}$ 为原数据向 $q$ 个主方向投影得到的矩阵,$Y_i$ 为一个细胞子集,$Y_i \subset Y$。

本方法首先对经过降维后的矩阵的前 $K-1$ 个主方向进行 $K$ 均值聚类,初始的聚类数目 $K$ 被设定为一个足够大的值,例如 30,来确保能包含大部分的细胞类型。

在初步聚类完成后,选取分别来自一对类 $(i,j)$ 的两个子集 $(Y_i,Y_j)$,计算这些观测值合并到一个类的概率,假设这个概率服从均值多元高斯分布,其协方差矩阵的计算方法为

其中 $(n_i,n_j)$,$(\mu_i,\mu_j)$ 和 $(\Sigma_i,\Sigma_j)$ 分别表示类 $i$ 和 $j$ 的大小、均值、协方差。对所有可能的对 $(i,j)$ 重复上述操作,然后选取两个类进行合并,方案有两种:

(1) 选取概率最高的一对类进行合并;

(2) 或者采样一对类以按其(归一化)合并概率成比例进行合并。

此时类数已降为 $K-1$,然后再对此结果进行主成分分析,保留前 $K-2$ 个主方向,去除第 $K-1$ 个主方向。以此方法重复进行降维和聚类,直到只剩一个类。如果使用了上述基于采样操作的合并方法,则整个过程可以重复进行来获得许多可选的聚类结果,这对于评估聚类结果的稳定性非常有用。

综上所述,pcaReduce 实际上是融合了主成分分析和 $k$ 均值聚类的具有层次性的聚类方法。

CIDR

PCA 等效于对从数据集导出的欧氏距离矩阵执行主坐标分析 (PCoA),假定只要能够可靠地估计存在缺失的情况下每对样本 (即单个单元格) 之间的差异,就无需明确估计缺失的值。

两个单细胞的基因表达样本 $C_{i}=\left(o_{1 i}, o_{2 i}, \ldots, o_{n i}\right)$ 和 $C_{j}=\left(o_{1 j}, o_{2 j}, \ldots, o_{n j}\right)$,其中 $o_{ki}$ 和 $0_{kj}$ 分别表示细胞 $C_i$ 和 $C_j$ 中基因 $k$ 的表达量,对于两个样本的欧式距离的平方有:

为简单起见,把基因表达数据中的所有 0 值作为 dropout 候选。通常来说,即使允许接近零的值作为 dropout 候选,本理论依然有效。对于上式中拆分成的三个部分,第一部分包括的样本都不为 0,对 dropout 不产生影响;第二部分都为 0 值,其结果也为 0 值(或接近 0 的很小值)。

对于 dropout 结果产生的主要影响在于第三部分,基因 $k$ 在一个细胞中表达值为 0 而在另一个细胞中不为 0。这时的 0 可以认为是真实情况下的无表达,也可以认为是非零表达基因的 dropout 事件。如果我们将所有的零值视作真实的无表达,也就是把零值认为是 dropout 事件的概率视为 0,这也就是直接将 PCA 直接应用于单细胞测序数据的情况,实际上是夸大了零值的真实性。

但是,已经观察到基因表达值丢失的可能性与真实表达水平成反比,也就是说,低表达水平的基因比高表达水平的更容易成为 dropout。基于此信息,假设可以通过给定的 dropou 概率分布在前式第三项中估算 dropout 候选者的表达值,来缩小 dropout 导致的夸大。这也就是 CIDR(Clustering through Imputation and Dimensionality Reduction) 算法的动机。

CIDR 算法分为五步:

(1) dropout 候选的识别

CIDR 使用 TPM(tags per million) 基因表达水平的对数分布作为判断是否为候选 dropout 的基础,这种分布在零值附近有一强峰值,在非零位置有一个或多个较小的峰值表示基因的表达水平。

CIDR 算法对于每个细胞 $C_i$ 都找出一个样本相关的阈值 $T_i$,表达水平低于此阈值的基因被视作候选的 dropout,这实际上将真实的 dropout 和真实的低水平表达样本都包含了进来。

(2) 估计 dropout 比率和基因表达水平之间的关系

设 $u$ 为细胞中某基因未观测到的真实表达水平,$P(u)$ 为该基因是 dropout 的概率,经验表明 $P(u)$ 为递减函数。CIDR 算法使用非线性最小二乘回归的方法将数据拟合为递减的对数函数(经验 dropou 比率与平均表达水平),将此函数作为 $P(u)$ 的估计,即 $\hat{P}(u)$。在此做出合理的假设,即数据集中的大多数 dropout 候选者实际上都是 dropout,这使得基因和细胞之间可以共享信息。

(3) 计算每对单细胞的预估基因表达谱之间的差异

上述 $\hat{P}(u)$ 用于在 CIDR 差异矩阵的计算中进行插补。考虑一对细胞 $C_i$ 和 $C_j$,相应地关于特征 $k$ 的表达水平为 $o_ki$ 和 $o_kj$,并设 $T_j$ 和 $T_j$ 为前述的 dropout 阈值。

插扑操作仅针对 dropout 候选者,即对 $o_i \geq T_i$ 和 $o_j \geq T_j$ 都不进行处理。对于二者中只有其一为 dropout 候选时,采取插补。例如当 $o_i \leq T_i$ 且 $o_j \geq T_j$ 时,对 $o_i$ 进行插补:

为了快速执行上述步骤,直接将式中的 $\hat{P}(u)$ 替换为一个更为简单的函数:

其中的 $T_W$ 默认为 0.5。此时,上述步骤的方法改进为:

具体的实验表明,这个改进确实能够提升速度,而且不会损害聚类精度。

然后,使用前述的 $\left[D\left(C_{i}, C_{j}\right)\right]^{2}$ 计算方法来计算 $C_i$ 和 $C_j$ 的差异矩阵。这种方法被称为隐式方法,因为对于每个细胞,其与不同细胞对比时,插补结果会有所不同。

(4) 使用上述步骤计算而来的 CIDR 差异矩阵进行 PCoA 方法降维

(5) 使用前几个主坐标进行聚类

CIDR 使用层次聚类方法进行聚类,并且使用 Calinski–Harabasz 指数确定聚类数。

参考文献

[1] Kiselev, V. Y., Andrews, T. S., & Hemberg, M. (2019). Challenges in unsupervised clustering of single-cell RNA-seq data. Nature Reviews Genetics, 20(5), 273–282. https://doi.org/10.1038/s41576-018-0088-9

[2] Xu, C., & Su, Z. (2015). Identification of cell types from single-cell transcriptomes using a novel clustering method. Bioinformatics, 31(12), 1974–1980. https://doi.org/10.1093/bioinformatics/btv088

[3] Jiang, H., Sohn, L. L., Huang, H., & Chen, L. (2018). Single cell clustering based on cell-pair differentiability correlation and variance analysis. Bioinformatics. https://doi.org/10.1093/bioinformatics/bty390/4996592

[4] žurauskiene, J., & Yau, C. (2016). pcaReduce: Hierarchical clustering of single cell transcriptional profiles. BMC Bioinformatics, 17(1), 1–11. https://doi.org/10.1186/s12859-016-0984-y

[5] Lin, P., Troup, M., & Ho, J. W. K. (2017). CIDR: Ultrafast and accurate clustering through imputation for single-cell RNA-seq data. Genome Biology, 18(1), 1–11. https://doi.org/10.1186/s13059-017-1188-0