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

编译:东方不赢,编辑:Tracy、江舜尧。

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

导读

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

论文ID

原名:Quantitative profiling of N6-methyladenosine at single-base resolution in stem-differentiating xylem of Populus trichocarpa using Nanopore direct RNA sequencing
译名:对毛果杨的茎干分化木质部中的N6-甲基腺苷使用纳米孔RNA测序技术进行单碱基分辨率的定量分析
期刊:Genomebiology
IF:10.806
发表时间:2021.01
通讯作者:顾连峰
通讯作者单位:福建农林大学

实验设计

实验结果

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

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

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

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

2.在哺乳动物和植物中使用已知的m6A位点验证了Nanom6A

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

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

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

图3 使用植物中已知的m6A位点验证Nanom6A
a.Nanom6A,miCLIP和EpiNano中所有m6A位点的比较。b.在Nanom6A,MINES和EpiNano中比较四个基序(AGACT,GGACA,GGACC和GGACT位点)。c.比较先前研究报告的差异错误位点(DES)和基于vir-1突变体和VIR互补系的Nanom6A鉴定的m6A的m6A位点。d.分别来自vir-1,VIRc和Col-0的m6A位点的分布。e.散点图显示vir-1和VIRc之间的m6A比。f.箱形图显示了在vir-1,VIRc和Col-0中的m6A比率。g.CCA1的条形图分别显示了vir-1,VIRc和Col-0的m6A修饰率。

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

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

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

