这篇6+文章告诉你m6A生信分析还能这么做！

导语

GUIDE ╲

卵巢癌(OC)是所有妇科肿瘤中导致死亡的主要原因。N6-甲基腺苷(m6A)相关调控因子在包括OC在内的多种肿瘤中发挥重要作用。然而，m6A RNA甲基化调控因子的表达及其与预后的相关性仍在很大程度上尚不清楚。

背景介绍

可能是由于近两年m6A的火热，20多个m6A相关调控蛋白已经被进行了多方面的分析，近期m6A生信分析的新鲜思路比较少有，今天小编给大家带来一篇很有意思的研究思路，作者基于几个m6A调节蛋白构建风险评分模型，并引入miRNA-m6A调控因子-m6A靶基因调控网络预测OC患者预后的文章。这篇文章于2021年发表在《frontiers in Cell and Developmental Biology》杂志上，最新影响因子6.684，题目为：A Risk Score Model IncorporatingThree m6A RNA MethylationRegulators and a Related Network ofmiRNAs-m6A Regulators-m6A TargetGenes to Predict the Prognosis ofPatients With Ovarian Cancer。

图1

数据介绍

GDC:374例OC患者RNA-seq转录组数据

GTEx:188例人正常卵巢组织RNA-seq数据

从最近发表的研究中收集21个关键的m6A甲基化调控因子，包括11个reader、7个writer和3个earser。

结果解析

01

OC组织中21个m6A甲基化调控因子的表达

作者首先利用GDC和GTEx数据库的数据集，比较了OC和正常人类卵巢组织中21个关键的m6A甲基化调控因子的mRNA表达水平。与正常对照组相比，21个m6A调控因子中有19个在OC组织中存在差异表达（图2A）。

根据OC患者的临床特征将其分为low-stage组(I+II期)和high-stage组(III+IV期)，发现21个m6A调节因子中有9个(METTL14、YTHDC2、FTO、ALKBH5、HNRNPA2BA2B1、VIRMA、IGF2BP1、RBM15、RBMX)在两组中表达差异显著，high-stage组与low-stage组相比均表达下调(图2B)。

图2

02

开发由三个m6A调节因子组成的风险signature

