GraphCpG: 基于位点感知近邻子图插补单细胞甲基化组.

0. TL; DR

GraphCpG 方法基于位点感知的邻近子图 (locus-aware neighboring subgraphs) 来插补甲基化矩阵。GraphCpG 仅使用 CpG 甲基化矩阵，通过提取每个待插补位点的邻近子图 (neighboring subgraph)，并利用关系图卷积网络 (R-GCN) 进行建模，有效捕捉了细胞间和位点间的复杂依赖关系。

在包含数百个细胞的大型、高稀疏性数据集上，GraphCpG 的性能超越了包括 DeepCpG, CaMelia 和 CpG Transformer 在内的所有先进方法。随着细胞数量的增加，其性能优势愈发明显，并且计算时间显著减少。

1. 背景介绍

DNA 甲基化是一种重要的表观遗传标记，与衰老、肿瘤发生等生物学过程密切相关。单细胞 DNA 甲基化测序技术，如 scRRBS-Seq, scBS-Seq 和 snmC-seq，使得我们能够在单个细胞层面探索 DNA 甲基化的异质性和动态变化。

然而，这些技术的一个主要瓶颈是测量覆盖率低，导致数据中存在大量的缺失值（稀疏性高达60%-99%）。这种高稀疏性严重阻碍了对潜在生物学过程的下游分析和理解。因此，开发有效的插补 (imputation) 技术来填补这些缺失值至关重要。

近年来，已经发展出多种用于单细胞甲基化状态插补的计算方法，它们大致可分为两类：

  1. 传统机器学习方法：
    • LightCpG：使用 LightGBM 模型，结合了 CpG 位置、DNA 序列和 CpG 岛等手工设计的特征。
    • Melissa：一个贝叶斯模型，通过聚类细胞来共享信息。
    • CaMelia：使用 CatBoost 模型，基于细胞间的局部甲基化模式相似性进行预测。
  2. 深度学习方法：
    • DeepCpG：基于 CNN 和 RNN，结合了 DNA 序列模式和甲基化状态。
    • CpG Transformer：受 Transformer 启发，使用二维滑动窗口自注意力机制处理甲基化矩阵。

尽管这些方法取得了一定的成功，但它们仍然面临挑战。传统机器学习方法依赖于复杂的手工特征工程；而基于深度学习的方法，如 CpG Transformer，虽然性能优越，但其模型复杂度随着细胞数量的增加呈二次方增长，限制了其在未来更大规模数据集上的可扩展性。

随着单细胞测序技术的发展，数据集的规模（特别是细胞数量）越来越大，数据也越来越稀疏。为了应对这一挑战，迫切需要一个既准确又高效、可扩展的插补框架。为此，作者提出了 GraphCpG，一个基于图深度学习的方法。

2. GraphCpG 方法

GraphCpG 的核心思想是将甲基化矩阵看作一个二分图，并将缺失值插补问题转化为图上的链接预测问题。它通过提取位点感知的邻近子图 (locus-aware neighboring subgraphs)，并利用图神经网络进行建模。

2.1 模型输入：二分图的构建

GraphCpG 的输入是一个由甲基化矩阵 $M$ 构建的无向二分图 $G$。图中有两类节点，一类是细胞 $u$（对应矩阵的行），另一类是 CpG 位点或称基因座 $v$（对应矩阵的列）。

细胞节点 $u$ 和基因座节点 $v$ 之间是否存在边，取决于 $M_{u,v}$ 的甲基化状态。作者将这个问题简化，定义了两种类型的边：甲基化 (methylated) 状态（$m=1$）和未甲基化 (unmethylated) 状态（$m=0$）。

2.2 邻近子图提取 (Neighboring Subgraph Extraction)

对于每一个待插补的目标 CpG 位点 $(u_t, v_t)$，作者不是对整个大图进行操作，而是提取其周围的一个邻近子图。以目标基因座 $v_t$ 为中心，通过一个宽度为 $W$ 的滑动窗口，水平地选取其邻近的基因座。这个子图包含了所有 $N$ 个细胞和这 $W$ 个邻近基因座，以及它们之间所有已知的甲基化状态（边）。

2.3 位点感知编码 (Locus-aware Encoding)

为了让图神经网络能够区分不同节点在子图中的角色，作者提出了一种位点感知编码方案。使用不同的整数标签来编码不同的节点角色，然后进行独热编码。

通过这种方式，模型能够学习到不同基因座位置的特定交互模式，以及细胞与基因座之间的交互模式。

