摘要

在这项研究中，我们提出了iDNA-ABF，一个多尺度的深度生物语言学习模型，能够仅基于基因组序列对DNA甲基化进行可解释的预测。基准比较表明，我们的iDNA-ABF在不同甲基化预测方面的性能优于最先进的方法。重要的是，我们展示了深度语言学习的力量，它可以从背景基因组中获取顺序和功能语义信息。此外，通过集成可解释的分析机制，我们很好地解释了模型所学到的内容，帮助我们建立了从发现重要的序列决定因素到深入分析它们的生物学功能的映射。

关键词：DNA甲基化、深度学习、多尺度信息处理、可解释性分析

背景

DNA甲基化对胚胎发育具有重要意义，在基因沉默、防止虚假重复元件活性、有丝分裂过程中基因组稳定性和亲本来源印记方面发挥重要作用[1]。此外，环境和衰老引起的DNA甲基化模式的改变可能有助于疾病的发展，特别是癌症[2，3]。目前，5-甲基胞嘧啶(5mC)、N6-甲基腺苷(6 MA)和4-甲基胞嘧啶(4mC)是三种主要的DNA甲基化类型，根据核苷酸的类型、添加的分子类型和核苷酸内的修饰位置来命名[4]。不同的甲基化有不同的作用机制。例如，其中，5mC是通过与胞嘧啶(C)的第五位甲基结合而产生的。它与转录抑制有关，因此与基因组印记和X染色体失活等经典表观遗传现象有关[5]。6 mA，通常在腺苷(A)的第六位甲基化，在染色体复制、细胞防御、细胞周期调节和转录中起着关键作用[6]。作为另一种重要的表观遗传修饰，4mC保护宿主DNA不受限制酶的降解，纠正原核生物的DNA复制错误，控制原核生物的DNA复制和细胞周期[7]。因此，DNA甲基化鉴定是揭示其作用机制的基础。

DNA甲基化可以通过下一代测序(NGS)方法来实验确定，例如全基因组亚硫酸盐测序(WGBS)[8]或简化代表性亚硫酸氢盐测序(RRBS)[9]。这些技术可以在核苷酸水平上确定DNA甲基化的全球基因组分布，并为DNA甲基化相关的下游任务分析提供黄金标准数据集。然而，使用传统实验技术检测DNA甲基化通常是昂贵和耗时的[10]。此外，由于短读测序，亚硫酸氢盐测序不能分析重复基因组区域的DNA甲基化[11，12]。因此，最近的研究更多地集中在开发计算方法，特别是基于机器学习的方法，以直接使用基因组序列来检测DNA甲基化。这些方法将DNA甲基化识别作为一项二元预测任务，并训练机器学习模型来区分真正的甲基化位点和非甲基化位点。

在过去的几十年里，一系列使用传统机器学习或深度学习的基于序列的方法被很好地开发用于DNA甲基化的预测。以4mC甲基化预测为例，Tang等人提出了DNA4mC-LIP，这是一种集成学习方法，通过线性集成策略将现有的6个预测因子组合起来进行预测[13]。DeepTorrent[14]是一个基于深度学习的预测器，它集成了初始模块、注意模块和迁移学习来提高4mC站点的预测性能。MM-6mAPred[15]作为6mA位点的预测，利用了相邻核苷酸之间的跃迁概率，基于马尔科夫模型。为了简化模型构建，SNNRice6mA[16]使用卷积神经网络(CNN)构建了一个简单而轻量级的深度学习模型，以识别水稻基因组中的6mA位点。后来，Li等人提出了Deep6mA[17]，这是一种CNN和Long Short- Term Memory (LSTM)的混合深度学习网络，具有更准确的6mA预测。BERT6mA[18]是一个类似的模型，但使用变压器来建立预测模型，展示了自然语言处理技术在6mA预测中的应用的有效性。对于5mC位点检测，iPromoter-5mC通过全连接层[19]融合了几个预测单热编码序列的模型的结果。BiLSTM-5mC主要利用双向长短期记忆(Bidirectional Long - Short-Term Memory, BiLSTM)提取由核苷酸性质和频率编码的序列特征，用于5mC预测[20]。然而，大多数现有的方法只能区分一种单一类型的DNA甲基化。它们很难推广到其他甲基化类型。iDNA-MS[21]是第一个机器学习预测器，它是为不同物种的不同甲基化的通用检测而设计的。iDNA-MS利用人工特征，如k元组核苷酸频率成分和单核苷酸二进制编码，以及支持向量机(SVM)和随机森林(RF)等传统机器学习算法。iDNA-MS的缺点是特征设计需要大量的先验知识，同时在不同甲基化预测任务之间缺乏适应性。为了解决这个问题，在我们之前的工作中，我们设计了一个深度学习模型，即iDNA-ABT[22]，它使用了来自Transformers [23] (BERT)的双向编码器表征架构，自动和自适应地学习可区分的特征，并对不同物种的不同甲基化类型做出相对准确的预测。

如上所述，越来越多的研究尝试探索深度学习在DNA甲基化预测中的潜力，在提高预测性能方面也取得了一定的进展。然而，现有的深度学习预测器还没有充分探索特征表示学习的力量，特别是在发现对阐明DNA甲基化机制很重要的关键序列模式方面。这也导致深度学习模型解释性差，无法挖掘出基于序列的模型在DNA甲基化预测中的重要影响。另一方面，现有的方法未能回答其他重要问题:(1)背景基因组序列是否包含额外的可区分的信息，可以指导DNA甲基化的发展;(2)DNA甲基化的发生是否存在跨物种或细胞系序列模式的守恒性和特异性，从计算的角度来看也是一个关键问题。

随着自然语言处理技术的发展，一些先进的技术如BERT[24]能够充分地挖掘和学习自然语言文本中的高潜在语境信息。受此启发，我们将基因组序列视为“生物学文本”，并将不同尺度的顺序决定因子视为不同的“生物学词”。因此，我们提出了一种多尺度的生物语言学习模型——iDNA-ABF，成功地构建了从自然语言到生物语言的映射，以及从甲基化相关序列决定因素到它们的功能的映射。具体而言，我们引入了一个经过大规模基因组序列预训练的模型来学习生物上下文语义，并提出了一种多尺度处理策略来从不同尺度捕获鉴别甲基化信息。我们进一步利用对抗训练和迁移学习来改善预测性能，增强模型的鲁棒性。在不同甲基化和物种的17个数据集上的基准测试结果表明，我们的模型显著优于最先进的基于序列的方法。重要的是，我们的模型通过探索基于注意机制的局部序列特征，在序列水平上提供了可解释的预测和分析。结果表明，我们的模型能够准确且自适应地定位与甲基化密切相关的序列区域，这表明可能存在参与细胞进程功能调节的“生物语言语法”。

