0. TL; DR

MethylQUEEN (Methylation Language Model based on Quintuple Bidirectional Transformer) 是一种能够感知甲基化状态并覆盖全基因组甲基化图谱的预训练 DNA 甲基化基础模型。将来自 BS-seq 的甲基化信息直接编码为五进制碱基（quintuple base）序列进行预训练，而不是简化为 beta 值，从而保留了read级别的丰富信息。

MethylQUEEN在多种甲基化相关的下游任务中表现卓越，如组织溯源、基因表达预测和疾病检测。通过分析注意力分数，能够识别出对模型决策至关重要的 CpG 位点，并揭示了潜在的新型转录因子基序。在液体活检场景中，成功区分了早期肝癌患者与健康人的血浆样本，展示了其在早期癌症诊断中的巨大潜力。

1. 背景介绍

DNA 5-甲基胞嘧啶 (5mC) 是最常见的表观遗传修饰，主要发生在哺乳动物基因组的 CpG 二核苷酸上。它通过影响染色质结构和转录活性，在调控基因表达、维持基因组稳定性和指导细胞分化等方面发挥着关键作用。异常的甲基化模式，如肿瘤抑制基因的高甲基化或全基因组的低甲基化，与癌症等多种疾病密切相关。

高通量测序技术，特别是亚硫酸氢盐测序 (BS-seq)，已经使我们能够在单碱基分辨率上表征 5mC，极大地丰富了我们对其动态作用的理解。然而，如何从海量的 BS-seq 数据中解锁隐藏的洞见，并高效地提取表观遗传信息，仍然是一个重大的计算挑战。

近年来，大规模预训练模型在自然语言处理领域取得了巨大成功。类似地，作为遗传信息的主要载体，核酸序列也完全可以被视为一种生物语言。尽管已有一些优秀的核酸序列模型，如用于 DNA 的 Evo 和 DNABERT，以及用于 RNA 的 Uni-RNA，但它们在处理表观遗传信息方面存在局限。现有研究要么仅在微调阶段引入表观遗传信息，要么将其简化为 beta 值进行预训练，这两种方式都忽略了单个生物样本中蕴含的丰富信息。因此，迫切需要一个能够直接在预训练阶段整合来自 BS-seq 数据的原始表观遗传信息的基础模型。

为此，作者开发了 MethylQUEEN，这是一个基于高质量全基因组 5mC 甲基化数据、能够感知甲基化状态的预训练模型。它采用了一种可扩展且易于推广的修饰编码方案，旨在捕捉序列特征及其背后隐藏的表观遗传信息。作为首个直接在 BS-seq 数据上训练的基础模型，MethylQUEEN 标志着甲基化分析的范式转变。

2. MethylQUEEN 方法

作者构建 MethylQUEEN 的核心在于其创新的数据编码方式和高效的预训练策略，旨在将原始的甲基化状态信息无损地融入模型学习中。

2.1 数据收集与预处理

作者从两大著名的 DNA 测序项目（ENCODE 和 NIH Roadmap Epigenomics）收集了超过500个生物样本的全基因组亚硫酸氢盐测序 (WGBS) 数据。这些数据经过了标准化的 BS-seq 比对流程处理，包括：

质量控制：使用 trim_galore 进行质量剪枝。
序列比对：使用 bwameth 将读段比对到 GRCh38 人类参考基因组。
去重与过滤：移除重复读段和低质量比对，并过滤掉亚硫酸氢盐转化不完全的读段。

经过处理后，作者开发了一种五进制碱基 (quintuple base) 编码方案：

未甲基化的胞嘧啶 (C)：在测序中被转化为胸腺嘧啶 (T)，因此在序列中仍表示为 C（基于其在参考基因组中的原始状态）。
甲基化的胞嘧啶 (mC)：在测序中保持不变，因此被编码为一个新的字符 M。

这样，原本由 “A”, “T”, “C”, “G” 组成的 DNA 序列就被扩展为由 “A”, “T”, “C”, “G”, “M” 组成的五进制序列。接着，作者将这些较短的测序读段（reads）通过重叠比对，组装成更长的序列（平均长度 > 600 bp），最终获得了超过5亿条高质量的、包含甲基化信息的序列。

2.2 模型架构与预训练

MethylQUEEN 的架构基于 MosaicBERT，并针对甲基化数据的特性进行了优化，以提高计算效率和内存使用率。

2.2.1 令牌化 (Tokenization)

作者使用字节对编码 (Byte Pair Encoding, BPE) 算法对五进制序列进行令牌化。训练 BPE tokenizer 时，作者特别关注了那些包含 “CG” 或 “M” 的 token，以及单个碱基 “C”，并将它们定义为甲基化相关令牌 (methyl-tokens)。

2.2.2 核心架构

MethylQUEEN 的核心架构包含以下优化特性：

ALiBi 位置编码：使用带有线性偏置的注意力 (Attention with Linear Biases, ALiBi) 代替传统的位置嵌入，使其能更好地处理不同长度的序列。
FlashAttention：采用 FlashAttention 和低精度层归一化 (Low Precision Layer Normalization) 来优化计算效率并减少内存占用。
GeGLU 激活函数：使用门控线性单元 (Gated Linear Unit) 配合 GeLU 激活函数，以提升 Transformer 模型的性能。

2.2.3 甲基化令牌倾向的预训练 (Methyl-token-prone Pre-training)

作者设计了一个两阶段的预训练策略：

