37、基于混合模型和神经网络的单细胞插补方法

基于混合模型和神经网络的单细胞插补方法

1. 引言

单细胞测序是在单细胞水平上对转录组或基因组进行扩增和测序，以检测来自多个组学（转录组学、表观基因组学和蛋白质组学）的单细胞数据。然而，单细胞分离并单独测序时，通量很低，主要受成本限制，随着待测单细胞数量增加，测序成本几乎呈线性增长。近年来采用条形码单细胞识别方法，为每个细胞添加独特的DNA序列，一次建库就能测量成百上千个单细胞的信息。

但每次建库后形成的细胞 - 基因矩阵会有很多零值，这些零值可能是真的未表达，也可能是表达值太小而未表达，这两种情况对结果的影响不同。单细胞插补技术应运而生，可有效解决上述问题。下游有很多分析方法可评估插补结果的质量，如聚类、归一化、差异表达等，其中聚类分析最常用。

当前的插补算法可分为基于模型的插补算法、基于平滑的插补算法和基于神经网络的插补算法。过去的算法更多考虑基因之间的相关性和细胞之间的相似性，但未考虑提取高表达基因对算法的作用，利用神经网络加快学习速度对插补结果也有益。因此，提出了一种新的插补方法，用统计模型筛选高表达基因，以高表达基因作为数据集训练模型，修正网络参数，训练后将低表达基因作为输入，输出网络预测的基因表达值，最后将结果插补到原始表达矩阵中。

2. 相关工作

• 混合模型 ：混合模型是几种不同模型的组合，在统计学中，它代表大群体中存在子群体的概率模型，通常拟合一个混合分布来表示大群体中观测值的概率分布。本实验中，将单细胞中的缺失部分表示为伽马分布，非缺失部分表示为正态分布。
• 正态分布 ：也称为高斯分布，在数学、物理和工程领域是非常重要的概率分布，对统计学的很多方面有重大影响。若随机变量X服从数学期望为μ、方差为∂²的正态分布，记为N(μ, ∂²)，其概率密度函数中，正态分布的期望值μ决定其位置，标准差∂决定分布的大小。当μ = 0，∂ = 1时为标准正态分布。正态分布在医学现象（如人群身高、血液中红细胞数量、血红蛋白量、胆固醇等）以及实验中的随机误差中呈现为正态或近似正态分布。在单细胞测序中，用正态分布曲线标记正常未缺失的细胞。
• 概率密度函数：(f (x) = \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{(x - \mu)^2}{2\partial^2} \right))
• 伽马分布 ：是统计学的连续概率函数，在概率统计中是非常重要的分布，解决的问题是“n个随机事件发生需要多长时间”，可看作n个指数分布的独立随机变量之和。参数α称为形状参数，β称为逆尺度参数，当α = 1时，变为指数分布。
• 概率密度函数：(f_x(x) = \begin{cases} \frac{\lambda^{\alpha}x^{\alpha - 1}e^{-\lambda x}}{\Gamma(\alpha)}, & x > 0 \ 0, & x \leq 0 \end{cases})
• EM算法 ：期望最大算法是一种通过迭代进行最大似然估计的优化算法，通常用作缺失数据的概率模型进行参数估计，广泛用于处理数据中的缺失值以及许多机器学习算法（如高斯混合模型和隐马尔可夫模型）的参数估计。标准计算框架由交替的E步和M步组成：
  1. E步 ：使用隐藏变量的现有估计值，计算它们的最大似然估计。
  • (Q(z) = p(z|x; \theta))
    2. M步 ：最大化E步中获得的最大似然值，计算参数的值。
  • (\theta = \arg\max_{\theta} \sum_{z} Q_i(z) \log \frac{p(z, x; \theta)}{Q(z)})

3. 方法