结果

提出的iDNA - ABF优于最先进的方法

为了评估我们提出的iDNA-ABF的性能，我们将其与四个最先进的预测器进行了比较，包括iDNA-ABT, iDNA-MS, BERT6mA和Deep6mA。在这四个预测因子中，前两个(iDNA-ABT和iDNA-MS)是不同甲基化预测的通用预测因子，而另外两个(BERT6mA和Deep6mA)最初是为6mA位点预测而设计的。之所以包含两个6mA预测因子进行性能比较，是因为它们是基于深度学习的最先进的预测因子。此外，他们的模型是灵活的，可以很好地扩展到其他甲基化预测，如5hmC和4mC，而不仅仅是6mA。所有比较的预测因子分别在不同物种和不同甲基化类型的17个训练数据集上进行训练，并在相应的独立测试数据集上进行评估(详情请参阅“数据集”部分)。ACC和MCC的评价结果分别如图1A和B所示。其他指标(如SN和SP)的详细结果见附加文件(附加文件1:表S1)。如图1A和B所示，我们的模型在17个数据集中的15个(附加文件1:表S1)上优于四个现有预测器，只有两个例外- 5hmc_m。肌肉和6mA_A。在Thaliana中，我们的模型实际上也与最好的预测因子相比较。具体而言，我们的模型在所有数据集上的平均ACC比排名第二的两个预测因子iDNA-ABT和BERT6mA分别高1.34%和3.73%。特别是在三个数据集上(4mC_C。equisetifolia 4 mc_s。酿酒酵母和6mA_S。而我们的iDNA-ABF在ACC上分别领先3.28-14.75%、1.88-3.59%和1.48-4.23%，比现有预测因子表现更好，幅度较大。在MCC方面也观察到类似的结果。为此，结果表明，我们的iDNA-ABF优于最先进的DNA甲基化通用预测方法。更重要的是，在三种甲基化类型下，它在物种间表现出了稳健的性能。

为了验证我们模型的稳健性，我们进一步说明了四个数据集(4mC_C)上预测因子的ROC和PR曲线(4mC_C.equisetifolia, 5hmC_M.musculus, 6mA_C. equisetifolia, and 6mA_F.vesca)，分别见图1C-F。我们可以看到，我们的iDNA-ABF在所有四个数据集中都具有最高的AUC和AP。具体而言，我们的模型在四个数据集上的平均AUC和AP值与其他预测因子相比分别增加了约1.39-2.81%和0.1-13.8%。结果进一步证明了我们的模型在DNA甲基化预测任务中的稳健性能。其他数据集上的ROC和PR曲线可以在附加文件中找到(附加文件1:图S1和图S2)。为了直观地讨论为什么我们的iDNA-ABF比其他方法表现得更好，我们进一步可视化了我们的iDNA-ABF和第二好的预测器iDNA-ABT在上述四个数据集(4mC_C。equisetifolia 5 hmc_m。, 6 ma_c。equisetifolia和6mA_F.vesca)使用统一流形近似和投影(UMAP)[25]，这是一种广泛使用的可视化工具，通过降维揭示基本数据特征。注意，其他数据集上的UMAP可视化结果可以在附加文件中找到(附加文件1:图S3)。图1G和H分别展示了我们的iDNA-ABF和iDNA-ABT的特征空间分布，其中每个点代表每个样本;甲基化位点(阳性样本)用红色标注，非甲基化位点(阴性样本)用蓝色标注。从图1G中可以看出，我们的模型将正负样本分离得很清楚，每一类都聚在一起而不是分散，而在图1H中，iDNA-ABT特征空间中的正负样本分布几乎是相连的，不容易围绕每一类的边界。通过对比图1G和H，我们发现这两类在我们的iDNA-ABF的特征空间中比最先进的iDNA-ABT分布更清晰。这表明，我们的模型从不同类别的样本中学习到更好的特征表示，这可能是由于我们在模型构建中预先训练好的模型，帮助我们从数百万个背景基因组序列中捕获更多高潜在的上下文语义信息。

图1 iDNA-ABF与现有方法的性能比较。A和B分别代表我们提出的iDNA-ABF和其他现有方法包括iDNA-ABT、dna - ms、BERT6mA和Deep6mA在17个基准独立数据集上的ACC和MCC值。C我们提出的iDNA-ABF和其他现有方法在5hmC_M.musculus中的ROC和PR曲线。D本文方法与现有方法在4mC_C.equisetifolia中的ROC和PR曲线。E在6mA_C.equisetifolia中，我们提出的iDNA-ABF和其他现有方法的ROC和PR曲线。F我们的iDNA-ABF和其他已有方法在6mA_F.vesca中的ROC和PR曲线。G和H为5hmC_M中iDNA-ABF和iDNA-ABT的特征空间分布(UMAP可视化)。骶4 mc_c。equisetifolia 6 ma_c。equisetifolia和6mA_F。分别vesca。阴性(红色)和阳性(蓝色)分别代表非甲基化和真甲基化样本。I在17个基准独立数据集上进行对抗训练和不进行对抗训练的模型的MCCs和ACCs;图中的每个点代表每个数据集。J在5hmC_M上使用对抗训练和不使用对抗训练的模型学习曲线。和6mA_F.vesca

对抗性训练增强了iDNA‐ABF的预测性能和稳健性

对抗训练是我们iDNA-ABF的重要组成部分。为了研究对抗训练的有效性，我们将原始的iDNA-ABF模型与没有使用对抗训练的模型进行了比较。17个独立数据集的结果如图1I所示，其中每个点代表每个数据集。可以看出，我们的原始iDNA-ABF(有对抗性训练)通常比没有对抗性训练的表现更好。具体而言，通过引入对抗性训练，可以在17个数据集中的14个和17个数据集中的15个上分别观察到ACC和MCC的性能改进。这表明对抗训练可以提高预测性能。其他指标(SN、SP和AUC)的结果可以在附加文件中找到(附加文件1:图S4)。此外，为了直观地展示对抗训练在模型优化中的重要性，我们进一步分析了训练过程中的学习曲线。图1J显示了在两个数据集(5hmC_M.)上进行对抗训练和不进行对抗训练的模型曲线。musculus和6mA_F.vesca)，从数据集中随机选择。从图1J中可以看出，进行了对抗训练的模型比没有进行对抗训练的模型获得了更低的测试损失，但损失减少率下降速度比没有进行对抗训练的模型慢。此外，使用对抗训练的模型在训练过程的后期保持较低的测试损失，而没有对抗训练的模型逐渐开始过拟合，这表明对抗训练增强了我们的模型在DNA甲基化预测中的鲁棒性。

