DeepCpG：使用深度学习准确预测单细胞DNA甲基化位点.

paper：DeepCpG: accurate prediction of single-cell DNA methylation states using deep learning

0. TL; DR

DeepCpG 是一个基于深度神经网络的计算方法，用于预测单细胞中的甲基化状态。DeepCpG包含三个模块：DNA 模块、CpG 模块和联合模块，分别用于学习 DNA 序列模式、邻近 CpG 位点的共甲基化模式，并将两者结合进行最终预测。

作者在五种不同细胞类型的单细胞甲基化数据上对 DeepCpG 进行了评估，这些数据由不同的测序方案生成。结果显示，DeepCpG 的预测比以前的方法要准确得多。此外，作者还证明了模型的参数是可解释的，从而为理解序列组成如何影响甲基化变异提供了新的见解。

1. 背景介绍

DNA 甲基化是研究最广泛的表观遗传标记之一，它与染色体不稳定性、X染色体失活、细胞分化、癌症进展和基因调控等多种生物学过程有关。

尽管已经有成熟的方案用于量化细胞群体的平均 DNA 甲基化水平，但最近的技术进步已经实现了在单细胞分辨率下对 DNA 甲基化进行分析，例如全基因组亚硫酸氢盐测序 (scBS-seq) 或简化代表性方案 (scRRBS-seq)。这些技术为理解单细胞中 DNA 甲基化的调控和动态变化提供了前所未有的机会。

然而，由于每个细胞的起始 DNA 量极少，单细胞甲基化分析固有地受限于中等的 CpG 覆盖度（scBS-seq 约为20-40%，scRRBS-seq 约为1-10%）。因此，预测缺失的甲基化状态以实现全基因组分析，成为一个至关重要的第一步。

尽管已有一些方法用于预测细胞群体的平均 DNA 甲基化谱，但这些方法并未考虑细胞间的变异性，并且通常需要预先定义的特征和基因组注释，这些注释往往只适用于有限的细胞类型和条件。

为此，作者提出了 DeepCpG，一个基于深度神经网络的计算方法，用于预测单细胞甲基化状态并建模 DNA 甲基化变异的来源。与以往的方法不同，DeepCpG 不将特征提取和模型训练分开，而是通过一个模块化的架构，以数据驱动的方式学习具有预测性的 DNA 序列和甲基化模式。

2. DeepCpG 方法

DeepCpG 的核心是一个模块化的深度学习架构，它通过训练来从局部的 DNA 序列窗口和观测到的邻近甲基化状态中预测二元的 CpG 甲基化状态。

2.1 DeepCpG 模型架构

DeepCpG 由三个核心模块组成：

2.1.1 DNA 模块

该模块旨在从 DNA 序列中检测信息丰富的模式。模块架构是一个包含多个卷积层、池化层和一个全连接隐藏层的卷积神经网络。

卷积层：使用多个卷积核扫描输入序列，以识别不同的序列基序。
池化层：对卷积层的输出进行最大池化，以总结邻近神经元的激活值并降低维度。
全连接层：最后，一个全连接层用于学习不同基序之间的相互作用。

输入以待预测的 CpG 位点为中心，长度为1001 bp的 DNA 序列。该序列通过独热编码 (one-hot encoding) 转换成一个二进制矩阵。

2.1.2 CpG 模块

该模块旨在捕捉邻近 CpG 位点之间在单个细胞内和跨细胞间的相关性。模块架构是一个双向门控循环网络 (bidirectional gated recurrent network, GRU)。GRU 是一种特殊的循环神经网络 (RNN)，它能够将可变长度的输入序列压缩成一个固定大小的特征向量，从而有效捕捉细胞间的依赖关系，而不受细胞数量的影响。

嵌入层：对输入的甲基化状态和距离进行变换，以学习它们之间的相互作用。
双向 GRU：从两个方向扫描细胞序列（即按细胞顺序排列的邻近 CpG 信息）。

对于一个目标 CpG 位点，其输入是该位点左右各25个邻近 CpG 位点的甲基化状态和距离，这些信息来自所有可用的细胞。

2.1.3 联合模块 (Joint Module)

该模块负责整合来自 DNA 模块和 CpG 模块的信息，并作出最终的预测。模块架构是两个全连接隐藏层，最后是一个包含 $T$ 个（$T$ 为细胞数量）sigmoid 神经元的输出层。

将 DNA 模块和 CpG 模块的输出特征向量进行拼接。通过两个全连接层学习这两种特征之间的相互作用。最后，输出层为每个细胞中的目标 CpG 位点预测一个0到1之间的甲基化概率。

\[\hat{y}_{nt}(x) = \text{sigmoid}(x) = \frac{1}{1+e^{-x}}\]

2.2 模型训练

通过最小化加权负对数似然 (negative log-likelihood) 损失函数来学习模型参数。该函数衡量了预测的甲基化率与观测到的二元甲基化状态之间的差异。