该方法基于深度学习，使用TensorFlow实现，有以下五个步骤：
- 预处理 ：
1. 过滤 ：原始输入矩阵是细胞 - 基因矩阵，遍历整个矩阵，若一个细胞样本中表达数值的基因少于200个，过滤该细胞样本（从二维矩阵中删除该列）；若一个基因有表达值的细胞数量少于3个，过滤该基因（从二维矩阵中删除该行）。删除不符合条件的基因和细胞后，得到计数矩阵。
2. 筛选线粒体基因 ：细胞分离过程中，单细胞测序可能导致细胞损伤或建库失败，引入低质量数据，表现为细胞表达、基因表达降低或线粒体基因比例增加，需在分析前去除这些低质量的行和列。
3. 归一化 ：使用TPM处理，先对基因长度进行归一化，再对测序深度进行归一化，公式为(TMP_i = 10^6 * \frac{n_i/l_i}{\sum_{j} n_j/l_j})，该方法更适合组内重复样本的比较和组内样本的比较。然后使用伽马 - 正态模型识别高表达基因。
- 识别高表达基因 ：
1. 聚类 ：将细胞定义为节点，细胞中的基因表达值作为特征值进行k - 均值聚类，用欧氏距离衡量细胞中基因表达的相似性，将特征值相似的细胞聚集为一个细胞子群，默认k为5。
- (c(i) := \arg\min_{j} \left\lVert x^{(i)} - \mu_j \right\rVert^2)
2. 区分高低表达基因 ：使用伽马 - 正态混合模型，伽马模型统计缺失值较多的细胞子集，正态模型统计缺失值较少的细胞子集，识别出的高表达基因将被保留。由于矩阵中的零值不一定都来自缺失现象，需构建统计模型判断这些零值是否真的缺失，提前统计缺失率以确定高表达基因。用两个统计模型确定表达内容，非零值可用伽马 - 正态混合模型统计，缺失值大的部分用伽马模型统计（即低表达基因），缺失值小的部分用正态模型确定（即高表达基因），每个聚类都有单独的混合模型，根据统计模型给出表达值内容，默认截断阈值为0.5。
- 伽马分布：(f (X, \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(a)}X^{\alpha - 1}e^{-\beta x})
- 正态分布：(f (X, \mu, \sigma^2) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{(X - \mu)^2}{2\sigma^2}})
- 假设一个聚类中低基因分布的概率为ρ，基因(x_i)的密度函数满足：(f (x_i, \alpha, \beta, \mu, \sigma^2) = (1 - \rho) \frac{\beta^{\alpha}}{\Gamma(a)}x_i^{\alpha - 1}e^{-\beta x} + \rho \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{(x_i - \mu)^2}{2\sigma^2}})
- 使用EM算法推导伽马 - 正态模型的最大期望，迭代更新直到结束，选择阈值0.5，若最终结果大于0.5，为高表达基因，否则为低表达基因。
- 网络架构 ：
每个聚类使用基于之前聚类结果的新神经网络，将每个聚类中选择的高表达基因分块输入神经网络训练。网络架构有五层：
1. 输入层 ：包含512个神经节点。
2. 全连接层 ：激活函数为ReLu，包含256个细胞节点。
3. Dropout层 ：选择20%的丢弃率，目的是减少参数之间的相互依赖，缓解过拟合问题。
4. 全连接层 ：包含512个细胞节点。
5. 全连接层 ：包含512个细胞节点，添加这两层是为了更好地增加细胞之间的连接，将学习到的特征映射到全局样本空间。
损失函数使用均方误差（MSE）来减少输入层和输出层之间的误差，表达值较高的基因权重更高。高表达基因中，95%的数据用作训练集，5%用作测试集，网络初始学习率设为0.0001，每个epoch周期会自动修改学习率，优化器为Adam。将筛选出的高表达基因输入神经网络的训练参数中，输出结果是低表达基因应反映的表达值，最后将结果替换原始矩阵中的值。

数据处理步骤	操作内容
过滤	遍历细胞 - 基因矩阵，删除表达数值基因少于200个的细胞样本列和有表达值细胞少于3个的基因行
筛选线粒体基因	去除细胞分离和建库过程中引入的低质量数据（表现为细胞表达、基因表达降低或线粒体基因比例增加）
归一化	使用TPM处理，公式为(TMP_i = 10^6 * \frac{n_i/l_i}{\sum_{j} n_j/l_j})
聚类	以细胞为节点，基因表达值为特征值进行k - 均值聚类，默认k为5
区分高低表达基因	使用伽马 - 正态混合模型，通过EM算法推导最大期望，阈值0.5区分高低表达基因
网络训练	高表达基因分块输入五层神经网络训练，损失函数用MSE，95%训练集，5%测试集，初始学习率0.0001，优化器Adam