我们的iDNA‐ABF在序列水平上揭示了物种间的甲基化守恒

为了研究不同物种甲基化序列模式是否保守，我们首先使用Lifemap[26]构建了相同甲基化类型的不同物种的进化树。对于4mC甲基化，图2A展示了四个物种的进化关系。可以清楚地看到，Fragaria vesca和木麻黄(Casuarina equisetifolia)属于进化分类学，属于常见的Fabids，而另外两个种属于Saccharomyces。一个有趣的观察是，我们的模型在具有进化分类的物种中表现出类似的性能。在蕨叶和木樨中，模型的ACCs分别为0.852和0.858;而另两组的ACCs分别为0.743和0.723。接下来，我们利用基于概率的motif可视化工具kplogo[27]进一步分析了四个物种的甲基化顺序模式。图2B显示了两个进化接近的物种(F. vesca和C. equisetifolia)的顺序模式，而图2C显示了另外两个物种的顺序模式。从图2B中，我们可以看到，该物种的甲基化序列区域非常相似，特别是CG含量丰富。从图2C可以观察到其他两个物种的相似结果。对于6mA甲基化，我们也发现了与4mC甲基化相似的结论(附加文件1:图S5)。总的来说，结果表明，甲基化序列模式在具有进化分类的物种中可能是保守的，从而有助于相似的预测性能;另一方面，在进化关系较远的物种中，甲基化模式会有很大的不同。

接下来，我们进一步研究了模型的跨物种表现，以研究它们之间的相互关系;也就是说，我们在一个物种上训练我们的模型，在另一个物种上进行测试。为了避免学习不足的问题，我们只在大数据集上训练模型，在小数据集上测试。从图2D可以看出，在进化分类学内的跨种性能明显优于没有进化分类学的跨种性能。结果进一步表明，序列水平的甲基化守恒与进化分类具有正相关关系。

图2跨物种甲基化守恒与模型精度的关系。4个物种在4mC数据集的分类树及其准确性。准确率较高的两种(F. vesca和C. equisetifolia)以红色表示，准确率较低的两种(S. cerevisiae和Tolypocladium)以蓝色表示。B鸢尾花和木犀花的图案logo分析。C酿酒酵母(S. cerevisiae)和Tolypocladium的motif logo分析。图2跨物种甲基化守恒与模型精度的关系4个物种在4mC数据集的分类树及其准确性。准确率较高的两种(F. vesca和C. equisetifolia)以红色表示，准确率较低的两种(S. cerevisiae和Tolypocladium)以蓝色表示。B鸢尾花和木犀花的图案logo分析。C酿酒酵母(S. cerevisiae)和Tolypocladium的motif logo分析。D 4mC数据集跨物种验证的准确性热图

多尺度顺序设计的选择更适合于阐明甲基化机制

在我们的模型中，我们提出了一种多尺度的信息处理策略，通过使用不同的k-mer来表示不同的“生物词”来进行特征表示学习。因此，我们首先验证了单尺度k-mers如何影响我们模型的预测性能。我们比较了不同的k-mer，从3-mer到6-mer。对比结果如图3A所示，我们可以看到不同的k-mer在不同的数据集上确实有各自的优势。没有观察到一致的结果。这可能是因为甲基化序列区域在不同物种和甲基化类型的长度上有所不同。因此，使用单尺度序列模式进行特征表示不能自适应和充分地捕捉甲基化的固有特征。为了解决这一问题，我们将不同尺度的k-mer作为模型输入，如3-mer + 6-mer, 4-mer + 6-mer和5-mer + 6-mer，并将其性能进行比较，如图3B所示。可以看出，与单尺度k-mer(即3-mer和6-mer)相比，多尺度k-mer集成(即3-mer + 6-mer)提高了模型性能。其中，3-mer和6-mer融合模型在所有数据集上的平均ACC达到了85.95%，分别比3-mer和6-mer模型提高了2.53和1.01%。这表明来自不同尺度的信息是相互补充的，以学习更好的特征表示。

图3多尺度信息处理的可解释分析。A单一尺度的比较，分别包括3-mer, 4-mer, 5-mer和6-mer。B多尺度组合的比较。C注意图，说明在一个随机选择的序列上以3个mer比例捕获的信息。两个子图分别显示了训练前后捕获的信息的变化。D注意图，说明在6 mer尺度捕获的信息。E-G我们的模型在三个物种中学习到的基元的可解释插图，涵盖了三种甲基化类型，包括4mC_Tolypocladium, 5hmC_H。和6mA_C。分别equisetifolia。左图通过从0到1的热图清晰地显示了模型更关注的区域。得分越接近1，颜色越深，模型考虑的区域越重要。通过比较我们的iDNA-ABF学习基序和STREME基序，使用TOMTOM计算p值。STREME中的p值通过单侧二项式检验计算。灰色虚线锚框内的图案被提取出来进行配对比较

接下来，我们进一步研究了为什么使用多尺度k-mer积分更适合于鉴别信息捕获。为此，我们利用注意机制来直观地解释我们的模型从3-mer和6-mer两个连续标度学习到的信息。我们将两个尺度的注意热图分别可视化在图3C和图D中。注意热图中的元素表示序列中两个位置的相关程度。图3C显示了我们的模型在3 mer尺度训练前后学习到的信息。我们可以看到，与初始模型相比，训练后的注意力机制更集中在热图的对角线上。这表明，与训练前相比，我们的模型学习了更多的局部判别信息。类似地，图3D展示了我们的模型在另一个连续尺度- 6-mer训练之前和之后学习到的信息。相比之下，该量表更侧重于培训后的全球信息。为此，我们可以得出结论，不同规模的顺序模式同时学习局部和全局信息，这可能是性能改进的补充。

