论文笔记：m6Acorr: an online tool for the correction and comparison of m6A methylation profiles

前言

论文发表期刊：《BMC BIOINFORMATICS》
期刊影响因子：3.242
期刊中科院分区：生物 3 区
m6Acorr服务器： http://www.rnanut.net/m6Acorr
校正函数及数据网址： https://github.com/emersON106/m6Acorr
补充材料网址： https://doi.org/10.1186/s12859-020-3380-6
作者简介：李建伟来自河北科技大学，黄燕来自北京大学

必备生物知识

*1、m6A定义
m6A中文名为N6-腺苷酸甲基化，指腺苷上的N6位发生甲基化所以称为N6-腺苷酸甲基化（m6A），目前科学家已经在RNA中鉴定了超过100种不同类型的碱基修饰行为。mRNA最常见的内部修饰包括了N6-腺苷酸甲基化（m6A）、N1-腺苷酸甲基化（m1A）、胞嘧啶羟基化（m5C）等。对于大热的m6A，截止当前，全球的科学家已经鉴定了参与m6A的许多酶，包括去甲基化酶、甲基化酶和甲基化识别酶等。N6-甲基化腺苷酸结构示意图如图一：

图1 N6-甲基化腺苷酸结构示意图
图一是已经甲基化的核糖核苷酸，确切地说叫N6-methyladenosine。一共分为2个大的结构。我们先来回顾下基础生物化学知识。左下角的是五碳糖，图2中a框部分也就是五碳糖的第二位C处的羟基发生脱氧就会变成脱氧核糖核苷酸（从RNA变成DNA）。图2中c框部分标注的，也就是第四位的C处通常会带有磷酸基，如果此处带有2个磷酸基团那么就叫ADP，如果带有3个磷酸基团那就是大名鼎鼎的ATP了。图2中b框部分通常就是我们所说的含氮碱基，这里特指腺苷酸（A）。三种与m6A有关的甲基化，如下图所示：

图2 m6A甲基化修饰和m6Am超甲基化修饰的概念

2、m6A甲基化加工过程
m6A这种甲基化修饰被证明是可逆化的，包括甲基化转移酶、去甲基化酶和甲基化阅读蛋白等共同参与。其中甲基化转移酶包括METTL3/14、WTAP和KIAA1429等，主要作用就是催化mRNA上腺苷酸发生m6A修饰。而去甲基化酶包括FTO和ALKHB5等，作用是对已发生m6A修饰的碱基进行去甲基化修饰。阅读蛋白主要功能是识别发生m6A修饰的碱基，从而激活下游的调控通路如RNA降解、miRNA加工等。，其基本过程如图三所示：

图3 m6A甲基化加工过程

当从DNA→RNA过程中，腺苷酸在甲基化酶METTL3、METTL14和WTAP等作用下在第六位N发生了甲基化修饰。这类酶我们称之为Writers，其中METTL3和METTL14形成杂络物（hetero complex），与WTAP（也叫做Fl(2)d）及其他factors如KIAA1429（也叫做Virilizer）一起对腺苷酸进行修饰。YTHDF家族包括YTHDF1、YTHDF2、YTHDF3以及酿酒酵母中的Mrb1基因、粟酒裂殖酵母中的Mmi1基因都是readers类蛋白。这些酶能够识别发生m6A甲基化的碱基，参与下游翻译、mRNA降解、加快mRNA出核速度等作用。具体信息如下表所示：

图4 writers和Readers信息表

3、m6A数据库
关于m6A测序的技术叫做meRIP-seq。这个测序结果的分析，类似于chip-seq。最后我们可以获得一种叫做peak的文件，这个代表m6A甲基化修饰在哪个位置。基于测序数据来创建的数据库一共有，四个。其中其中18年相同时间发表了两篇(MeT-DB, RMBase)，这两篇被引次数也是很高。说明两个数据库挺好用的。同时其实也能说明m6A在这几年研究的还是很多的呀。另外的一个Whistle是19年发表的，被引8次也是可以的。而最后一个REPIC则是今年刚刚表法的一个数据库。

图5 测序数据的m6A数据库

4、SE/DC和SC/DE的理解

通过7个实验17个组织或者细胞得到36个人样本这些样本来自MeT-DB v2.0这个数据集，样本被分成了两种变现形式：

1. same cell type across different experiments（SC/DE）
2. experiment but different cell types（ SE/DC）

