【原】科研 | 福建农林大学：对毛果杨的茎干分化木质部中的N6-甲基腺苷使用纳米孔RNA测序技术进行单碱基分辨率的定量分析（国人佳作）

编译：东方不赢，编辑：Tracy、江舜尧。

原创微文，欢迎转发转载。

导读

目前尚没有全面方法以单碱基分辨率鉴定单个转录物N⁶-甲基腺苷（m⁶A），但这对于估算m⁶A丰度又很重要，因此本文中研究者开发了一个名为Nanom6A的新方法，用于基于XGBoost模型的Nanopore RNA直接测序技术以单碱基分辨率鉴定和定量m⁶A的蛋白质修饰。研究者使用甲基化的RNA免疫沉淀测序（MeRIP-Seq）和m⁶A敏感的RNA-内切核糖核酸酶促测序（m⁶A-REF-seq）验证了该方法，证实了其准确性。研究提供了在茎干分化木质部中m⁶A修饰的转录组范围内的定量分析，并揭示了不同的替代聚腺苷酸化（APA）用法显示了不同的m⁶A比值。

论文ID

原名：Quantitative profiling of N⁶-methyladenosine at single-base resolution in stem-differentiating xylem of Populus trichocarpa using Nanopore direct RNA sequencing

译名：对毛果杨的茎干分化木质部中的N⁶-甲基腺苷使用纳米孔RNA测序技术进行单碱基分辨率的定量分析

期刊：Genomebiology

IF：10.806

发表时间：2021.01

通讯作者：顾连峰

通讯作者单位：福建农林大学

实验设计

实验结果

1.使用基于纳米孔DRS读取信号的XGBoost算法开发出的m⁶A预测模型

本文开发了一种新的测序流水线Nanom⁶A，基于XGBoost（极端梯度增强）模型，通过直接使用m⁶A位点周围的原始信号，以单核苷酸分辨率鉴定和可视化meA位点（图1）。研究作者分析了来自修饰和未修饰的RRACH k-mers的已公开Nanopore原始信号，发现修饰或未修饰位点的信号均值，中位数，标准差和信号宽度均不同。因此提取了所有读段的上述特征，并将它们按4：1的比例分为训练和测试数据集，以训练XGBoost模型（图1）。然后，研究者执行了10倍交叉验证，以评估模型的性能，其中，曲线下面积（AUC）为97％。之后将XGboost与其他几种分类算法进行了比较，结果发现XGBoost表现出最好的性能，因此，在后续研究中使用了XGBoost模型。                           

图1 Nanom⁶A方法用于以单碱基分辨率鉴定m⁶A

对于XGBoost训练模块，按需修改m⁶A的原始序列。使用Tombo将校正后的序列分配给原始信号后，研究者根据RRACH特征（包括信号长度，平均值，标准和宽度）提取了不同碱基的信号。之后用软件训练了相应的12个k-mers，以分别构建XGBoost模型。对于预测模块，将碱基文件与基因组序列进行比对。然后用基因组序列校正该序列。将校正后的序列分配给原始信号后，提取各个碱基的信号与对应的RRACH k-mer的相应特征，来使用训练的模型来预测修饰碱基。 

为了验证Nanom⁶A是否可以从单个转录本中识别出m⁶A修饰，研究者使用已知的m⁶A与非m⁶A比率的合成mRNA已知样品评估了此方法。首先，研究者使用Nanom⁶A模型分别预测已知的修饰和未修饰读段，总共可以成功鉴定出单个转录本中91-96％的修饰和未修饰位点；其次，研究者将已知的修饰和未修饰读段混合在一起，以模拟m⁶A与非m⁶A比率不同的数据集，涵盖了从0（所有未修饰的DRS读段）到1（所有修饰的读段）修饰比率的整个区间。研究者使用Nanom⁶A分析了DRS模拟数据集，发现基于Nanom⁶A预测的m⁶A与非m⁶A比率，与已知的m⁶A与非m⁶A比率具有很强的相关性。研究者还模拟了不同的序列深度，通过将各10000次对应于不同m⁶A与非m⁶A比率的20、50和100个DRS读数进行混合来生成三种不同的水平。这些结果表明，较高的深度可以生成较高的相关性。