为了明确哪个序列区域对甲基化预测最重要，我们从3个不同DNA甲基化类型的物种中随机选取了3个序列，并应用注意机制从这些序列中识别关键区域。从图3E-G(左)可以看出，对于每个序列，我们的模型在不同的顺序尺度下识别出了不同的区域。这进一步证实了不同的尺度捕获不同的重要信息。对于那些被识别的区域，我们进一步提取并使用注意力分数将相应的主题可视化。图3E-G(右)分别显示了我们的iDNA-ABF学习到的基序和传统工具- streme[28]发现的基序。正如我们所看到的，我们学习的母题(用灰色窗口突出显示)几乎与每个物种的STREME母题相匹配。为了定量比较motif相似度，我们采用TOMTOM[29]来计算两个motif的相似度，用p值来衡量。p值越低，motif一致性越高。如图3E-G所示，我们的motif与STREME的motif高度相似，这表明我们的模型可以学习保守的序列特征。

我们的iDNA‐ABF充分探索了人类细胞系5mC预测中的基因组信息

在本节中，我们分析了我们的iDNA-ABF如何在人类细胞系中执行甲基化预测。由于5mC是人类基因组中研究最充分的甲基化类型之一，我们选择了5mC甲基化来执行我们的方法。因此，我们构建了三个新的5mC数据集，分别对应于三个人类细胞系GM12878、K562和HepG2。数据集的详细信息可以在“数据集”一节中看到。

首先，我们讨论了甲基化序列区域长度对5mC甲基化预测的影响。因此，我们为每个细胞系构建了4个5mC数据集，每个5mC序列的长度分别为11、41、71和101 bp(碱基对)。数据集的详细信息汇总在附加文件中(附加文件1:表S2，表S3和表S4)。图4A显示了三种细胞系中不同序列长度下模型性能的变化。首先，随着序列长度的增加，模型的性能显著提高，说明较长的序列为模型带来了额外的基因组上下文信息。当长度为71 bp时达到峰值。之后，模型性能逐渐下降。值得注意的是，使用11 bp长的序列训练的模型表现出极差的性能，ACC约为55%。原因是，在11个碱基范围内的甲基化中心区域在阴性和阳性样品之间非常相似。这进一步说明甲基化与甲基化区上下游高度相关。

图4我们的模型在三种人类细胞系上的5mC预测性能。A不同序列长度的人类细胞系数据集的ACC和AUC结果。三种人类细胞系不同序列长度下序列数据、ChIP-seq数据和ChIP-seq数据集成的B-D性能。E我们的模型预测的5mC分布和随机选择的基因组区域(Chr1: 187000 - 192000, GRCh38)的WGBS注释。注意，5mC分布来源于HepG2细胞系

众所周知，5mC甲基化是目前研究最深入的甲基化类型之一，有很多NGS数据支持，如ChIP-seq数据、ATAC-seq数据等。一个有趣的问题是，将NGS数据与序列数据集成是否有助于更准确的预测。为此，我们选择了两个组蛋白修饰数据H3k4me3和H3k36me3，这两个数据被报道与5mC[31]密切相关。我们分别使用(1)序列数据，(2)ChIP-seq数据和(3)序列+ ChIP-seq数据在三个细胞系上训练和测试模型。对比结果如图4B-D所示。我们可以看到，序列数据训练的模型比ChIP-seq数据训练的模型有明显的提高，在不同序列长度的三种细胞系中，平均ACC、AUC和MCC分别领先10.4、10.9和21.1%。当将ChIP-seq数据与序列数据相结合进行模型训练时，所有性能指标都得到了进一步的改进，与使用序列数据训练的模型相比，平均ACC、AUC和MCC分别提高了3.8、5.2和8.1%，表明ChIP-seq数据和序列数据对于改进的5mC预测是相互补充的。

iDNA‐ABF在5mC甲基化检测中的应用

考虑到实际应用场景，我们的iDNA-ABF在全基因组尺度上检测5mC分布的性能非常重要。因此，我们基于HepG2训练的iDNA-ABF模型，预测了人类基因组(GRCh38)中一个5k-bp长的基因组区域(Chr1: 187000 - 192000)的甲基化概率。预测过程如下。首先，我们使用一个71 bp长的窗口来筛选区域。其次，选取满足以下两个要求的序列:(1)以基C为中心，(2)以CPG模式为中心。最终，得到的序列被提交给我们的iDNA-ABF进行预测。我们的模型给出了每个候选站点的预测置信度。

图4E说明了两种HM数据分布，分别是我们模型预测的5mC分布和WGBS标注的真实5mC分布。从图4E可以看出，我们预测的5mC分布与真实的5mC站点分布大致重叠。预测的5mc与两个HM数据基本吻合，说明我们的预测具有功能意义。值得注意的是，我们发现我们的模型识别了一些WGBS无法识别的区域(图4E中蓝色框)，但它们与H3K4me3数据的信号匹配得很好。这意味着我们的模型可能会发现潜在的新功能区域。虽然我们的模型也产生了一些假阳性，但从顺序仓的角度来看(这里，我们将100-bp区域视为一个仓)，预测的5mC区域分布与真实的5mC区域分布几乎相同。结果至少表明，我们的模型可以很好地定位5mC区域。这对甲基化研究也有帮助。

我们的iDNA‐ABF在未见过的人类细胞系的5mC预测中具有强大的性能