使用在腺嘌呤甲基化,去甲基化和m6A读数中起作用的MTA70,WTAP和YTH蛋白质的保守结构域,研究者鉴定了所有主要成分(7个m6A甲基转移酶,15个m6A脱甲基酶和18个m6A结合蛋白),说明了P. trichocarpa保守的m6A修饰机制。接下来研究者从茎分化木质部(SDX)中提取了poly(A)+ RNA,并使用GridION和MinION平台进行了直接RNA测序(图4a)。具体如下:首先,在使用训练好的模型之前,将信号正确分配给相应的碱基(图1);然后通过Guppy(V3.6.1)进行碱基检出后,使用minimap2(V3)将纳米孔长读段映射到参考序列,Tombo(v1.5)的重采样模块用于将原始信号分配给每个单独的基准。最后使用经过训练的XGBoost模型(图1)提取每个单独序列的原始信号以进行m6A预测。分析可知,Nanopore DRS读段显示了共计22,953个基因的表达,其中,从repeat1和repeat2分别识别出总共12,338和29,380个唯一的m6A位点,每个转录本的平均m6A位点数为4。m6A的分布表明修饰位于终止密码子和3'UTR区域附近(图4b),进一步验证了研究者方法的可靠性。
图4 毛果杨的茎分化木质部中m6A的定量分析
a.DRS库构建和测序流程图。b.转录本中m6A位点的分布。c.比较两个独立的生物学重复的m6A位点。d.方框图显示了DRS读取重叠的m6A位点和重复2特异性m6A位点的覆盖范围。e.纳米孔DRS支持的所有RRACH基序甲基化均读取。f.所有修饰基因和与木材形成相关的基因的m6A比。g.m6A比率与基因表达之间的相关性。h.Venn图显示了用Nanom6A和Epinano预测的m6A位点之间的重叠。i.m6A位点在重叠位点,Nanom6A特定位点和EpiNano特定位点的分布。j.重叠位点和Nanom6A特定位点的m6A比。k.来自DRS,EpiNano,MINES和MeRIP-Seq的ARK2的修饰位点。MeRIP-Seq中的蓝色和红色分别表示输入library和IPlibrary。
研究者通过使用不同的序列平台(GridION和MinION)和不同的试剂盒(SQK-RNA001和SQK-RNA002)进行两次生物学重复, 1,087个重叠的m6A位点的Pearson相关系数为0.96,这表明Nanom6A即使使用不同的文库制备试剂盒和不同的测序仪,也可以提供良好的可重复性。来自第一个重复序列的总共98%m6A位点与第二个重复序列重叠,表明两个独立的生物学重复序列之间具有很高的可重复性(图4c)。由于较高的测序深度,第二个重复序列鉴定出17,293个新的m6A位点,这些位点在第一个重复序列中缺失。显然,这些repeat2特异的修饰位点主要源自低丰度的转录本(图4d),这表明提高测序深度可以增加鉴定m6A位点的机会。SDX中m6A位点的平均比例为0.44,每个基序支持的DRS读数的分布如图4e所示,而含有m6A的外显子的长度比对照的长,与人类一致。
在与木材形成相关的76个基因中,分别来自17、26和11个半纤维素,木质素和纤维素基因包含m6A的转录本,这表明富含m6A修饰的转录本编码的蛋白质参与木材形成。但是,研究者发现与木材形成相关的基因的修饰率(~0.20)低于平均值的修饰率(图4f)。在调查了最低m6A水平的前1000个基因后,研究者发现涉及木质素生物合成过程(P=5.54E15),木聚糖生物合成过程(P=2.78E06)和次生细胞壁生物发生(P=5.08E-17)相关基因被富集。此外,研究者发现基因表达水平与m6A水平呈负相关(图4g)。因此,SDX中与木材形成相关的基因中较低的m6A水平可能是有助于这些转录物的稳定性及其高表达水平的因素之一。在哺乳动物系统中,已证明m6A修饰会增加mRNA的衰变,需要进一步的研究来探究m6A修饰是否调节与木材形成相关的基因的表达水平。
最后,研究者提供了使用EpiNano(图4h-j)和MINES使用相同cutoff(>20 DRS读取支持)的方法的比较。Nanom6A与EpiNano之间的维恩图描绘了共计2326个重复m6A位点,10012个Nanom6A特异位点和2838个EpiNano特异位点的重叠(图4h)。m6A位点的分布表明,与2838 EpiNano特异性位点相比,2326个重叠且10,012个Nanom6A特异的m6A位点在终止密码子和3'UTR区域附近更富集(图4i),尤其是研究者发现来自10,012个Nanom6A特异性位点的m6A比率低于2326个重叠位点的比率(图4j),这表明Nanom6A在检测低m6A / A比率方面具有优势,因为Nanom6A可以在单转录本水平检测到m6A。MINES仅为AGACT,GGACA,GGACC和GGACT生成了m6A位点。因此,研究者仅使用这三种方法与MINES比较了上述四个基序。Nanom6A与MINES之间的比较也显示了Nanom6A与EpiNano相似的趋势。纤维素ARK2包含三个基于Nanom6A的m6A位点。通过所有方法可确定远端m6A部位(图4k)。
4. 使用MeRIP-Seq和m6A-REF-seq验证Nanom6A

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

图5 使用MeRIP-Seq和m6A-REF-seq / MAZTER-seq验证Nanom6A

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

m6A-REF-seq或MAZTER-seq ACA基序中的m6A位点检测提供了单碱基分辨率,为了进一步验证Nanom6A的可靠性,研究者还进行了m6A-REF-seq验证基于直接RNA测序的方法。研究者首先按照先前的m6A-REF-seq文库构建方法在毛果杨中进行m6A-REF-seq;然后,使用IlluminaNova6000平台对m6A-REF-seq库进行测序(图5e)。研究者发现,在测序读段结束时,ACA基序的富集,以及来自Illumina末端修复模型的多于60%左右读段分别在ACA位点分别开始和终止(图5f)。研究者在测序读取末尾的ACA百分比结果与先前一致,表明m6A-REF-seq具有很高的可靠性,而从m6A-REF-seq结果得到的ACA基序中的m6A位点在终止密码子附近富集(图5g)。总体而言,m6A-REF-seq验证了基于Nanom6A80%的m6A位点(图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)位置的m6A比率。f.分别具有远端或近端poly(A)位点的转录本的不同m6A修饰的代表性示例。
6.毛果杨茎分化木质部中m6A的定量分析

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

结论

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

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

原文链接:
https://pubmed.ncbi.nlm.nih.gov/33413586/
(0)

相关推荐