第一阶段：进行标准的掩码语言模型 (Masked Language Modeling, MLM) 预训练。以15%的概率随机掩盖 token，让模型学习恢复被掩盖的内容。
第二阶段：采用甲基化令牌倾向的掩码策略。具体来说，非甲基化相关令牌的掩码概率设为0%，而甲基化相关令牌的掩码概率被大幅提高（约97%）。这一策略迫使模型更加关注和学习与甲基化相关的信息。

整个预训练过程在8个 NVIDIA A800 80GB GPU 上完成。

3. 实验分析

作者通过一系列精心设计的实验，从不同维度验证了 MethylQUEEN 的性能和应用价值。

3.1 预训练效果与零样本分析 (Zero-shot Analysis)

作者首先探究了预训练是否使模型学到了有意义的表观遗传学知识。

图 a, b: UMAP 可视化显示，在预训练后，所有的甲基化相关令牌在嵌入空间中被清晰地聚为一类，与其他 token 分离开来，表明模型成功识别并区分了与甲基化相关的序列模式。
图 c, d: 作者进一步分析了来自不同组织差异甲基化区域 (Differentially Methylated Regions, DMRs) 的读段嵌入。预训练后，来自免疫相关组织（如胸腺、脾脏）的读段明显地聚集在一起，并与其他组织的读段分开。这表明 MethylQUEEN 具备了零样本（zero-shot）区分组织来源的能力。
图 e, f: 作者评估了模型在甲基化状态插补 (imputation) 任务上的性能。利用 PacBio HiFi 测序数据作为基准，MethylQUEEN 能够以极高的准确率（AUROC > 0.92）恢复被掩盖区域的甲基化状态（甲基化 vs. 未甲基化），显示了其强大的上下文理解和预测能力。

3.2 在序列分类任务中的卓越性能

作者定义了一种名为甲基化到组分 (Methylation-to-Components, M2C) 的新任务，即直接将携带甲基化信息的测序读段（reads）归类到其来源组分（如癌组织 vs. 正常组织）。

图 b, d: 作者在六个不同的 M2C 任务上对 MethylQUEEN 与其他模型（MethylBERT, DISMIR, AutoGluon）进行了基准测试。这些任务涵盖了癌症检测（HCC vs. 正常肝脏）、液体活检（胎盘 vs. 白细胞）、淋巴瘤分类（DLBCL vs. 正常B细胞）、细胞衰老和白细胞解卷积等多种场景。结果显示，MethylQUEEN 在所有任务中均取得了最佳性能，AUROC 最高提升了6.03%。
图 c: 以 DLBCL 与正常B细胞的分类任务为例，UMAP 可视化清晰地展示了在微调后，来自两个组分的读段被完美地分离开来。
图 e：在 DLBCL 分类中，一个高注意力 CpG 位点位于 MSX1 基因启动子区域，其上下文序列恰好匹配 c-Myc 转录因子（TF）的基序（motif），暗示该位点可能通过影响 c-Myc 的结合来调控 MSX1 的表达。在 FOXQ1 基因内，另一个高注意力 CpG 位点的上下文序列不匹配任何已知的 TF 基序，这可能是一个潜在的新型调控基序。
图 f, g: 对模型注意力分数的统计分析发现，甲基化相关令牌 (methyl-tokens) 受到的关注度显著高于其他 token，这直接证明了模型对甲基化信息的敏感性。此外，作者还发现，即使一个读段只包含一个 CpG 位点，模型也足以进行有效分类，增加 CpG 位点数量对性能的提升有限。
图 h: 在鲁棒性测试中，对非 CpG 位点引入随机突变后，模型的性能仅有轻微下降，显示了其在 M2C 任务中的稳健性。
图 i: 最引人注目的是，作者将模型应用于临床肝癌（HCC）血浆样本的分析。结果显示，MethylQUEEN 能够显著区分早期（Stage I）HCC 患者与健康人的血浆样本，预测的肿瘤负荷（tumor burden）存在显著差异。这突显了模型在极具挑战性的早期癌症液体活检中的巨大应用潜力。

3.3 准确预测基因表达水平并识别关键调控位点

作者进一步探索了利用甲基化信息预测基因表达（Methylation-to-Expression, M2E）的可行性。

图 a, b: 作者使用基因转录起始位点 (TSS) 上下游 ±2 kb 区域的甲基化序列来预测基因的表达水平（TPM）。在测试集上，预测的 TPM 值与真实的 RNA-seq 数据高度相关，Pearson 相关系数（PCC）达到0.966，R-squared 值为0.934。
图 c, d: 分析发现，平均表达量更高、表达变化更大的基因更容易被准确预测。
图 e, f: 为了提高模型的泛化能力，作者整合了来自 ENCODE、NIH Roadmap 和 GTEx 的多组织数据来扩充训练集。经过扩充训练后，模型能够准确地预测与 DNA 甲基化相关的关键基因的表达水平，并且预测的基因表达谱与 RNA-seq 检测到的真实图谱高度相似。
图 g: 作者通过分析模型在预测过程中的注意力分数，来识别对基因表达调控至关重要的 CpG 位点。以与阿尔茨海默病相关的 ABCA7 基因和与癌症相关的 KLHL6 基因为例，模型在进行表达预测时，会将更高的注意力集中在少数几个特定的 CpG 位点上。这些高注意力位点很可能是在近端甲基化调控中起决定性作用的关键位点，未来可能成为精准表观遗传编辑的治疗靶点。