\[\mathcal{L}(w) = \text{NLL}_w(\hat{y}, y) + \lambda_2 \|w\|^2\]

其中，$\text{NLL}_w(\hat{y}, y) = -\sum_{n=1}^N \sum_{t=1}^T o_{nt}[y_{nt}\log(\hat{y}_{nt}) + (1-y_{nt})\log(1-\hat{y}_{nt})]$，$o_{nt}$ 是一个指示符，当位点 $(n,t)$ 有观测值时为1。

使用带动量和自适应学习率的小批量随机梯度下降和 Adam 优化器进行训练。同时，使用 L2 权重衰减和 Dropout 来防止过拟合。

采用预训练-微调的策略。首先，独立地预训练 DNA 模块和 CpG 模块。然后，在训练联合模块时，固定预训练好的 DNA 和 CpG 模块的参数，只优化联合模块的隐藏层和输出层。

3. 实验分析

作者在多种单细胞甲基化数据集上对 DeepCpG 进行了全面的性能评估。

3.1 单细胞甲基化状态的准确预测

作者将 DeepCpG 与多种基准方法进行了比较，包括：

WinAvg 和 CpGAvg：简单的局部或跨细胞均值插补。
RF：使用与 DeepCpG 类似特征的随机森林分类器。
RF Zhang：一个使用了全面的基因组和表观组学注释的先进随机森林模型。

结果分析：

图 a: 在 scBS-seq 数据集上，DeepCpG 在全基因组范围内的预测性能（以 AUC 衡量）显著优于所有其他方法，包括依赖于大量外部注释的 RF Zhang。这表明 DeepCpG 能够从原始 DNA 序列中自动学习到更高级别、更具预测性的特征。
图 b: 作者对 DeepCpG 和 RF 模型的不同模块进行了分解评估。结果显示邻近的 CpG 甲基化状态比 DNA 序列特征更具预测性。结合两种特征的模型性能最佳。DeepCpG 的 CpG 模块性能优于结合了两种特征的 RF 模型，这得益于其循环网络架构能有效传递跨细胞信息。在只使用 DNA 序列信息时，DeepCpG DNA 模块相对于 RF DNA 模块的性能提升最大（AUC 0.83 vs 0.80），证明了 CNN 在从长序列窗口中提取预测性特征方面的强大能力。
图 c: 在不同 CpG 覆盖度的基因组区域中，DeepCpG 的性能优势在低覆盖度区域最为明显。
图 d: 在不同的基因组上下文中，所有模型在 GC 丰富区域表现最好。而 DeepCpG 的优势在 GC 贫乏区域（如启动子、增强子）中最为突出，这些区域通常与更高的细胞间甲基化变异性相关。
图 e: DeepCpG 在其他多种细胞类型（包括低甲基化的人类癌细胞）和不同的测序方案（scRRBS-seq）的数据上，同样表现出一致的优越性能，证明了其广泛的适用性。

3.2 DNA基序和单核苷酸突变对甲基化状态影响的评估

通过对 DNA 模块第一个卷积层学习到的卷积核进行可视化，作者发现了128个与甲基化相关的序列基序。这些基序根据其对甲基化水平的影响（促进或抑制）以及核苷酸组成，可以被清晰地分为两大类。

与甲基化水平降低相关的基序富含 CG，主要出现在启动子等区域；而与甲基化水平增加相关的基序富含 AT，主要出现在 CG 贫乏区域。其中20个基序与 CIS-BP 和 UniPROPE 数据库中的已知基序显著匹配，包括 CTCF, E2f 等与 DNA 甲基化和细胞分化密切相关的转录因子。

作者还使用一种基于梯度的方法，高效地估计了单核苷酸突变对 CpG 甲基化的影响。结果显示，距离目标 CpG 位点越近的突变，其影响越大。在 CG 密集的区域（如 CpG 岛），突变的影响相对较小，表明这些区域的甲基化状态对单核苷酸突变更具鲁棒性。作者还发现，已知的甲基化数量性状基因座 (mQTLs) 变异比随机变异具有显著更大的影响效应，这证明了模型预测的真实性。

3.3 发现与甲基化变异性相关的DNA基序

作者进一步训练了一个新的神经网络，用于联合预测每个 CpG 位点的平均甲基化水平和细胞间的变异性（方差）。

图 d：该模型能够以相当的准确性预测细胞间的变异性（Pearson’s R = 0.44）。
图 a, b, c：通过分析不同基序对平均甲基化水平和变异性的不同影响，作者识别出了22个主要与细胞间变异性相关的基序。这些基序大多富含 AT，并倾向于出现在进化上不保守的区域，如活跃的增强子。相比之下，主要影响平均甲基化水平的基序则富含 CG，并富集在进化上保守的区域，如基因启动子。

这些结果表明，DeepCpG 不仅是一个高精度的预测工具，还是一个强大的探索性工具，能够为理解 DNA 甲基化的调控机制提供深刻的生物学洞见。