为了分析iDNA-ABF在未见细胞系中的预测性能，我们进行了跨细胞系验证。具体来说，我们在一个细胞系上训练我们的模型，并在另一个细胞系上评估它。图5A显示了四个指标的热图结果，分别包括ACC、MCC、SN和SP。纵轴为训练细胞系，横轴为测试细胞系。如图5A所示，我们的模型在跨细胞系验证下获得了相对稳定的ACC和MCC。此外，我们还可以看到，在K562上进行评估时，我们在GM12878上训练的模型获得了最高的SN，与在HepG2上训练的模型相比，获得了16%的相对改进。为了更好的解释，我们介绍了三种人类细胞系甲基化中心区域的概率分布分析。图5B和C分别显示了三种细胞系中阳性和阴性样本的概率分布。一方面，从图5B可以看出，K562的正motif logo在−1到1的位置上与GM12878的相似度高于HepG2。另一方面，我们从图5C中观察到，K562的负motif标志与GM12878的正motif标志在−1到1的位置相同，这可以解释我们在GM12878上训练的模型在K562上测试时SP最低。此外，我们在图5A中发现，我们在跨细胞系验证热图中的模型在SP方面表现不太好。这可能是由于三种细胞系之间的负基序标识存在很大差异。为此，通过跨细胞系的验证结果，我们可以得出结论，我们的模型具有稳健的性能，即使是看不见的细胞系。进一步探讨了模型的应用价值。

图5未见细胞系的5mC预测性能。A不同指标的跨细胞系验证热图，分别包括ACC, MCC, SN和SP。在三种人类细胞系阳性数据集的中心序列区域分别有B Motif标识。C分别在三个人类细胞系的阴性数据集中中心序列区域的Motif标识

我们的iDNA‐ABF具有良好的转移学习能力，可以捕获甲基化序列模式的特异性

人类基因组中5mC甲基化主要发生在CpG模式的序列中;实际上，在少数情况下，甲基化也在CHH和CHG模式中被检测到(其中H = a, C，或T)。为了找出不同的甲基化序列模式是否相互关联，我们分别为这三种细胞系构建了额外的CHG和CHH数据集。值得注意的是，CpG数据集的序列数量远远多于CHG或CHH数据集。数据集的详细信息见附加文件(附加文件1:表S5)。

为了检验我们的模型在检测不同甲基化模式时是否具有良好的迁移学习能力，我们首先在CpG数据集上预训练一个模型，并在CHG或CHH数据集上对其进行微调，得到另一个模型，称为“迁移学习模型”。此外，我们还直接用CHG或CHH数据集训练了一个模型进行比较，称为“基线模型”。然后使用CHG或CHH数据集的相同测试数据集对这两个模型进行评估。在两个数据集上的性能分别如图6A和B所示。我们可以看到，“迁移学习模型”的性能始终优于基线模型，三个细胞系的平均AUC和AP分别提高了3.1和3.3%。结果表明，该模型具有良好的迁移学习能力;预训练机制可以从特定模式中带来额外的判别信息，有利于目标模式的预测，从而提高预测性能。

图6迁移学习结果及模型分析。A HepG2细胞系CHG数据集上基线模型和迁移学习模型的ROC和PR曲线。请注意，基线模型是用CHG数据集训练的，而迁移学习模型是用CpG数据集预训练的，并在CHG数据集上进行微调。B HepG2细胞系CHH数据集中基线模型和迁移学习模型的ROC和PR曲线。请注意，基线模型是用CHH数据集训练的，而迁移学习模型是用CpG数据集预训练的，并在CHH数据集上进行微调。C基元分别从基线模型、预训练模型和迁移学习模型三种模型中学习。D使用我们学到的母题对topenrichment进行搜索的结果。MM数据库

为了深入解释迁移学习提高性能的可能原因，我们进一步分析和比较了从基线模型、迁移学习模型和仅在CpG数据集上训练的预训练模型中学习到的主题。在HepG2细胞系中motif比较结果如图6C所示。以第1、5、9个motif图为例，我们观察到迁移学习模型不仅保留了一些从预训练模型继承的CpG模式，而且还捕获了基线模型学习到的CHG模式的特异性。此外，迁移学习模型还可以发现一些新的模式，例如第7个motif图，它与基线模型在第3个motif图中没有相似的模式。

从我们的模型中学到的母题在生物学上是有意义的

接下来，我们进一步探讨了从上述三个模型(基线模型、预训练模型和迁移学习模型)中学习到的基序(或顺序模式)是否具有生物学意义。因此，我们对topenrichment进行了学习母题的搜索。MM，一个公共甲基化数据库[32]。有趣的是，从图6D中，我们发现我们的模型学习到的基序可以与数据库中的一些功能基序显著匹配，这些功能基序之前被报道与甲基化机制密切相关。结果表明，该模型能够准确地挖掘函数序列特征;另一方面，新发现的序列基序也具有生物学意义，表明该模型具有较强的学习不同序列模式之间功能语义的能力。

讨论

我们提出了一种仅基于基因组序列的生物语言学习来识别DNA甲基化的新方法——iDNA-ABF。我们的iDNA-ABF不仅能够相对准确地预测跨物种和跨细胞系的甲基化，而且还使用可解释的注意机制建立了从顺序水平到功能水平的映射，以深入研究DNA甲基化机制。

首先，我们研究了模型的预测性能，看看它的表现有多好，有多稳定。在17个涵盖多种物种的三种甲基化类型(4mC, 5hmC和6mA)的基准数据集上的实验结果表明，与目前最先进的基于序列的方法相比，我们的模型表现出始终优越和稳健的性能。烧蚀分析揭示了对抗训练在模型性能中的重要性。特别是，训练过程中的对抗训练减轻了大规模参数的影响，特别是对一些小数据集的影响，提高了模型在不同物种和甲基化类型之间的泛化能力。此外，我们还研究了序列长度对甲基化预测的影响。结果表明，随着输入序列长度的增加，模型的性能普遍提高，说明甲基化序列区域周围的上游和下游区域对于识别DNA甲基化位点至关重要。从序列的角度来看，它们可能包含一定程度的特异性信息，以帮助我们的模型区分甲基化位点和非甲基化位点;另一方面，我们采用DNABERT[33]进行模型构建，这是一个强大的自然语言学习模型，使用百万级基因组序列数据进行预训练。它使我们的模型能够从背景基因组中捕获更多的序列语义。特征空间可视化分析结果证明，与现有预测器相比，我们的模型学习了更多可区分的特征表示。有趣的是，通过将组蛋白修饰数据(如H3K4me3和H3K36me3)等ChIP-seq数据集成到我们的模型中，我们观察到性能进一步提高，这表明生物信号和序列数据对于改进的预测是互补的。该结果显示了序列数据在DNA甲基化预测和其他基因组功能分析方面的巨大潜力。我们可以想象，对于新的细胞系数据，NGS数据是有限的，我们无法训练出有效的模型。至少，这些分析为整合序列数据建立更精确、更稳健的模型提供了一种新方法。