2.在哺乳动物和植物中使用已知的m⁶A位点验证了Nanom⁶A

在毛果杨（P. trichocarpa）中使用此方法之前，研究者首先使用哺乳动物和植物中已发布的DRS数据来进一步验证Nanom⁶A的可靠性（图2a）。首先，研究者使用来自野生型和METTL3敲除突变体的原始信号文件（shMETTL3），通过Nanom⁶A识别m⁶A位点，并使用连接辅助提取的方法和薄层色谱法与ACTB位点（1216）中的已知位点进行比较（猩红色）。研究者使用Tombo的重新波形算法提取了信号分配后的平均值，标准差，中位数和停留时间（图2a）。对于两个从野生型和shMETTL3读取的DRS，研究者在ACTB位点（1216）的修改概率分别为0.974和0.003（图2a），结果与使用SCARLET报道的ACTB修饰位点（1216）一致。通过使用相同的策略，研究者为其他DRS读取计算了ACTB站点（1216）的修改概率，根据来自ACTB的m⁶A读数/总DRS读数的公式，ACTB站点的m⁶A比值（1216）分别为0.62和0.33。研究者还使用这种方法来计算所有其他基因的比率，这些比率已被DRS读取所覆盖，结果表明，在METTL3敲低样品中，修饰率（图2b）和m⁶A位点数量（图2c）均降低。

图2 使用哺乳动物中已知的m⁶A位点验证Nanom⁶A

a.分别使用野生型和shMETTL3的两个DRS读数中的离子电流信号从ACTB位点（1216）识别m⁶A位点的策略。b.条形图分别显示了来自野生型和shMETTL3的预测m⁶A比率。c. Venn图分别为野生型和shMETTL3预测的m⁶A位点数量。d.条形图显示了基于SCARLET方法和Nanom⁶A预测的已知m⁶A位点之间的一致性。x轴表示TUG1中的修饰位点。y轴显示已修改读段的数量。e.三个柱状图分别显示了ACTB和Nanom⁶A预测中已知修饰位点之间的一致性，输入读取中MeRIP-Seq IP的富集以及Nanom⁶A修饰的比例。f.条形图分别显示了BSG和Nanom⁶A预测中已知修饰位点之间的一致性，输入读取时MeRIP-Seq IP的富集以及Nanom⁶A修饰比例之间的一致性

为了进一步评估这个方法的可靠性，研究者基于SCARLET计算了所有已知的m⁶A位点。总的来说，科研者们已经报道了两个人类lncRNA（MALAT1和TUG1）和三个人类mRNA（TPT1，ACTB和BSG）的m⁶A修饰，这些结果代表了从体内的转录本中可以得到可靠的已知修饰位点。研究者将基于DRS数据，通过Nanom⁶A识别的m⁶A位置与基于SCARLET方法的已知m⁶A修饰进行了比较。在两个lncRNA（MALAT1和TUG1）中，仅三个DRS读数支持MALAT1（NR_002819.2），而这些读段太少，无法覆盖已知的m⁶A位点。因此，研究者从下游分析中排除了MALAT1的比较。在TUG1 RNA中，Nanom⁶A鉴定了十个已知的修饰位点，与SCARLET方法一致（图2d）。特别之处在于，研究者发现位于3'末端区域附近的这些位点包含更多的DRS读数支持。对于TPT1 mRNA（NM_003295.1），基于Nanom⁶A的DRS读数在687、694和703个修饰位点也与SCARLET呈现一致的结果。基于SCARLET方法，ACTB mRNA（NM_001101.5）和BSG mRNA（NM_198591.1）具有三个已知的m⁶A修饰，该方法支持这些修饰位点并与SCARLET结果一致（图2e，f）。值得注意的是，基于Nanom⁶A，DRS数据显示了在ACTB的1216位点和BSG的1335位点为两个高修饰率位点的。因此，研究者选择了这些m⁶A修饰位点，以进行野生型和METTL3组合（shMETTL3）的定性比较，发现使用Nanom⁶A，野生型和shMETTL3中ACTB在1216位的修饰率分别为61.6％和33.1％（图2e）；在野生型和shMETTL3中，BSG转录物在1335位点的修饰率分别为60％和47％（图2f）。已发布的MeRIP-Seq数据集也报告了类似的下降趋势，这进一步验证了基于Nanom⁶A的结果。