2.4 图神经网络的模块化架构

GraphCpG 的训练部分是一个序列化的模型，依次包含 R-GCN 模块、CNN 模块和 MLP 模块。

2.4.1 关系图卷积网络 (R-GCN) 模块

R-GCN 是一种能够处理多关系图的 GNN。在这里，它被用来学习由甲基化（$m=1$）和未甲基化（$m=0$）两种边类型引入的丰富子图模式。

在每一层 $l$，节点 $i$ 的特征向量 $x^{l+1}_i$ 是通过聚合其邻居节点 $j$ 的信息来更新的。

\[x^{l+1}_i = \sigma \left( \sum_{m \in \mathcal{M}} \sum_{j \in \mathcal{N}_m(i)} \frac{1}{|\mathcal{N}_m(i)|} W^l_m x^l_j + W^l_0 x^l_i \right)\]

其中，$W^l_m$ 和 $W^l_0$ 是可学习的参数矩阵，$\mathcal{N}_m(i)$ 是节点 $i$ 在关系 $m$ 下的邻居集合。

通过堆叠 $L$ 层 R-GCN，并将每一层的输出拼接起来，得到每个节点的最终表示 $h_i$。

2.4.2 CNN 模块和 MLP 模块

将从 R-GCN 得到的、代表整个子图的节点表示矩阵 $H_{u_t, v_t}$，通过一个二维 CNN 和最大池化层，压缩成一个图级别的向量 $g$。

最后，将图级别向量 $g$ 送入一个多层感知机（MLP），并通过 sigmoid 激活函数输出最终的预测甲基化状态 $\hat{m}_{u_t,v_t}$。

\[\hat{m}_{u_t,v_t} = \text{sigmoid}(w^T \text{ReLU}(g W^T))\]

3. 实验分析

作者在五个公开的单细胞甲基化数据集上，将 GraphCpG 与 DeepCpG, CpG Transformer 和 CaMelia 进行了全面的性能比较。

3.1 性能比较

在五个数据集上的比较结果显示：在细胞数量较少的数据集（HCC, MBL）上，GraphCpG 取得了与最先进方法相当的性能。在细胞数量较多、数据更稀疏的数据集（Hemato, Neuron-Mouse, Neuron-Homo）上，GraphCpG 的性能（AUROC 和 MCC）全面超越了所有其他方法。

3.2 GraphCpG 的研究

作者通过消融实验研究了位点感知编码的重要性。与完整模型相比，移除对具体位点位置的感知（即所有基因座节点编码相同）后，模型性能在所有数据集上都有明显下降。如果完全不使用任何编码（即所有细胞和基因座节点编码都相同），性能会进一步轻微下降。这证明了通过位点感知编码来区分不同基因座的相对位置对于模型的性能至关重要。

作者可视化了 HCC 数据集中一个甲基化和一个未甲基化位点的邻近子图。结果直观地显示，甲基化和未甲基化的子图具有截然不同的模式。甲基化的子图通常具有更高的细胞和基因座平均甲基化水平。

通过计算节点嵌入间的余弦相似度，可以发现具有相似链接模式的节点（无论是细胞还是基因座）在嵌入空间中也更相似。这表明 GraphCpG 能够学习到有意义的节点表示，并利用这些表示来过滤相似信息，从而实现更准确的预测。

3.3 计算加速

作者比较了不同模型随着细胞数量增加的训练时间。结果显示，GraphCpG 的时间消耗增长最慢，是最有效率的模型。相比之下，CpG Transformer 的时间消耗随着细胞数量呈二次方增长，而 CaMelia 在细胞数量超过100后，其预处理时间急剧增加。

作者还发现，在大数据集上，GraphCpG 无需遍历整个训练集，仅使用1%甚至0.01%的训练位点，就能在第一个周期内达到接近最优的性能，这极大地节省了训练时间。

3.4 下游分析

为了展示 GraphCpG 在下游分析中的应用价值，作者在 Hemato 数据集上进行了细胞聚类分析。结果显示，使用经过 GraphCpG 插补后的数据进行层次聚类 (hierarchical clustering)，能够比使用原始稀疏数据更准确地识别出不同的细胞类型（如 CLP, CMP, GMP 等）。此外，插补后的数据在差异甲基化分析中也表现出与批量测序（bulk）数据更高的一致性。

这些结果表明，GraphCpG 的高质量插补能够显著增强下游生物学分析的准确性和可靠性。