其次，我们的iDNA-ABF的主要特点是我们为DNA甲基化预测提供了可解释的分析。现有的基于深度学习的方法的主要问题是，它们不能很好地解释为什么它们的模型对甲基化预测是有效的，因为深度学习是作为“黑箱”工作的。为了解决这个问题，我们对模型构建做了两个主要的改进。一是提出了用于模型训练的多尺度序列处理策略，二是引入了用于模型分析的注意机制。受自然语言学习中的分词技术的启发，我们利用多尺度序列处理策略，将序列模式的不同尺度(3mer和6mer)的DNA序列分割为“生物词”。此外，我们采用注意机制来解释我们的模型从不同尺度的“生物词”中学习到的信息。分析表明，多尺度策略能够从局部和全局两个层面带来更多的语义判别信息，有效地克服了单一尺度上的信息缺失和各尺度上的信息冗余。重要的是，不同的顺序尺度使我们的模型学习不同的主题。结果表明，我们从不同尺度学习到的motif与传统的motif发现工具streme的motif高度一致，表明我们的模型具有发现保守序列模式的能力。接下来，自然的问题是，我们的模型学习到的序列模式是否在生物学上有意义或与甲基化相关。为了回答这个问题，我们将我们的模型应用于人类细胞系5mC甲基化的预测。选择人类5mC甲基化的原因是它具有保守的甲基化顺序模式，如CpG、CHH和CHG;另一方面，得到了大量NGS数据的支持，便于进一步的功能验证分析。我们研究了模型的迁移学习能力，实验结果表明，从CpG甲基化模式中学习知识可以帮助改进对其他两种甲基化模式(即CHH和CHG)的预测。这表明我们的模型在学习不同甲基化模式的特异性方面有很强的能力。这些结果也暗示了我们的模型在发现其他很少发生的甲基化模式方面的潜力。重要的是，我们发现通过使用迁移学习，我们的模型可以学习一些新的基序，同时保持基序在原来的甲基化模式中。通过搜索我们学习到的基序与著名的甲基化数据库topenrichment。MM，我们发现我们的基序与一些已知的甲基化相关功能基序显著相似(见图6D)。这也证明了我们的模型可以在不同的甲基化模式下学习不同的生物语义信息。

最终，为了验证模型在实际应用场景中的性能，我们进一步将模型应用于人类基因组中5mC甲基化的检测。在随机选择的基因组区域的实验结果表明，我们的模型能够准确检测真实的DNA甲基化区域(由WGBS注释)。重要的是，我们的模型发现了一些潜在的甲基化区域，这些区域没有被WGBS检测到，但与甲基化相关的组蛋白修饰数据高度重叠(例如H3K4me3)。这证明了我们的模型在发现生物学上有意义的序列区域方面的强大能力。深度预训练模型可能帮助我们从数百万个背景基因组序列中学习功能语义。

结论

总之，我们提出的深度生物语言学习模型在DNA甲基化预测方面取得了令人满意的效果。重要的是，我们展示了深度语言学习在从背景基因组中捕获序列和功能语义信息方面的力量。此外，通过整合可解释分析机制，我们很好地解释了我们所学到的东西，帮助我们建立了从发现重要序列决定因素到深入分析其生物学功能的映射。然而，仍有很大的改进空间。例如，在甲基化预测模型的构建中，我们只考虑了甲基化位点周围的局部序列区域，其中的判别信息可能在一定程度上受到限制。研究[34]表明基因调控在基因组中存在一种长期的相互作用，如增强子-启动子相互作用。因此，探索远程序列整合信息如何影响DNA甲基化水平可能是未来工作的一个重要方向。

方法

数据集

不同物种数据集

严格的数据集对于训练有效和有前途的预测因子是至关重要的。为了用最先进的方法进一步评估我们提出的方法，我们选择了最初由dna - ms[21]提出的相同的基准数据集。数据集包括三种主要的DNA甲基化类型，包括总共17个数据集。在17个数据集中，C.equisetifolia (4mC_C.equisetifolia)、F.vesca (4mC_F.vesca)、S.cer - evisiae (4mC_S.cerevisiae)和Ts.SUP5-1 (4mC_Ts.SUP5-1)属于4mC。6mA包括拟南芥(Arabidopsis thaliana, 6mA_A.thaliana)、秀丽隐杆线虫(Caenorhabditis elegans, 6mA_C.)。e- gans)、木麻黄(Casuarina equisetifolia, 6mA_C.equisetifolia)、果蝇(Drosophila melanogaster, 6mA_D.;melanogaster)， Fragaria vesca (6mA_F.vesca)，智人(6mA_H.sapiens)， Rosa chinensis (6mA_R.chinensis)， Saccharomyces cerevisiae (6mA_S.cerevisiae)， Tolypocla- dium sp SUP5-1 (6mA_Tolypocladium)，嗜热四膜虫(6mA_T. sapiens)。嗜热菌)和米白单胞菌PV。Oryzicola (Xoc) BLS256 (6mA_Xoc.BLS256)。此外，还获得了来自两个物种的两个5hmC数据集，包括H.sapiens (5hmC_H.sapiens)和M.musculus (5hmC_M.musculus)。需要注意的是，采用CD-HIT[35]程序，阳性和阴性数据集的长度均为41个碱基对(bp)，序列同一性小于80%，如图7A所示。来自17个物种的训练数据集和验证数据集的详细信息在附加文件(附加文件1:表S6)中给出。

人类细胞系数据集

从ENCODE (ENCSR765JPC, ENCSR890UQO和ENC- SR786DCL)[36]中收集了三种人细胞系(K562, GM12878, hepG2)的5mC甲基化数据，提供了通过全基因组亚氢盐测序(WGBS)实验的三种甲基化模式(CpG, CHG和CHH)的位置信息。为了构建高质量的数据集，阳性样本保留100%甲基化和10 - 200倍测序覆盖率的甲基化位点，阴性样本选择0%甲基化和0测序覆盖率的甲基化位点。利用GENECODE GRCh38的注释进一步绘制了位于启动子和基因体区域的加工甲基化位点。启动子区域被定义为基因转录起始位点(TSS)上游1000 bp的区域。每个单元行处理数据集的编号显示在附加文件(附加文件1:表S7，表S8和表S9)中。为了评估序列长度对甲基化预测的影响，分别从GRCh38中提取甲基化位点两侧长11、41、71和101 bp的DNA序列。与不同物种的数据集相似，使用CD-HIT程序的人类细胞系数据集的序列同一性也小于80%。