图3 使用植物中已知的m⁶A位点验证Nanom⁶A

a.Nanom⁶A，miCLIP和EpiNano中所有m⁶A位点的比较。b.在Nanom⁶A，MINES和EpiNano中比较四个基序（AGACT，GGACA，GGACC和GGACT位点）。c.比较先前研究报告的差异错误位点（DES）和基于vir-1突变体和VIR互补系的Nanom⁶A鉴定的m⁶A的m⁶A位点。d.分别来自vir-1，VIRc和Col-0的m⁶A位点的分布。e.散点图显示vir-1和VIRc之间的m⁶A比。f.箱形图显示了在vir-1，VIRc和Col-0中的m⁶A比率。g.CCA1的条形图分别显示了vir-1，VIRc和Col-0的m⁶A修饰率。

研究者进一步将本方法与拟南芥DRS数据的几种公开方法进行了比较，发现用Nanom⁶A方法鉴定的3770个m⁶A位点与使用EpiNano或miCLIP检测到的修饰重叠（图3a）。尤其是Nanom⁶A报告了用miCLIP检测到1667个m⁶A修饰位点，而EpiNano未检测到这些位点（图3a）。对于AGACT，GGACA，GGACC和GGACT基序，研究者用Nanom⁶A方法检测到的大约40％m⁶A位点与MINES或EpiNano重叠（图3b），而进一步比较前人文献中17,491个微分错误位点（DES），发现其中包括总共4936个具有RRACH基序的位点，约有66％（3289/6936）来自文献的RRACH基序，Nanom⁶A也结果也类似（图3c）。与VIR互补系相比，使用DRS数据基于Nanom⁶A的m⁶A位点数量表明VIRILIZER（vir-1）突变体中m⁶A位点明显减少（图3c）。有趣的是，研究者发现VIR的消耗减少了终止密码子和3'UTR附近m⁶A位点的富集（图3d）。该结果与先前报道的VIRMA功能一致，该功能介导终止密码子和3'UTR中的甲基化。

在vir-1中，基于Nanom⁶A预测的m⁶A位点修饰比例（图3e，f）也降低了，这与VIR（保守的m⁶A编写器复杂组件）的功能一致，进一步验证了研究者方法的可靠性。有趣的是，昼夜节律的调节剂CCA1在拟南芥中的mRNA中具有m⁶A修饰；Nanom⁶A也鉴定出该修饰位点，在vir-1突变体中m⁶A修饰位点的比例降低（图3g）。

3. 毛果杨茎干分化木质部中m⁶A的定性分析

使用在腺嘌呤甲基化，去甲基化和m⁶A读数中起作用的MTA70，WTAP和YTH蛋白质的保守结构域，研究者鉴定了所有主要成分（7个m⁶A甲基转移酶，15个m⁶A脱甲基酶和18个m⁶A结合蛋白），说明了P. trichocarpa保守的m⁶A修饰机制。接下来研究者从茎分化木质部（SDX）中提取了poly（A）+ RNA，并使用GridION和MinION平台进行了直接RNA测序（图4a）。具体如下：首先，在使用训练好的模型之前，将信号正确分配给相应的碱基（图1）；然后通过Guppy（V3.6.1）进行碱基检出后，使用minimap2（V3）将纳米孔长读段映射到参考序列，Tombo（v1.5）的重采样模块用于将原始信号分配给每个单独的基准。最后使用经过训练的XGBoost模型（图1）提取每个单独序列的原始信号以进行m⁶A预测。分析可知，Nanopore DRS读段显示了共计22,953个基因的表达，其中，从repeat1和repeat2分别识别出总共12,338和29,380个唯一的m⁶A位点，每个转录本的平均m⁶A位点数为4。m⁶A的分布表明修饰位于终止密码子和3'UTR区域附近（图4b），进一步验证了研究者方法的可靠性。

