深度学习插补单细胞中的DNA甲基化位点并增强精神分裂症表观遗传改变的检测.
0. TL; DR
scMeFormer 是一个基于 transformer 的深度学习模型，用于插补单细胞中的 DNAm 状态。模型包含一个 DNA 模块和一个 CpG 模块。DNA 模块利用 CNN 和 transformer 学习序列基序和长程依赖；CpG 模块则通过两个独立的 transformer 分别捕捉位点间和细胞簇间的甲基化关系。
通过在五个人类和小鼠的单细胞 DNAm 数据集上进行综合评估，scMeFormer 展现了优于其他模型的卓越性能，即使在 CpG 位点覆盖率降低到原始数据的10%时，仍能实现高保真度的插补。
1. 背景介绍
DNA 甲基化（DNAm）是一种基本的表观遗传机制，在基因调控、哺乳动物发育和多种人类疾病中扮演着关键角色。单细胞技术使得我们能够在单个细胞内分析 DNA 序列上胞嘧啶的甲基化状态，这为理解细胞异质性、发育过程和疾病状态提供了有力的补充。然而，由于单个细胞内的 DNA 物质有限以及技术本身的局限，当前技术往往面临 CpG 位点覆盖率稀疏的问题，通常只能测量到单个细胞中不到10%的 CpG 位点。这限制了在单细胞分辨率下全面揭示表观遗传图景的潜力。
为了应对单细胞 DNAm 数据集覆盖率稀疏的挑战，已经开发了一些计算方法来插补单个细胞中 CpG 位点的甲基化状态。
- 传统机器学习方法：如 Melissa（贝叶斯模型）、LightCpG（LightGBM）和 CaMelia（CatBoost），这些方法通常依赖于手工设计的特征或局部相似性模式。
- 基于深度学习的方法：如 DeepCpG（RNN）和 CpG Transformer，它们能够自动从 DNAm 谱和 DNA 序列中提取特征。然而，这些模型在扩展到数千个细胞的（这在当前研究中越来越普遍）大规模数据集时，面临着计算上的巨大挑战。最新的 GraphCpG 模型虽然性能先进，但它仅利用了 DNAm 谱，未考虑 DNA 序列特征，且一次只能预测一个细胞中的一个 CpG，难以扩展。
为此，作者提出了 scMeFormer，一个基于 transformer 的深度学习模型。它能够在一个多任务预测框架内，利用局部 DNA 序列和跨细胞的 DNAm 谱信息，高效地插补数千个单细胞中每个 CpG 位点的 DNAm 状态。
2. scMeFormer方法
scMeFormer 的模型架构包含两个主要模块：一个 DNA 模块和一个 CpG 模块，其输出最终被一个全连接网络整合以进行预测。
2.1 DNA 模块
DNA 模块学习 DNA 序列中的甲基化基序（motifs）。该模块由一个三层的卷积神经网络和一个八层的 transformer 编码器组成。
输入是一个长度为 $L$（实际模型中为2000 bp）的 DNA 序列，被编码成一个 $L \times 4$ 的独热矩阵。CNN 用于从 DNA 序列中提取局部特征。然后，这些局部特征被传递给 transformer，以学习局部 DNA 特征之间潜在的远程相互作用。
2.2 CpG 模块
CpG 模块利用预定义细胞簇内部和之间的邻近 CpG 位点的 DNAm 信息。该模块由一个三层 CNN 和两个独立的八层 transformer 组成。
首先，所有单细胞根据其 DNAm 水平被聚类成 $N$ 个细胞簇。对于一个待预测的目标 CpG，模型会取其周围的 $K$ 个（实际模型中为100）邻近 CpG 位点在每个细胞簇中的平均 DNAm 水平，形成一个 $K \times N$ 的矩阵作为输入。
CNN 通过一维卷积核识别每个簇内 CpG 位点的局部 DNAm 模式。CNN 输出的跨簇平均特征矩阵被送入位置感知的 Transformer (position-wise transformer)，以捕捉不同 CpG 位置之间的关系。位于目标 CpG 位置的特征矩阵（红色框）被传递给簇感知的 Transformer (cluster-wise transformer)，以捕捉不同细胞簇之间的 DNAm 关系。
2.3 全连接网络
全连接网络整合来自 DNA 和 CpG 模块的信息，并做出最终预测。
DNA 模块和 CpG 模块的输出特征向量被拼接起来，然后送入一个包含隐藏层、dropout 层和输出层的全连接网络。输出层的单元数等于细胞总数，每个单元通过 sigmoid 函数输出对应细胞中目标 CpG 位点的甲基化状态（0或1）。
2.4 模型训练
CpG 位点按染色体被划分为训练集（1-20号染色体）、验证集（21号染色体）和独立的测试集（22号染色体）。
为每个数据集单独训练一个模型。为了提高计算效率，作者还实现了一种策略，即在一个数据集上预训练一个模型，然后在其他新数据集上进行微调。
3. 实验分析
作者在五个人类和小鼠的单细胞 DNAm 数据集上，对 scMeFormer 的性能进行了全面的评估。
3.1 scMeFormer 的预测性能
在对22号染色体上的所有独立测试 CpG 位点进行的比较中，无论是从零开始训练（scMeFormer(scratch)）还是经过微调（scMeFormer(fine-tune)），scMeFormer 的性能（平均 AUPRC 分别为0.871和0.869）都显著优于仅使用 DNA 模块、仅使用 CpG 模块、纯 CNN 模型以及基于簇平均值的基线模型。
这表明 DNA 和 CpG 两个模块提供了互补的信息，并且 transformer 架构对于性能提升至关重要。值得注意的是，预训练-微调策略的性能与从零开始训练相当，但极大地降低了计算成本。该策略甚至可以成功地跨物种应用。
作者根据 CpG 位点在所有细胞中的变异水平，对测试集进行了分层评估。对于变异性最低的 CpG 位点（如变异性最低的10%），所有模型都表现良好（AUPRC > 0.97）。随着 CpG 位点变异性的增加，所有模型的预测性能都有所下降。
然而，在所有变异性水平上，scMeFormer 都始终优于其他模型。例如，在变异性最高的10%的 CpG 位点上，scMeFormer(fine-tune) 的平均 AUPRC 比 CNN 模型和簇模型分别高出0.083和0.039。
3.2 scMeFormer 在低 CpG 覆盖率下的性能
作者通过下采样 (downsampling) 模拟了 CpG 覆盖率较低的场景（从50%到1%）。结果显示，即使在仅有1%原始 CpG 覆盖率的极端情况下，scMeFormer 仍然保持了良好的性能（平均 AUPRC = 0.821），并始终优于基线的簇模型（平均 AUPRC = 0.728）。
3.3 插补质量指标提升插补性能
作者定义了一个插补质量得分 (imputation quality score)，用于衡量插补结果的可靠性。该分数计算的是目标 CpG 位点与其邻近 CpG 位点预测值之间的平均绝对差异。
结果显示，使用这个分数进行过滤，可以显著提高预测性能。分数阈值越严格（值越小），保留的插补数据质量越高，AUPRC 也越高。例如，在未下采样的 sn-m3C-seq 数据集上，使用0.1的阈值过滤后，AUPRC 从0.855提升到了0.935。
3.4 scMeFormer 在低覆盖率下恢复细胞簇
作者评估了在低覆盖率下，插补后的数据恢复原始细胞类型簇的能力。结果显示，即使在低至5%的下采样率下，使用 scMeFormer 插补（特别是经过滤后）的数据仍然能够很好地恢复原始的细胞簇结构（ARI > 0.6）。相比之下，未经插补的数据或使用簇模型插补的数据，其恢复细胞簇的能力则要差得多。
3.5 scMeFormer 增强细胞类型特异性 DMR 的检测
作者进一步评估了插补数据在识别差异甲基化区域 (DMRs) 方面的能力。在原始未插补的数据中，平均只能检测到几百个 DMRs。而在使用 scMeFormer 插补后，能够检测到数万到数十万个 DMRs。scMeFormer 插补后的数据能够以较高的召回率（recall rates）恢复在原始数据中发现的 DMRs。
通过分层LD得分回归 (S-LDSC) 分析，作者发现，从插补数据中识别出的细胞类型特异性 DMRs 显著富集了与该细胞类型相关的脑部疾病的遗传力，而与非脑部性状（如身高、糖尿病）无关。
3.6 scMeFormer 增强精神分裂症相关 DMR 的检测
最后，作者将 scMeFormer 应用于一个包含精神分裂症（SCZ）患者和对照组的单细胞 DNAm 数据集。
- 图 A: 在未插补的数据中，没有检测到任何与 SCZ 相关的 DMRs。而在插补后，发现了数千个 DMRs，其中大部分在 SCZ 患者中呈上调状态（即甲基化水平升高）。
- 图 B, C, D: 进一步的分析表明，这些 DMRs显著富集了 SCZ 的 GWAS 信号。其关联的基因也显著富集了 SCZ 的 GWAS 信号。与在 SCZ 患者中下调的差异表达基因（DEGs）显著富集。这些富集现象主要由上调的 DMRs 驱动，并且主要发生在兴奋性神经元中。
- 图 E: GO 富集分析显示，这些 DMRs 关联的基因富集于突触信号传导等通路。
这些发现共同指向兴奋性神经元中上调的 DMRs 在 SCZ 发病机制中的关键作用，并突显了 scMeFormer 在增强复杂疾病表观遗传学研究方面的巨大潜力。