图7所提出的iDNA-ABF概述。A显示了DNA甲基化数据集集合，其中属于三种主要DNA甲基化类型的不同数据集被重组为它们的训练数据集和独立数据集。我们的iDNA-ABF的整体架构在B-E中给出。B多尺度信息处理模块，分别利用分词器的3-mer和6-mer两个尺度对输入序列进行处理，自适应获得相应的嵌入。C BERT编码模块，利用BERT编码器提取高潜伏特征表示。D多尺度提取模块，生成基于多尺度嵌入的最终输出特征表示。E分类模块，集成二元分类概率值进行预测。F可解释分析的工作流程。简而言之，我们的模型使用注意机制从查询序列中提取和学习顺序母题

拟议iDNA‐ABF的描述

图7说明了我们的iDNA-ABF的整体架构。图7A显示了数据集收集过程，该过程在“数据集”一节中描述。iDNA-ABF的工作流程如图7B-E所示，主要由四个模块组成:(B)多尺度数据处理模块，(C) BERT编码器模块，(D)特征融合模块，(E)分类模块。预测过程描述如下。在模块B中(见图7B)，我们分别利用两种尺度的标记器(3-mer和6-mer)来处理输入序列并自适应地学习相应的嵌入。由于输入序列包含不同尺度的多面特征，我们设计了一个多尺度的体系结构，而不是使用一个简单的标记器，这可能会导致信息丢失。然后，在模块C中(见图7C)， iDNA-ABF分别使用BERT编码器提取经过标记化处理的不同嵌入。然后，iDNA-ABF在模块D中结合基于BERT输出的多尺度嵌入，生成最终的进化输出特征。之后，在模块E(图7E)中，模型使用全连接层来预测输入序列是否甲基化。值得注意的是，我们采用对抗性训练来增强模型的鲁棒性和防止早期过拟合，这可以分为两个部分:(1)对抗性每扰动，利用这种传播的交叉熵损失作为对抗性扰动返回网络;(2)对抗性优化，获得对抗性损失，用于反向传播和优化我们的模型。请注意，我们将如下描述这四个模块的详细信息。

多尺度信息处理模块

在我们的模型中，我们用k-mer表示标记DNA序列。通过这种方式，每个标记都由k个碱基表示，从而为每个核苷酸集成了更丰富的上下文信息。例如，给定的DNA序列“ATGGCTG”可以被标记为两个6-mers序列:ATGGCT和TGGCTG。不同的k会导致不同的令牌表示。在我们的工作中，我们将k设为3或6，从而获得两个尺度的标记表示。整个令牌表有4k+5个令牌，由k-mer的所有排列以及5个特殊令牌组成:[CLS]、[PAD]、[UNK]、[SEP]和[MASK]，分别代表分类令牌、填充令牌、未知令牌、分离令牌和掩码令牌。

基于BERT的编码器模块

BERT模型的预训练

BERT是[24]提出的第一个基于转换器的双向语言表示模型。由于BERT在多种大型语料库中具有强大的语言理解能力，它已被广泛应用于大量的自然语言处理任务中。为了更好地发挥BERT的作用，它通常会首先在一个大型背景相关语料库上进行训练，并进行两个预训练任务，即蒙面语言模型和下一句预测。这里我们使用预训练BERT模型dnbert[33]，使用与BERT基础相同的架构，它由12个变压器层组成，每层有768个隐藏单元和12个注意头。值得注意的是，由于DNA序列之间没有直接的语义逻辑，这个域预训练调整了序列长度，并使模型能够预测适应DNA序列的连续k个标记。此外，它还使用了类似于原始BERT的掩码语言模型技术。

BERT的编码过程

BERT是一种基于变压器的情境化语言表示模型，它已被应用于生物学的许多方面，并取得了许多杰出的成绩。BERT的基本组成部分包括多头注意机制、前馈网络和残差连接技术。为了捕获上下文信息，BERT在自我注意机制的基础上执行多头注意机制，其描述如下:

其中 $X\epsilon R^{L\times d_{m}}$ 为序列嵌入模块的输出。并分别根据线性层W Q, W K, W V∈Rdm ×dk，将X变换为查询矩阵Q∈RL×dk，关键矩阵K∈RL×dk，值矩阵V∈RL×dk，其中L为输入蛋白质序列的长度，dm为初始嵌入维数，dk为矩阵Q、K、V的维数。

由上述基本单元，多头注意机制可表示为:

其中W Q, W K, W V∈Rdm×dk为第i个头的查询、键和值线性变换层iii, h为头的个数。然后，将h个独立头与不同WiQ、WiK、WiV集合的结果进行拼接，利用线性转换层WO将多头注意的输出维数映射到嵌入模块的初始嵌入维数。整个过程执行L次，其中L表示层数。

特征融合模块

为了得到两个BERT部分的最终输出hM，我们将第一个比例输入格式层的输出hkmer1和第二个比例输入格式层的输出hkmer2通过一个维度融合门F结合起来。F由sigmoid激活函数完成，对两部分表示进行编码:

其中W1和W2是融合门的可训练参数。然后通过F生成特定分子hM的最终矢量表示输出:

分类模块

对抗训练[37]是一种新的正则化方法，用于提高分类器对小的，近似最坏情况扰动的鲁棒性。在这里，由于BERT具有较大的参数，我们使用这种策略来防止模型过拟合。在众多对抗训练方法中，我们使用了一种针对文本分类[38]的快速梯度方法(FGM)。交叉熵损失函数LCE用于训练输出模块，以提高预测性能。我们定义p为预测概率，y为真实标签，x为输入，θ为模型的参数，ε为一个附加参数。在应用这种方法时，对抗性训练在代价函数中增加了以下项:

其中r是对输入的扰动，θ是模型当前参数的常数集。反向传播算法不应该通过对抗实例构建过程来传播梯度，这意味着θ与Eq.(6)中的θ不一致。然后，在训练过程中，我们将Eq.(6)中的θ最小化，以获得对当前模型的最坏情况扰动radv。

