m6A甲基化调节因子影响口腔鳞状细胞癌生存预后的生物信息学分析*

徐 莉， 余 程， 杜习金

华中科技大学同济医学院附属同济医院口腔医学中心，武汉 430030

口腔癌是最常见的恶性肿瘤之一，其5年生存率仅41%～59.2%[1-2]，2018年全球新增口腔癌患者超过35万[3]。口腔鳞状细胞癌(oral squamous cell carcinoma，OSCC)作为口腔癌最主要的病理类型，研究它的发生、发展和预后，有助于降低口腔癌整体的发病率并改善患者的预后。OSCC发病的分子机制仍不确定，研究有效的预后标志物和治疗靶点一直是OSCC研究的热点。m6A甲基化修饰在多种肿瘤的发生和发展中起重要作用[4-5]，由甲基转移酶、去甲基化酶和识别蛋白等m6A甲基化调节因子调控[6]，调节因子参与乳腺癌、白血病和胰腺癌等多种癌症的进展过程[7-9]。m6A甲基化调节因子可能参与OSCC的进展过程，影响其预后。目前国内外尚未见m6A甲基化调节因子与OSCC相关性研究的报道。

生物信息学方法能分析癌症基因组图谱(The Cancer Genome Atlas，TCGA)数据库中m6A甲基化调节因子的基因，利用相关临床数据获得基因的表达水平与患者生存预后的关系。本研究通过TCGA数据库中口腔鳞状细胞癌level 3的RNA-seq数据，应用生物信息学方法分析m6A甲基化调节因子对OSCC患者生存预后的影响，探讨相关的重要信号通路，筛选出与预后相关的基因，为OSCC的治疗和预后评估提供有价值的分子靶点，为OSCC的研究开辟新的思路。

1 材料与方法

1.1 TCGA数据收集

从TCGA数据库选择OSCC病例，日期截止2019年7月17日，选择头颈部肿瘤中颊、口底、舌、牙龈、唇、牙槽嵴和硬腭部位的肿瘤。下载病例的RNA-seq数据，生成矩阵文件。接下来通过ensembl数据库，将基因名由ensembl id转换成gene symbol的矩阵。根据样本信息，拆分出258个样本，其中包括241个肿瘤样本和17个正常组织样本。剔除临床信息不完善及重复测序的样本，本研究对239个OSCC样本进行研究。

1.2 选择m6A甲基化调节因子及其相关性分析

从gene symbol矩阵提取m6A相关基因的表达，得到13个m6A甲基化调节因子的表达量，包括：甲基转移酶样因子3(methyltransferase 3，METTL3)、甲基转移酶样因子14(methyltransferase 14，METTL14)、Wilm’s肿瘤1相关蛋白(Wilm’s tumor 1-associated protein，WTAP)、KIAA家族基因1429(KIAA1429)、RNA结合基序蛋白15(RNA binding motif protein 15，RBM15)、CCCH型锌指结合蛋白13(zinc finger domain-containing protein 13，ZC3H13)、脂肪量和肥胖相关蛋白(fat mass and obesity-associated protein，FTO)、AlkB同源物5(alkB homolog 5，ALKBH5)、YTH结构域结合蛋白1(YTH domain-binding protein 1，YTHDF1)、YTH结构域结合蛋白2(YTH domain-binding protein 2，YTHDF2)、YTH结构域包含蛋白1(YTH domain-containing 1，YTHDC1)、YTH结构域包含蛋白2(YTH domain-containing 2，YTHDC2)以及异质核糖蛋白C(heterogeneous nuclear ribonucleoprotein C，HNRNPC)。然后，corrplot R包对这些基因进行相关性分析，探讨它们间的相互作用。

1.3 样本分型

以“ConsensusClusterPlus”R包对13个m6A甲基化调节因子的基因进行样本聚类分析，参数设定如下：iterations是50，resample rate是80%，分析方法是Pearson相关性分析。选择k=2对样本分型，将样本分为2组。主成分分析算法(principal component analysis，PCA)比较不同聚类分型m6A甲基化调节因子的基因表达谱。