采用LASSO Cox回归算法分析GDC数据集中的21个调控因子，更好地预测m6A甲基化调控因子异常表达患者的临床结果，获得了基于三个m6A调节因子(VIRMA、IGF2BP1和HNRNPA2B1)作为变量的风险评分(图3A，B)。以中位风险评分（-0.03)作为cut-off，将所有患者(train）分为两组，即高风险组和低风险组。还构建了一个由GDC数据库中随机选择的252名OC患者组成的测试数据集。随后进行KM生存分析，以评估三个选定的m6A调节因子的有效性，并构建训练和测试数据集中OC患者生存率的风险模型(图3C，D)。

图3

生存率为50%、VIRMA低表达或HNRNPA2B1高表达的OC患者的生存时间长于对照组(图4A，B)。然而，IGF2BP1的表达水平与OC患者的5年生存率无关(图4C)。然而，当在train数据集中使用基于三个m6A调节器的风险模型时，低风险组生存率为50%的OC患者的生存时间明显高于高危组患者(p=4e-05)，表示其潜在的临床意义(图4D)。

对所有变量应用Cox比例风险模型的单变量分析表明，风险模型、分期和肿瘤残留疾病均是预测OC患者不良预后的独立因素(图4E)。对于测试数据集，低风险组生存率为50%的OC患者的生存时间也明显高于高危组患者(p=0.0058，图4F)。此外，对OC患者应用Cox比例风险模型的单因素分析也表明，风险模型(P=0.026)也是预测OC患者不良预后的独立因素(图4G)，这与train数据集的结果一致。

图4

03

基于临床病理特征建立列线图

作者建立了一个列线图，以帮助了解生存时间与临床病理特征之间的关系。通过计算列线图对应的total points，能够估计OC患者在3年和5年的生存概率(图5A)。列线图对OC患者的3年和5年生存率具有良好的预测能力(图5B)。曲线下面积(AUC)计算为0.6,3年面积为0.62,5年计算为0.64(图5C)，提示预后风险评分具有中等的准确性。

图5

04

预测靶向三个选定的m6ARNA甲基化调节因子的miRNAs

从GDC数据集中获得358个与三个选定的m6A调控因子的表达呈负相关的miRNA。基于miRbase数据库的生物信息学分析，预测了40个具有HNRNPA2B1潜在结合位点的miRNAs，46个具有VIRMA结合位点，68个具有IGF2BP1结合位点的miRNAs(共133个非重复miRNAs)。将上一步分析结果取交集，只剩下两个miRNAs，即hsa-miR-196b-5p和hsa-miR-98-5p，它们在IGF2BP1中具有潜在的结合位点，与OC患者的预后显著相关(图6A)。

通过进一步的表达和生存相关分析发现与low-stage组相比，high-stage组hsa-miR-miR-196b-5p表达增加，而hsa-miR-196b-5p表达的hsa-miR-196b-5p表达的OC患者生存率明显降低(p=0.0016)，表明其在OC中有致癌功能(图6B、C)。同时，与low-stage组相比，high-stage组hsa-miR-98-5p表达量降低，hsa-miR-miR-98-5p表达量高时，hsa-miR-98-5p表达量高的OC患者生存率提高(p=0.044)，提示其在OC中具有潜在的抑癌功能(图6B、D)。

图6

05

miRNA-m6A调控基因-m6A靶基因网络的构建

380个m6A甲基化相关基因来自m6Avar数据库，并进一步根据person相关系数(p<0.05)与3个m6A调控因子(VIRMA、IGF2BP1和HNRNPA2B1)进行过滤，构建miRNA、m6A调控因子和m6A靶基因之间的调控网络。获得123个VIRMA相关基因，110个IGF2BP1相关基因，148个HNRNPA2B1相关基因，通过GO和KEGG通路富集分析，发现其中47个可能共同调控三个m6A调控因子的基因(图7A)。GO分析显示，这48个基因的生物过程主要富集于调控蛋白质转运和定位等生物过程(图7B，C)。KEGG通路分析表明，这些基因主要富集于病毒癌变和MAPK信号通路中，这可能表明它们与OC的发生进展有关(图7C)。

然后，构建了一个miRNAs-m6A调控因子-m6A靶基因的调控网络，该网络由两个miRNA、3个m6A调控因子和47个mRNA组成(图7D)。还构建了两个选定的miRNAs(miR-196b-5p和miR-98-5p)、IGF2BP1和110个预测的m6A靶基因或10个先前报道的m6A靶基因的网络(图7E，F)。

图7

06

选择一个miRNA-m6A调因子-m6A靶基因通路初步验证

为了确定生物信息学分析的结果的可靠性，作者首先检测了三个选定的m6A调控因子(VIRMA、IGF2BP1和HNRNPA2B1)在OC组织中的表达。PCR和IHC检测结果显示，与low-stage OC组织相比，high-stage OC组织中VIRMA和IGF2BP1表达下调，而与low-stage OC组织相比，HNRNPA2B1中HNRNPA2B1表达上调(图8A-D)。此外还发现，根据Human Protein Atlas数据库，所选择的三种m6A调控因子(VIRMA、IGF2BP1和HNRNPA2B1)在OC组织中的一些IHC数据与作者的结果一致(图8E)。

为了进一步验证构建网络的可信度，作者选择了 miRNA-m6A 调节因子-m6A 相关基因通路之一，即 miR-196b-5p-IGF2BP1-PTEN 进行初步验证。 由于miR-196b-5p的功能及其与OC中IGF2BP1的关系尚不清楚，作者首先检测了OC细胞和正常卵巢上皮细胞中miR-196b-5p和IGF2BP1 mRNA的表达水平，发现miR-196b 与对照细胞相比，卵巢癌细胞中的miR-196b-5p上调且 IGF2BP1 下调（图 8F、G）。

还检测了之前描述的miR-196b-5p在30个OC组织中的表达，发现与IGF2BP1负相关的miR-196b-5p(r=-0.583，P<0.05)，在高水平组的OC组织中比对照组上调(图8H，I)。最后基于一系列生物学实验验证了miR-196b-5p-IGF2BP1-PTEN在OC中的调控通路，从而证明构建的调控网络具有一定的可靠性。

图8

小编总结

作者建立了一个基于三个m6A相关基因(VIRMA、IGF2BP1和HNRNPA2B1)的风险评分模型，作为OC的潜在预后生物标志物。同时还进行了ROC曲线、列线图、单变量和多变量Cox回归分析，以确认风险评分模型的准确性，表明该模型可以可靠地用于预测OC患者的预后。还构建了miRNAs-m6A调控基因-m6A靶基因的调控网络，其中一条通路miR-196b-5p-IGF2BP1-PTEN首次在OC中得到验证，表明其作为未来OC中潜在的生物标志物或治疗靶点的实用性。