在FGM方法中，我们将对抗性扰动应用于提取的序列嵌入，而不是直接应用于输入。为了定义词嵌入上的对抗性扰动，我们将k-mer的相关嵌入表示为s。然后我们定义s上的对抗性扰动radv为

为了训练鲁棒模型，我们基于式(6)中定义的对抗性扰动定义了新的对抗性损失，其公式如下:

其中N为批量大小的数量。在我们的工作中，对抗性训练是最小化基于随机梯度下降的交叉熵损失的Ladv。

性能指标

在这项研究中，我们用以下四个常用的指标来评估我们的iDNA - ABF和其他现有方法的性能:准确度(ACC)、马修斯相关系数(MCC)、灵敏度(SN)和特异性(SP)。这些指标的公式描述如下:

其中TP、FN、TN、FP分别代表真阳性、假阴性、真阴性、假阳性样本数。ACC和MCC都被用来衡量模型的整体性能。SN是指预测模型正确预测的真实甲基化样本的比例，SP衡量模型正确预测的非甲基化样本的比例。并利用ROC (receiver operating characteristic)曲线和PR (precision-recall)曲线[39]直观地评价了模型的整体预测性能。AUC为ROC曲线下面积，AP为PR曲线下面积，分别为[39]。它们被进一步用于定量度量模型的整体性能。总之，这些指标越高，模型就越好。

补充信息

在线版本包含补充材料，可在https://doi.org/10.1186/s13059-022-02780-1获得。

附加文件1:表S1。iDNA-ABF的性能和最先进的方法在17个基准数据集跨物种和甲基化类型。表S2。在三种人类细胞系中，我们的iDNA-ABF在不同甲基化模式下的性能随序列长度的不同而不同。表S3。利用ChIP- seq数据，我们的iDNA-ABF的性能随三种人类细胞系的不同序列长度而变化。表S4。使用ChIP-seq数据+序列数据的iDNA-ABF性能在三种人类细胞系中随着序列长度的不同而不同。表S5。在三种人类细胞系中，我们的模型的迁移学习性能随序列长度的不同而不同。表S6。17个不同物种甲基化类型的基准数据集的统计。表S7。三种人类细胞系中不同甲基化模式的原始数据的统计。表S8。用CD-HIT方法对三种人类细胞系进行序列相似性降阶后的数据统计。表S9。三种人类细胞系在不同甲基化模式下的训练和测试数据的统计。表S10。以不同尺度的性能比较作为模型在不同物种中的训练输入。表S11。与5mC方法在癌细胞系百科全书(CCLE)上的性能比较。表S12。iPromoter-5mC在三种人细胞系中的表现。表向。在17个基准数据集上训练模型参数。表S14系列。我们的iDNA-ABF用于SNP分类的性能。表S15。多任务模型的性能。表S16。基于ChIP-seq数据建立的回归模型的性能。图S1。基准数据集上的ROC曲线。图S2。基准数据集上的PR曲线。图S3。基准数据集上的UMAP可视化结果。图S4。在17个基准数据集上进行独立测试，计算了经过对抗训练和未经过对抗训练的模型的SN、SP和AUC。图S5。6mA数据集中11个物种的分类树及其准确性。图S6。在GM12878和K562细胞系的另外两个模式数据集(即CHG, CHH)上，我们的基线模型和我们的模型的ROC和PR曲线从CpG模式转移。图S7。基于ChIP-seq数据建立的回归模型中信号预测的回归结果。补充方法。

附加文件2.回顾历史

iDNA-ABF：DNA甲基化可解释预测的多尺度深度生物语言学习模型相关推荐

IDNA-ABF: DNA甲基化可解释预测的多尺度深度生物语言学习模型
IDNA-ABF:multi‑scale deep biological language learning model for the interpretable prediction of DN ...
DSML：深度子空间相互学习模型(用于癌症亚型预测)
目录背景描述模型架构数据表示 Single level omics data representation learning model(单级组学数据表示学习模型) Multi-level om ...
IJCAI 2019 | 为推荐系统生成高质量的文本解释：基于互注意力机制的多任务学习模型...
编者按:在个性化推荐系统中,如果能在提高推荐准确性的同时生成高质量的文本解释,将更容易获得用户的"芳心".然而,现有方法通常将两者分开优化,或只优化其中一个目标.为了同时兼顾二者, ...
Facebook：易于解释的神经元可能会阻碍深度神经网络的学习
选自 Facebook AI 博客作者:Matthew Leavitt.Ari Morcos 机器之心编译编辑:张倩.杜伟易于解释的神经元对于提升神经网络的性能来说是必要的吗?Facebook ...
易基因｜RRBS单碱基绘制580种动物的基因组规模DNA甲基化谱：Nature子刊
大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因. 2023年01月16日,奥地利科学院分子医学研究中心(CeMM)研究团队在<Nat Commun>杂志发表了题为" ...
易基因｜糖尿病大血管和微血管病变的表观遗传学调控机制（DNA甲基化+组蛋白修饰）：综述
大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因. 2022年12月09日,<Frontiers in Endocrinology>杂志发表了题为"Epigeneti ...
Nat. Mach. Intell. | 基于深度学习预测DNA甲基化位点
研究人员开发了一种预测DNA甲基化位点的机器学习算法可以帮助识别致病机制.该论文2020年8月3日发表在"Nature Machine Intelligence"上. 研究人员通过 ...
易基因｜综述：单细胞DNA甲基化分析方法全介绍及未来发展前景预测
大家好,这是专注表观组学十余年,领跑多组学科研服务的易基因. 2021年07月10日,<Biomolecules>杂志上发表一篇关于单细胞表观测序的综述文章,详细介绍了单细胞DNA甲基化 ...
精准医学：粪便DNA甲基化标志物在大肠癌分期及其先兆诊断中的应用｜早期筛查
易点评在精准医学系列文章中我们介绍过有关血液中甲基化标志物在大肠癌复发监测中的应用,今天则为大家带来一篇关于粪便甲基化标志物在大肠癌早期筛查中应用的综述.在这篇综述中作者调研了2020年1月之前发表 ...

iDNA-ABF：DNA甲基化可解释预测的多尺度深度生物语言学习模型

摘要

背景

结果