图4 毛果杨的茎分化木质部中m⁶A的定量分析

a.DRS库构建和测序流程图。b.转录本中m⁶A位点的分布。c.比较两个独立的生物学重复的m⁶A位点。d.方框图显示了DRS读取重叠的m⁶A位点和重复2特异性m⁶A位点的覆盖范围。e.纳米孔DRS支持的所有RRACH基序甲基化均读取。f.所有修饰基因和与木材形成相关的基因的m⁶A比。g.m⁶A比率与基因表达之间的相关性。h.Venn图显示了用Nanom⁶A和Epinano预测的m⁶A位点之间的重叠。i.m⁶A位点在重叠位点，Nanom⁶A特定位点和EpiNano特定位点的分布。j.重叠位点和Nanom⁶A特定位点的m⁶A比。k.来自DRS，EpiNano，MINES和MeRIP-Seq的ARK2的修饰位点。MeRIP-Seq中的蓝色和红色分别表示输入library和IPlibrary。

 

研究者通过使用不同的序列平台（GridION和MinION）和不同的试剂盒（SQK-RNA001和SQK-RNA002）进行两次生物学重复， 1,087个重叠的m⁶A位点的Pearson相关系数为0.96，这表明Nanom⁶A即使使用不同的文库制备试剂盒和不同的测序仪，也可以提供良好的可重复性。来自第一个重复序列的总共98％m⁶A位点与第二个重复序列重叠，表明两个独立的生物学重复序列之间具有很高的可重复性（图4c）。由于较高的测序深度，第二个重复序列鉴定出17,293个新的m⁶A位点，这些位点在第一个重复序列中缺失。显然，这些repeat2特异的修饰位点主要源自低丰度的转录本（图4d），这表明提高测序深度可以增加鉴定m⁶A位点的机会。SDX中m⁶A位点的平均比例为0.44，每个基序支持的DRS读数的分布如图4e所示，而含有m⁶A的外显子的长度比对照的长，与人类一致。

在与木材形成相关的76个基因中，分别来自17、26和11个半纤维素，木质素和纤维素基因包含m⁶A的转录本，这表明富含m⁶A修饰的转录本编码的蛋白质参与木材形成。但是，研究者发现与木材形成相关的基因的修饰率（~0.20）低于平均值的修饰率（图4f）。在调查了最低m⁶A水平的前1000个基因后，研究者发现涉及木质素生物合成过程（P=5.54E⁻¹⁵），木聚糖生物合成过程（P=2.78E⁻⁰⁶）和次生细胞壁生物发生（P=5.08E^-17）相关基因被富集。此外，研究者发现基因表达水平与m⁶A水平呈负相关（图4g）。因此，SDX中与木材形成相关的基因中较低的m⁶A水平可能是有助于这些转录物的稳定性及其高表达水平的因素之一。在哺乳动物系统中，已证明m⁶A修饰会增加mRNA的衰变，需要进一步的研究来探究m⁶A修饰是否调节与木材形成相关的基因的表达水平。

最后，研究者提供了使用EpiNano（图4h-j）和MINES使用相同cutoff（>20 DRS读取支持）的方法的比较。Nanom⁶A与EpiNano之间的维恩图描绘了共计2326个重复m⁶A位点，10012个Nanom⁶A特异位点和2838个EpiNano特异位点的重叠（图4h）。m⁶A位点的分布表明，与2838 EpiNano特异性位点相比，2326个重叠且10,012个Nanom⁶A特异的m⁶A位点在终止密码子和3'UTR区域附近更富集（图4i），尤其是研究者发现来自10,012个Nanom⁶A特异性位点的m⁶A比率低于2326个重叠位点的比率（图4j），这表明Nanom⁶A在检测低m⁶A / A比率方面具有优势，因为Nanom⁶A可以在单转录本水平检测到m⁶A。MINES仅为AGACT，GGACA，GGACC和GGACT生成了m⁶A位点。因此，研究者仅使用这三种方法与MINES比较了上述四个基序。Nanom⁶A与MINES之间的比较也显示了Nanom⁶A与EpiNano相似的趋势。纤维素ARK2包含三个基于Nanom⁶A的m⁶A位点。通过所有方法可确定远端m⁶A部位（图4k）。