5 、“本文中有关的参数”

1.p-value：是一种概率，在原假设为真的前提下，出现该样本或比该样本更极端的结果的概率之和。
参考文章：
https://blog.csdn.net/hahohehehe/article/details/79206002
https://www.zhihu.com/question/23149768
2. T检验用于小样本（样本容量小于30）的两个平均值差异程度的检验方法。它是用T分布理论来推断差异发生的概率，从而判定两个平均数的差异是否显著。
参考文章：
https://blog.csdn.net/hahohehehe/article/details/79206002
3.Cosine Similarity：余弦相似性通过测量两个向量的夹角的余弦值来度量它们之间的相似性。0度角的余弦值是1，而其他任何角度的余弦值都不大于1；并且其最小值是-1。从而两个向量之间的角度的余弦值确定两个向量是否大致指向相同的方向。
4.雅卡尔指数（Jaccard index）：又称为并交比（Intersection over Union）、雅卡尔相似系数（Jaccard similarity coefficient），是用于比较样本集的相似性与多样性的统计量。雅卡尔系数能够量度有限样本集合的相似度，其定义为两个集合交集大小与并集大小之间的比例，如下图所示：

一、论文创新点

作者论文的创新点主要有以下三点
1、m6Acorr在公共数据库中有效地纠正模拟数据集和真实m6A配置文件中的实验室偏差
2、m6Acorr能够更好地保存不同的甲基化信号
3、建立的服务器可以消除m6A甲基化谱中可能存在的实验室偏差，并根据修正的甲基化谱对高(低)甲基化基因进行基因谱比较和功能分析。

二、背景

基因甲基化水平的比较方式：
基因的相对甲基化水平可以描述为甲基化读计数(m6A-IP文库)与总读计数(输入文库)比较的富集分数。样本高甲基化(低甲基化)基因可以根据富集分数（enrichment score）很容易地确定。
甲基化图谱（methylation profiles）存在的缺陷：
1、MeT-DB数据的初步分析未能验证当前甲基化谱的相互可比性
2、SC/DE比 SE/DC的相对相似性更高，但是在人类数据库和老鼠的数据库中这个结论并不成立。具体的实验效果如下所示。
Cosine correlation的理解：余弦相关被用来证明同一组内的两个甲基化剖面之间的相似性，以避免由零分造成的伪影（artifacts）

图6 在人的基因数据集下的实验，低偏移的甲基化谱在SC/DE组应该比SE/DC组有显著更高的相关性，可以通过ComBat 方法和分位数规范化（quantile normalization），P值可以通过t-test方法得到

图7 在小鼠的数据集上，较少偏向的甲基化谱在SE/DC组应该比SC/DE组有显著更高的相关性

三、应用

（一）、模拟数据集

R (v3.6.1) 的库 Splatter (v1.10.0)被用于实验室偏置模拟RNA m6A甲基化谱splatEstimate 和 splatSimulate函数分别用于获得参数s和生成人工数据集实验室偏差。主要存在两个参数，通过这两个参数对网格搜索进行优化，以匹配来自真实数据集的m6A配置文件的分布最后，一个人工数据集，包含20个配置文件，包括4个批(实验)和两种细胞类型，以测试m6Acorr的性能。：

• batch.facLoc：表示位置，设置为0.3的效果最好
• batch.facScale：表示批量影响因子的对数-正态分布尺度，设置为0.2的效果最好
  

（二）、甲基化剖面校正方法

• Quantile normalization：
• Z-score 同一实验样本之间的归一化
• 每个实验变异系数的分割
• limma软件包中基于贝叶斯的经验批量回归方法
• 基于贝叶斯的经验批量回归方法
• RUVg方法
• 考虑轮廓之间可能存在的差异这些基因被调整到相同的分布与分位数归一化，以配合Combat model校正。

四、结果

（一）、本文方法的优点

• SE/DC在没有经过任何方法的处理下表现出比SC/DE更高的相似性，这与客观事实不符
• 通过Fig. 1b–d参考只有Z-score能够有效的纠正偏差
• Fig. 1 f ComBat 方法有效的纠正扭转了高的相似性
• Fig. 1h和Additional file 1: Figure S1在人和老鼠的数据集上可以看出将 ComBat和quantile normalization结合能够有效的消除SE/DC的高联系问题。
• ComBat和quantile normalization结合在人工数据集能够有效的解决SE/DC的高相关问题
• ComBat和quantile normalization组合的方法就叫做m6Acorr
  