1.4 生存分析及热图

绘制生存曲线，便于观察患者的生存期。使用survival R包对不同聚类分型的样本做Kaplan-Meier分析绘制生存曲线，采用Log-rank方法进行不同分型间差异显著性检验。pheatmap R包对m6A甲基化相关基因按照聚类分型进行分析，绘制成热图。同时，在热图上方加入临床性状，并使用卡方检验计算临床性状在不同聚类分型间差异的显著性。

1.5 单因素Cox分析

为了能更方便且迅速地评估OSCC的预后，我们筛选出预后相关基因。使用Survival包对m6A甲基化调节因子的基因做单因素Cox分析，筛选条件P<0.05，获得与OSCC预后显著性相关的基因。

1.6 基因集富集分析

对聚类分型的两组样本做基因集富集分析(GSEA)，因为分型是根据13个m6A甲基化调节因子的基因得到的，所以相当于是基于以上调节因子的基因进行富集分析。首先，我们选取基因集“c5.all.v6.2.symbols.gmt”做GO富集分析，得到GO富集分析的结果。然后选取基因集“c2.cp.kegg.v6.2.symbols.gmt”做KEGG富集分析，得到KEGG富集分析的结果。

1.7 统计学分析

采用Rv3.5.2软件进行统计学处理；GSEA富集分析采用GSEA 4.0.0软件处理；生存分析采用Kaplan-Meier分析方法，其组间差异比较采用Log-rank检验；临床性状的组间差异比较采用卡方检验；调节因子的相关性分析以P<0.001为差异具有统计学意义，其它数据的分析以P<0.05为差异具有统计学意义。

2 结果

2.1 样本分型

对239个OSCC样本聚类分型，得到cluster1(n=84)和cluster2(n=155)两个分组。每个样本的分型保存后，得到PCA的结果，从图1可以看出，cluster1和cluster2可以明显分开。使用m6A甲基化调节因子可以将OSCC样本分为两型。

图1 cluster1和cluster2的PCA分析结果Fig.1 PCA analysis results of cluster1 and cluster2

2.2 生存分析及热图

从生存曲线(图2)可以看出，两组生存期的差异显著(P<0.05)。cluster1和cluster2组患者5年生存率分别为36.0%(95%CI：23.6%～54.9%)和55.7%(95%CI：44.8%～69.1%)。从热图中可以看出，OSCC的肿瘤级别与聚类分型显著相关(图3)。这说明分型影响OSCC的预后并与其恶性程度相关，也就是m6A甲基化调节因子影响OSCC的预后并反映其恶性程度。

图2 cluster1和cluster2的生存曲线Fig.2 Overall survival curves of cluster1 and cluster2

红色：基因高表达；绿色：基因低表达；基因后面的上标“a”代表0.001≤P<0.01，无上标代表P≥0.05图3 m6A甲基化调节因子热图及组间性状比较Fig.3 m6A methylation regulators’ heatmap and comparison between two clusters

2.3 m6A甲基化调节因子基因间的相关性

调节因子基因间的相关性结果如图4。WTAP与HNRNPC、RBM15、YTHDF2、METTL14及YTHDC1相关；HNRNPC、YTHDF2和METTL14互相相关。ALKBH5与FTO和YTHDF1均正相关，与其它基因不相关；WTAP、HNRNPC、YTHDF2、METTL3与YTHDF1不受FTO影响。

2.4 生物学功能

GO功能分析显示，m6A甲基化调节因子基因所富集的功能有1331个，根据校正后P值(FDR)选取富集程度最高的前10条(图5，表1)，表明其主要参与mRNA结合、核运转、核蛋白复合体的生物合成等过程。