4. 使用MeRIP-Seq和m⁶A-REF-seq验证Nanom⁶A

为了进一步验证基于纳米孔DRS的预测m⁶A位点，研究者还用抗m⁶A抗体进行了MeRIP-Seq（图5a）。使用KAPA链式mRNA-Seq试剂盒对甲基化的RNA和输入的RNA进行cDNA文库构建，并在IlluminaNovaseq平台上测序。Hisat2（v2.0.3）用于将MeRIP-Seq读段与基因组比对，独特的定位读段用PEA R软件包（v1.1）来调用的m⁶A峰，来自MeRIP-Seq的免疫沉淀读数也显示了终止密码子和3'UTR区附近的富集（图5b）。来自两个重复序列的重叠m⁶A位点，比来自第二个重复序列的新颖m⁶A位点具有更高的覆盖率，这表明MeRIP-Seq的更高测序深度可以从低丰度转录物中检测到更多的m⁶A位点（图5c），而研究者通过MeRIP-Seq总共检测到81％（2626/3253）由DRS预测的m⁶A修饰基因。在单碱基分辨率下，MeoreIP-Seq峰覆盖了Nanopore DRS预测的m⁶A位点的49％（图5d）。

图5 使用MeRIP-Seq和m⁶A-REF-seq / MAZTER-seq验证Nanom⁶A

a.概述MeRIP-Seq的流程图。b.MeRIP-Seq中的m⁶A分布跨转录本读段。c.箱形图显示了重叠的m⁶A位点和repeat2特异位点的读数覆盖率。d.通过MeRIP-Seq验证的m⁶A位点的百分比。e. m⁶A-REF-seq / MAZTER-seq的流程图。f.条形图显示了序列读取结束时ACA基序的百分比。g.跨转录本的ACA修饰位点的分布。h.通过m⁶A-REF-seq/ MAZTER-seq验证的m⁶A位点的百分比。

m⁶A-REF-seq或MAZTER-seq ACA基序中的m⁶A位点检测提供了单碱基分辨率，为了进一步验证Nanom⁶A的可靠性，研究者还进行了m⁶A-REF-seq验证基于直接RNA测序的方法。研究者首先按照先前的m⁶A-REF-seq文库构建方法在毛果杨中进行m⁶A-REF-seq；然后，使用IlluminaNova6000平台对m⁶A-REF-seq库进行测序（图5e）。研究者发现，在测序读段结束时，ACA基序的富集，以及来自Illumina末端修复模型的多于60％左右读段分别在ACA位点分别开始和终止（图5f）。研究者在测序读取末尾的ACA百分比结果与先前一致，表明m⁶A-REF-seq具有很高的可靠性，而从m⁶A-REF-seq结果得到的ACA基序中的m⁶A位点在终止密码子附近富集（图5g）。总体而言，m⁶A-REF-seq验证了基于Nanom⁶A80％的m⁶A位点（图5h）。

5.Populus trichocarpa分化木质部的poly（A）尾长度分析