图S3 模拟数据集中SE/DC组与SC/DE组的组内相关性比较， 通过quantile normalization and Combat处理的人工数据集能够有效的解决SE/DC的高相关问题

（三）、本文方法存在的问题

1.m6Acorr一个突出的问题是校正后生物信号的消除
解决方法：为了解决这个问题将腺苷酸在甲基化酶用p007_HeLa1_KO_M14和p007_HeLa2_KO_M3替换METTL3 and METTL14。这些图谱与上面对校正管道的评估无关，因为它们来自m6A酶突变细胞，而不是野生型细胞。通过计算前20%差异甲基化基因之间的Jaccard指数，比较校正前后共有的差异甲基化基因。

图S3 Jaccard指数描述了三种方法校正前后鉴定出的差异甲基化基因之间的共享分数。差异：前20%差异甲基化的基因；随机：随机选择相同数量的基因(重复100次，误差条显示标准误差

2.随机选择的基因不能获得共享的差异甲基化基因的比例。
解决方法：通过检查校正后鉴定的差异甲基化基因是否与功能性m6A靶基因表现出良好的一致性来检验这些方法。考虑了两类典型的功能性m6A靶基因：

• 第一类是其翻译效率受m6A修饰强烈调控的基因,根据GEO数据集GSE63591的记录，这些基因在METTL3或METTL14被敲除后，翻译效率显着降低。
• 第二类是其mRNA稳定性受m6A修饰强烈调控的基因,第二类是其mRNA稳定性受m6A修饰强烈调控的基因比较结果如图3所示，校正后鉴定的差异甲基化基因与任何一类功能性m6A靶基因都表现出良好的一致性，即使与从未经校正的甲基化图谱鉴定的差异甲基化基因相比也是如此。 这些结果表明，m6Acorr对甲基化图谱的校正也有助于寻找m6A调控的重要功能靶点。这些结果表明，m6Acorr对甲基化图谱的校正也有助于寻找m6A调控的重要功能靶点。
  
  

差异甲基化基因(校正前后)与功能性m6A靶基因的一致性。

这些结果表明，m6Acorr对甲基化图谱的校正也有助于寻找m6A调控的重要功能靶点。先前使用未经校正的甲基化图谱进行的比较分析表明，整个样本中的m6A甲基化广度与基因重要性相关的特征相关，包括dN/dS的比值、组织表达特异性、PPI网络。 这些显著相关在修正后被保留，如图所示：

图S4 M6A调控广度与各种基因重要性相关特征的相关曲线，给出了(A-C)校正前和(D-F)校正后的相关关系，拟合曲线采用黄土光滑法绘制。阴影表示置信区间。A.使用未经校正的剖面，将m6A调节宽度剖面与DN/DS比率相关。B.利用未经校正的剖面，将m6A调节宽度剖面与PPI网络进行关联。C.使用未校正的图谱，m6A调节宽度图谱与组织表达特异性的相关性。

四、讨论

• 这些方法对m6A甲基化图谱校正的适用性和有效性并不能自然得到保证。
• 甲基化水平是通过将甲基化读取计数与总读取计数进行比较从哪个样本来导出。
• 并不是所有众所周知的方法都能很好地用于甲基化剖面校正。尽管实验室的偏差非常的严重。
• 这项研究的新颖性集中在为什么以及如何将校正管道应用于m6A甲基化图谱
• m6Acorr服务器可以根据用户提供的批量(实验)任务执行甲基化配置谱校正。
• 如果没有分配批次，则整个数据集将被视为一个实验。此外，用户还可以分配试验组样本(例如，患病和健康)
• 可以根据m6Acorr服务器中经过整理的基因集注释来分析它们丰富的功能（见下图）
  

五、结论

• 本文针对公共数据库中rnam6A甲基化谱存在的实验室偏差，开发了基于分位数归一化和经验贝叶斯批量回归方法的m6Acorr。
• M6Acorr在真实数据集和人工数据集上都取得了良好的效果。而m6Acorr消除了偏差同时保存了生物关系。
• M6Acorr服务器还可以用来比较m6A图谱，并根据校正的甲基化图谱进行高(低)甲基化基因的功能分析。

本文数据集介绍