蓝色：负相关；红色：正相关；颜色越深，则相关程度越高；打叉的圆圈：相关性不显著(P>0.001)图4 m6A甲基化调节因子基因间的相关性Fig.4 Correlation of m6A methylation regulators

上半部分是富集得分富集的过程，曲线上升，直到到达最大的富集得分值，此时的得分值就是该功能富集的富集得分值；下半部分是基因在每个功能的分布情况；不同的颜色代表不同的功能。图5 m6A甲基化调节因子基因的GO功能富集图Fig.5 GO function enrichment of m6A RNA methylation regulator genes

表1 GO功能富集程度最高的前10条Table 1 The top 10 of GO function enrichment

2.5 KEGG信号通路

KEGG信号通路富集发现，m6A甲基化调节因子基因所富集的通路有31个，根据校正后P值(FDR)选取富集程度最高的前10条(图6，表2)，表明其主要涉及剪接体、细胞周期和同源重组等KEGG通路。

上半部分是富集得分富集的过程，曲线上升，直到到达最大的富集得分值，此时的得分值就是该功能富集的富集得分值；下半部分是基因在每个功能的分布情况；不同的颜色代表不同的功能。图6 m6A甲基化调节因子基因的KEGG信号通路富集图Fig.6 KEGG pathway enrichment of m6A methylation regulator genes

表2 KEGG信号通路富集程度最高的前10条Table 2 The top 10 of KEGG pathway enrichment

2.6 预后相关基因

根据单因素Cox分析结果(表3)，有3个m6A甲基化调节因子的基因与生存期显著相关(均P<0.05)，可以作为OSCC的预后相关基因。这3个基因是YTHDF2、HNRNPC和METTL14，它们都是OSCC的风险基因，其表达越高则OSCC预后越差。

表3 OSCC中m6A甲基化调节因子的预后价值Table 3 The prognosis value of m6A methylation regulators in OSCC

3 讨论

OSCC的分子机制尚未明确，原癌基因的活化或多个肿瘤抑制基因的失活是促发OSCC的主要原因[10-11]，目前没有评估其预后的最佳方法。TCGA数据库包括人类50种肿瘤(含亚型)的基因组变异图谱。本研究通过生物信息学技术挖掘TCGA数据库中的信息，探讨m6A甲基化调节因子对OSCC预后的影响。我们对13个m6A甲基化调节因子聚类分型，得到cluster1和cluster2两组OSCC样本，它们在生存期和肿瘤级别间差异显著(均P<0.05)，说明m6A甲基化调节因子影响OSCC的预后并反映其恶性程度。m6A调节因子间具有复杂的相关性，相互作用，影响OSCC的进程。本研究筛选出了3个预后相关的调节因子YTHDF2、HNRNPC和METTL14，其表达增高是OSCC预后不良的因素。目前尚未见它们在OSCC中的作用相关研究的报道，而它们对其它恶性肿瘤的作用已被部分证实，譬如：METTL14在急性髓细胞性白血病中起到癌基因的作用[12]；YTHDF2在胰腺癌组织中上调，有望用于诊断胰腺癌并评估其预后[9]；HNRNPC能促进乳腺癌细胞的增殖和肿瘤的生长[13-14]。