纳米孔DRS在检测poly（A）尾巴长度方面具有巨大优势，这在拟南芥体外转录的RNA 和GM12878细胞中已有报道。在这项研究中，研究者采用了用于天然RNA测序的标准文库制备方案，该方案保留了完整的poly（A）尾巴，以鉴定同工型特异性poly（A）尾巴长度。研究者在nanopolish 安装包（0.11.1）中用nanopolish-polya模块确定了poly（A）的尾长，他们发现P. trichocarpa mRNA的平均和中位数poly（A）长度分别为92 nt和82 nt（图6a）。基因表达与poly（A）尾巴长度在Populus中呈负相关（r=−0.26）（图6b）。总共有2421个备选的聚腺苷酸化基因，包括近端和远端poly（A）位点，这些位点显示poly（A）长度发生了两倍变化（图6c）。研究者的DRS数据表明，该次生壁生物合成基因包括两个具有不同poly（A）长度的聚腺苷酸化位点（图6d），具有远端聚腺苷酸化的转录物比近端转录物具有更长的poly（A）。这项分析提供了替代聚腺苷酸的初步poly（A）尾长，以表明poly（A）尾在降解和翻译中的潜在调控作用。 

图6 poly（A）尾的不同长度和不同的聚腺苷酸化位点的不同m6A修饰

a.基于nanopolish-polyA的poly（A）尾部长度的密度图。b.散点图，显示了基因表达与poly（A）尾巴长度之间的相关性。c.散点图显示了poly（A）尾巴长度，DRS从近端和远端转录本丰度。d. PARVUS-L-2的基因结构和DRS读数，其具有两个带有不同poly（A）尾长的聚腺苷酸化位点。e.上图中的密度图分别显示了远端（D）和近端（P）polyA部位甲基化率的倍数变化。下部面板中的图分别显示了远端和近端poly（A）位置的m⁶A比率。f.分别具有远端或近端poly（A）位点的转录本的不同m⁶A修饰的代表性示例。

6.毛果杨茎分化木质部中m⁶A的定量分析

本文方法的优点是研究者可以从原始的纳米孔DRS数据中量化m⁶A，因此，研究者的方法根据修饰的和未修饰的转录本计算每个m⁶A位点的比率，可以在单碱基分辨率下检测每个转录物的m⁶A位点，纳米孔DRS可以识别远端和近端poly（A）位置，从DRS分析得到的Poly（A）位点相邻24 nt之内会被分为一个poly（A）群。先前的研究表明，受损的m⁶A甲基化酶复合物可以改变poly（A）位点的使用。在所有经过m⁶A修饰的基因中，研究者发现使用3152个替代性聚腺苷酸化（APA）的基因，而且，这项研究中的方法可以区分m⁶A和非m⁶A转录本。基于甲基化/非甲基化DRS转录本的甲基化比率的分布显示，转录本具有远端或近端poly（A），修饰百分比不同（图6e），例如，来自Potri.019G083300的具有远端poly（A）的转录本比邻近的poly（A）转录本具有更少的m⁶A位点（图6f）。这项研究提供了初步数据，以进一步研究m⁶A修饰的同工型与替代的聚腺苷酸化。

结论

最近，纳米孔直接RNA测序被用于检测RNA中的碱基修饰信号，这已由EpiNano和MINES 研究得到，然而，EpiNano无法将m6A与其他类型的RNA修饰区分开，例如m1A，因为该方法基于碱基质量和缺失频率来预测m6A位点。

MINES仅在某些序列上下文（AGACT，GGACA，GGACC和GGACT）中检测到m6A位点。在这项研究中，研究者开发了Nanom6A，这是一种用于以单碱基分辨率鉴定m6A修饰的新方法，以克服使用基于XGBoost模型的Nanopore直接RNA读数的现有局限性，与以前使用SVM 或随机森林的方法不同。重要的是，这与已发布的DRS数据方法不同，可以通过研究者的Nanom⁶A方法对m6A位点的数量进行量化，因为该方法可以在转录本分辨率下提供m⁶A修饰。为了验证基于Nanopore RNA直接测序的方法，研究者使用了MeRIP-Seq和m6A-REF-seq，这表明Nanom⁶A可以在检测m6A的定性和定量分析中实现高精度。使用此方法，研究者以基于转录物的分辨率提供了SDX中m⁶A修饰的转录组范围的鉴定和定量，并揭示了不同的APA用量显示了不同的m⁶A比例。研究者的方法大大扩展了纳米孔直接RNA测序在探索m⁶A调控机制中的应用。

原文链接：  

https://pubmed.ncbi.nlm./33413586/