graph LR
    A[预处理] --> B[过滤]
    A --> C[筛选线粒体基因]
    A --> D[归一化]
    D --> E[识别高表达基因]
    E --> F[聚类]
    E --> G[区分高低表达基因]
    G --> H[网络架构]
    H --> I[训练高表达基因]
    I --> J[预测低表达基因]
    J --> K[插补结果到原始矩阵]

基于混合模型和神经网络的单细胞插补方法

4. 插补

在单细胞测序中，插补的作用是补充因建库过程中意外片段化而未表达的基因表达值，便于后续对单细胞测序数据进行分析，提高后续分析结果的准确性。具体操作步骤如下：
1. 划分数据集 ：将筛选出的高表达基因分为训练样本数据集和测试样本数据集。
2. 训练模型 ：使用训练样本数据集对神经网络模型进行训练，得到预测模型。
3. 预测结果 ：将聚类中的低表达基因输入到训练好的预测模型中，得到预测结果。
4. 替换值 ：用预测结果替换原始矩阵中低表达基因的表达值，并且是整行基因进行替换。替换后，原始矩阵中大量的零值显著减少，方便进行下游分析。

5. 评估

为了评估插补效果，采用了以下两种指标：
- PCC（皮尔逊相关系数） ：广泛用于衡量两个变量之间的相关程度，其值介于 -1 和 1 之间。两个变量的皮尔逊相关系数定义为两个变量之间的协方差与标准差的商，可用于衡量插补基因与原始基因之间的相关性。公式如下：
- (Corr = \frac{Cov(X (gene), Xmethod(gene))}{\sqrt{Var(X (gene)) * Var(Xmethod(gene))}})
其中，(X) 是输入矩阵，(Cov) 是协方差矩阵，(Var) 是方差。
- MSE（均方误差） ：是一种更方便的衡量“平均误差”的方法，可评估数据的变化程度。MSE 的值越小，预测模型描述实验数据的准确性越高。公式如下：
- (MSE = \sum (X (gene) - Xmethod(gene))^2)

6. 实验

6.1 数据集

选择了四个数据集，包括两个小鼠数据集和两个人类数据集，这些都是公开数据集。小鼠数据集可从 GEO 数据库下载，分别为 GSE67602（小鼠表皮单细胞）和 GSE60361（小鼠大脑皮层单细胞）；人类数据集为 Jurkat 和 293T，可从 10X Genomes 官方网站下载。

6.2 预处理结果

对数据集进行预处理过滤后，原始矩阵中许多不相关的基因和死细胞被筛选掉。以下是四个数据集在筛选前后细胞和基因数量的对比：
| 数据集 | 筛选前 | 筛选后 |
| ---- | ---- | ---- |
| Jurkat | 3258 × 26024 | 3258 × 15269 |
| 293T | 2885 × 32738 | 2885 × 16316 |
| GSE60361 | 3005 × 16382 | 3005 × 15602 |
| GSE67602 | 1442 × 16383 | 1442 × 13001 |

从表中可以看出，细胞数量基本保持不变，而基因数量显著减少。经过筛选后的细胞基因矩阵被称为计数矩阵。

6.3 插补结果

本实验与 DeepImpute 算法进行了比较，对比了皮尔逊相关系数和均方误差这两个评估指标。结果表明，使用该方法后，皮尔逊相关系数显著提高，插补结果的相关性明显增强；均方误差也相应降低。具体数据对比可参考以下表格：
| 评估指标 | 本方法 | DeepImpute 算法 |
| ---- | ---- | ---- |
| 皮尔逊相关系数 | 显著提高 | 较低 |
| 均方误差 | 降低 | 较高 |

综上所述，基于混合模型和神经网络的单细胞插补方法在处理单细胞测序数据时，能够有效识别细胞中的缺失值，增强细胞群体之间的聚类效果，提高差异表达分析的能力。通过对多个数据集的实验验证，该方法在插补效果上优于传统的 DeepImpute 算法。

graph LR
    A[实验数据集] --> B[预处理]
    B --> C[插补]
    C --> D[评估]
    D --> E[与 DeepImpute 算法对比]