GSEA富集分析探讨m6A甲基化调节因子的生物学功能和信号通路。GO富集分析发现它们在1331个功能中富集，m6A甲基化基因的功能多且复杂，主要参与mRNA结合、核运转、核蛋白复合体的生物合成等，影响遗传物质的复制、转录和翻译过程，参与OSCC细胞发生和发展的恶化过程。KEGG通路富集分析显示m6A甲基化调节因子在31个通路中富集，主要包括剪接体、细胞周期和同源重组等相关的KEGG信号通路。剪接体在前体mRNA的选择性剪切过程中发挥作用[15]，同一个前体mRNA经过选择性剪切可以形成不同的蛋白分子。发生在剪切位点的同义突变使剪接体无法识别内含子，导致内含子被保留，肿瘤抑制因子失活[16-17]，肿瘤发生。有研究证明剪接体在癌症细胞中发生明显的突变[18]。细胞周期调节因子主要有肿瘤蛋白53(P53)和细胞周期蛋白依赖性激酶(CDKs)酶家族。P53蛋白活化后通过调节基因的表达，调控细胞的增殖、衰老和死亡[19]，大约50%的肿瘤P53蛋白存在功能缺陷[20]。CDKs受细胞周期素的正调控和细胞周期蛋白依赖性激酶抑制蛋白的负调控；当CDKs发生异常时，会使细胞周期发生紊乱，细胞增殖和凋亡间的平衡失控，导致癌症的发生[21]。同源重组是修复DNA双链断裂损伤的主要方式，影响肿瘤细胞对放疗和基因药物治疗的敏感性，其功能缺陷会增加正常细胞癌变的可能性[22]。筛选出的这些通路在OSCC的形成与发展中起重要作用，进一步研究m6A甲基化调节因子的富集结果，有助于阐明OSCC发生和发展的分子机制。

本研究发现METTL14、YTHDF2和HNRNPC是OSCC的预后相关基因，有望成为OSCC治疗和预后评估的分子靶点，它们的表达越高则OSCC预后越差。m6A甲基化调节因子通过甲基转移酶“写入”、去甲基化酶“去除”以及识别蛋白“识别”甲基化修饰，发挥对甲基化修饰的调控作用。当m6A甲基转移酶存在时，发生m6A甲基化修饰，并由特异性识别蛋白识别后发挥生物学功能；若去甲基化酶催化则m6A甲基化修饰被去除[23]。通过对m6A甲基化调节因子的相关性进行分析，明确基因的相互作用。预后相关基因间两两相互作用；YTHDF2与YTHDC1、YTHDF1、YTHDC、KIAA1429、ZC3H13、WTAP和RBM15均正相关；HNRNPC与KIAA1429、METTL3、WTAP、RBM15和YTHDC1均正相关；METTL14与除ALKBH5外的其它基因均正相关。ALKBH5与预后相关基因均无相关关系，FTO仅与METTL14正相关，而FTO还参与N6，2-o-dimethyladenosine(m6Am)的去甲基化[24]，推测在OSCC的发生和发展过程中去甲基化酶作用不大，主要是甲基化酶和识别蛋白发挥作用，促使m6A甲基化修饰增高。另外，预后相关基因中只有METTL14一种识别蛋白，可见它在OSCC中的作用很关键。METTL14可以作为癌基因，通过m6A甲基化修饰来调控它的靶mRNA(如MYB和MYC)，它则受SPI1的负调控，形成SPI1-METTL14-MYB/MYC信号轴，参与骨髓细胞和白细胞的生成[12]。YTHDF2参与胰腺癌细胞的增殖、抑制癌细胞的迁移和浸润；敲除YTHDF2基因使YAP的表达增加，抑制TGF-β/Smad信号通路，调节上皮-间充质转化过程[9]。HNRNPC在RNA剪切[25]、序列特异性的RNA输出[26]和RNA的翻译[27]中起调节作用；抑制HNRNPC会导致乳腺癌中内源性双链RNA上调，活化下游的干扰素反应，抑制癌细胞的增殖[13]。继续研究YTHDF2、HNRNPC和METTL14，探讨它们间的相互作用，将为寻找治疗OSCC等恶性肿瘤的新方法提供帮助。

综上所述，本文应用生物信息学方法揭示了m6A甲基化调节因子对OSCC患者的生存预后产生影响，其中YTHDF2、HNRNPC和METTL14可以作为预后相关基因，它们的表达越高则OSCC预后越差，有望成为OSCC治疗和预后评估的分子靶点。后期将进一步通过体内或体外实验验证YTHDF2、HNRNPC和METTL14的表达增高与OSCC发生和发展的相关性，探讨它们在预后评估中的价值。