肝癌m6A RNA甲基化调控因子风险预测模型的构建及其临床意义

王进，毛德利，赵鑫，汤权，董晓强

1.苏州市独墅湖医院(苏州大学附属独墅湖医院)普外科，江苏 苏州 215000；2.苏州大学附属第一医院普外科，江苏 苏州 215006

肝癌是临床上常见的恶性肿瘤之一，发病率居全球第5位，恶性程度高，预后差，全球每年新发肝癌病人超过74万，每年因肝癌而死亡病人约70万[1]。因此，寻找可以判断肝癌的独立预后因子，对于肝癌病人的后期治疗寻找新的靶点及改善病人总体生存率有着至关重要的作用。RNA修饰是一种转录后水平的调控方式，其中RNA甲基化修饰是主要形式之一。mRNA是遗传信息到蛋白质的中间体，mRNA修饰直接导致蛋白质改变从而影响细胞功能,并且研究显示mRNA甲基化是可逆、动态的修饰过程。其中N6-甲基腺嘌呤(N6-methyladenosine，m6A)修饰是mRNA含量最丰富一种甲基化修饰形式[2]。m6A广泛参与调控生命周期的各个阶段，包括mRNA剪接、加工、翻译和降解等，在恶性肿瘤的进展中发挥重要作用[3-7]。

近年来，多项研究表明m6A调控因子突变与癌症的发生发展密切相关，其中甲基化调控因子METTL3缺失增加了胰腺癌细胞对化疗和放疗敏感性，但是对胰腺癌细胞增殖无明显影响[8]，Lin等[9]研究显示METTL3在结直肠癌与肺腺癌中高表达，通过促进表皮生长因子受体(EGFR)和河马信号效应器TAZ等mRNA翻译，并促进肺腺癌细胞生长、侵袭。甲基化修饰中沉默负责“擦除”的调控因子ALKBH5后可明显抑制恶性胶质瘤干细胞增殖。而ALKBH5通过使FOXM1 mRNA去甲基化来增强FOXM1表达[10]。Chen等[11]研究发现，在肝癌中METTL3被激活后甲基化修饰SOCS1 mRNA通过YTHDF2依赖性机制，并影响肝癌病人预后[11-12]。目前，尚无m6A调控因子与肝癌病人临床病理及预后关系的综合性报道。本研究从癌症基因组图谱(TCGA)数据库(n=424)和国际癌症基因组联盟(ICGC)数据库(n=445)中获得基因表达谱数据及相应的临床病理数据进行分析，证实了m6A调控因子在肝癌病人中存在差异性表达。利用单因素Cox比例风险回归分析方法筛选6个关键的m6A调控因子并建立风险预测模型，为肝癌病人的治疗及改善预后提供了重要的参考依据。

资料与方法

一、数据资料

50例正常组织和374例肝癌组织的基因表达谱数据和相应的临床资料获取自TCGA数据库(https://portal.gdc.cancer.gov/)，243例正常组织和202例肝癌病人基因表达谱数据和相应的临床资料获取自ICGC数据库(https://icgc.org/)。

二、m6A RNA甲基调控因子的挑选

修饰调控因子中主要负责“书写”的甲基转移酶是最主要的调控因子，包括KIAA1429[13]，METTL3[14-15]，METTL14[14-15]，RBM15B[16-17]，RBM15[16]，Wtap[13, 15]，ZC3H13[18]；负责“读取”的调控因子包括HNRNPC、HNRNPA2B[19]，YTHDF、YTHDF2、YTHDF3[20]，YTHDC1[21]，YTHDC2[22]，LRPPRC[23-24]，FMR1[24]；负责“擦除”的调控因子主要由FTO[25]，ALKBH5[26]组成。

三、m6A RNA甲基调控因子差异性表达及相关作用的分析

利用edgeR函数包对从数据库获取的正常组织和肝癌组织基因表达谱数据进行归一化处理，提取18种m6A RNA甲基化调控因子表达量，分析m6A RNA甲基化调控因子在肿瘤组织及正常组织表达的差异性，其中log2fold change(log2FC)绝对值>1，P<0.05的基因为差异表达基因。Corrplot R包分析调节因子在肝肿瘤组织中表达相关性，并计算相关系数。使用基因集富集分析(Gene Set Enrichment Analysis，GSEA)来探索调节因子在肝癌中的相关功能。

四、风险预后模型的建立及评估

为了确定m6A RNA甲基化调控因子在肝癌中的预后价值，使用单因素Cox回归分析方法筛选出12个与总体生存率相关基因(P<0.05)。使用 LASSO(the least absolute shrinkage and selection operator)算法从12个基因中确定了6个基因来建立风险预测模型。每位病人风险值等于6个基因中每个基因表达量乘于系数之和(Coefi是风险系数，Xi是基因表达量)。计算公式如下：

根据风险值将TCGA和ICGC数据库中病人分为高风险组与低风险组，分析不同风险组病人临床病理特征及总体生存率关系，并绘制Kaplan-Meier受试者工作特征(ROC)及校准曲线。为了评估风险值是否可以作为独立预后指标，我们对TCGA数据库进行单因素和多因素Cox回归分析，并在ICGC数据库进行验证。

五、统计学分析

Wilcoxon秩和检验比较正常组织与肿瘤组织基因表达差异性；皮尔逊(Pearson)检验用于计算两组基因相关系数；使用LASSO Cox回归算法来建立风险预后模型；高风险组和低风险组的总体生存率由Kaplan-Meier方法完成； 由ROC曲线检验风险评分预测的准确性；使用单变量和多变量Cox回归分析来评估风险值与临床病理特征在肝癌预后中的危险程度；Wilcoxon秩和检验用于比较不同风险组及临床分期之间风险基因表达差异性。所有统计分析由SPSS(19.0版) and R v3.6.0软件完成。P<0.05为差异有统计学意义。

结 果

一、m6A RNA甲基化调控因子在肝癌组织及正常组织中的表达差异性

TCGA数据库的基因表达谱数据分析结果显示：与正常组织相比，18个m6A RNA 甲基化调控因子在肿瘤组织中均高表达(log2FC>1)，除METTL14和ZC3H13之外，其余调控因子表达差异均有统计学意义。其中以下因子差异表达倍数较大(图1A)：METTL3(log2FC=1.352，P=7.37E-24)；KIAA1429(log2FC=1.233，P=5.59E-24)；RBM15B(log2FC=1.013，P=6.26E-23)。ICGC数据库的基因表达谱数据结果显示：与正常组织相比，肿瘤组织ZC3H13调控因子显著低表达(log2FC<-1)，其余调控因子均高表达，除METTL14之外其余表达差异均具有统计学意义。其中KIAA1429 (log2FC=1.153，P=4.14E-58)以及RBM15B (log2FC=1.01，P=3.96E-48)差异表达倍数较大(图1B)。ICGC 数据集基因表达量结果与TCGA数据集基本一致。

图1 癌症基因组图谱(TCGA)数据库(A)和国际癌症基因组联盟(ICGC)数据库(B)中m6A RNA甲基化调控因子在肝癌组织及正常组织的表达差异性

二、m6A RNA甲基化调控因子相互作用分析

通过Corrplot R包对TCGA数据库中m6A RNA甲基化调控因子进行相关性分析，结果显示调节因子之间有着复杂的联系，其中负责“擦除”的RBM15、METTL14、METTL3、KIAA1429和Wtap与其余调节因子之间彼此显著相关(图2A)。最后GSEA显示HNRNPC、 HNRNPA2B1和METTL3与肿瘤发生及进展有关，比如癌症[标准化富集得分(NES)=1.908, 标准化P<0.001]，ERBB信号通路(NES=1.925，标准化P<0.001)，RIG-I-like-受体信号通路(NES=2.073，标准化P<0.001)，Wnt信号通路(NES=2.070，标准化P<0.001)(图2B)。以上分析结果提示m6A RNA甲基化调控因子在肝癌进程中和生存预后发挥着重要作用。

图2 癌症基因组图谱(TCGA)数据库中m6A RNA甲基化调控因子相互作用分析 A.Corrplot R包分析m6A RNA甲基化调控因子的相关性；B.肿瘤发生和进展相关调控因子的基因集富集分析(GSEA)

三、建立预后风险模型

利用单因素Cox回归分析m6A甲基化调控因子对病人生存预后的影响。结果显示YTHDF2、LRPPRC、YTHDF1、METTL3、RBM15B、KIA-A1429、HNRNPA2B1、HNRNPC、WTAP、YTHDC1、RBM15、ZC3H13显著影响病人的总体生存率(图3A)，差异均有统计学意义(P<0.05)。其中YTHDF2、LRPPRC、YTHDF1、METTL3、RBM15-B、KIAA1429、HNRNPA2B1、HNRNPC、WTAP、YTHDC1、RBM15这些调控因子的风险比(HR)>1，属于高风险基因；而ZC3H13调控因子的HR小于1，属于低风险基因。使用LASSO算法，从以上12个基因中筛选出6个基因(KIAA1429、ZC3H13、METTL3、YTHDF1、YTHDF2、LRPPRC)，建立风险预测模型并计算风险值。根据中位风险值将肝癌病人分为高风险组与低风险组(图3B～D)。GSEA显示风险值与WNT信号通路及泛素介导的蛋白质水解密切相关(图3E、F)。

图3 在癌症基因组图谱(TCGA)数据库中筛选6个肝癌预后相关基因建立风险预测模型 A.单因素Cox回归分析m6A甲基化调控因子对病人的生存预后的影响；B、C.LASSO Cox回归分析筛选6个预后相关基因建立预后模型；D.被用来建立预后模型的6个基因系数；E、F.基因集富集分析(GSEA)显示风险值与Wnt信号通路及泛素介导的蛋白质水解密切相关

四、高、低风险组肝癌病人预后生存分析以及ROC曲线的绘制

使用Kaplan-Meier曲线对高风险组与低风险组预后生存分析显示，在TCGA与ICGC数据库中高风险组与低风险组总体生存时间差异有统计学意义，高风险组病人的5年生存率明显低于低风险组(P<0.001)(图4A)。同时使用ROC曲线检测风险值预测肝癌病人5年生存率的准确性，结果显示：曲线下面积(AUC)分别为61.9%(TCGA)与76.5%(ICGC)，提示风险预测模型可以有效预测肝癌病人5年生存率。两组数据AUC值相差较大可能是不同区域来源的病人存在的较大异质性导致。TCGA肝癌数据主要来源于美国，ICGC肝癌数据来自于日本，统计发现TCGA数据肝癌病人5年死亡比例占39%，ICGC数据肝癌病人5年死亡比例占18%(图4B)。校准曲线显示风险模型对肝癌病人5年总体生存率预测效果较好(图 4C)。

五、单因素与多因素Cox回归分析

对TCGA和ICGC数据库中肝癌病人临床病理特征及风险值进行单因素与多因素Cox回归分析。在TCGA数据库中，单因素分析显示分期(P<0.001)、T分期(P<0.001)、N分期(P<0.05)、M分期(P<0.05)、风险值(P<0.001)与病人总体生存率显著相关。将分期、T分期、N分期、M分期、风险值纳入多因素Cox回归分析，结果显示风险值(P<0.001)与病人总体生存率显著相关(P<0.001)(图5A)。在ICGC数据库中，单因素分析显示分期(P<0.001)、性别(P<0.05)、风险值(P<0.001)与病人总体生存率显著相关。将分期、性别、风险值纳入多因素Cox回归分析，显示分期(P<0.001)、性别(P=0.006)，风险值(P=0.001)与病人总体生存率显著相关。以上结果提示风险值可以作为肝癌的独立预后因子(图5B)。

注：AUC.曲线下面积；TCGA.癌症基因组图谱；ICGC.国际癌症基因组联盟。图4 风险值在肝癌在诊断及预后中的价值 A1、A2.不同风险组肝癌病人总体生存率；B1、B2.受试者工作特征(ROC)曲线判断风险预测模型对肝癌病人5年生存率预测的准确性；C1、C2.校准曲线检验风险预测模型对肝癌病人5年生存率预测的准确性

图5 在癌症基因组图谱(TCGA)数据库和国际癌症基因组联盟(ICGC)数据库中不同临床病理特征及风险值的单因素及多因素Cox回归分析 A1.TCGA数据库中单因素分析；A2.TCGA数据库中多因素分析；B1.ICGC数据库中单因素分析；B2.ICGC数据库中多因素分析

六、不同风险组肝癌病人临床病理特征及预后的差异性

分析在不同风险组中，临床病理特征、生存预后以及6个调节因子表达的差异性。结果显示YTHDF2、LRPPRC、YTHDF1、METTL3、KIAA1429在高风险组高表达，而ZC3H13在高风险组低表达。进一步分析显示：在TCGA数据库中，不同风险组之间分期(P<0.05)、生存时间(P<0.001)存在着显著性差异(图6A)。在ICGC数据库中，不同风险组之间分期(P<0.01)、分级(P<0.01)、年龄(P<0.05)、生存状态(P<0.05)存在着显著性差异(图6B)。

图6 在癌症基因组图谱(TCGA)数据库(A)和国际癌症基因组联盟(ICGC)数据库(B)中，不同风险组临床病理特征及预后的差异性

讨 论

肝癌是全球第五大癌症，具有高死亡率及侵袭性特点，占全球癌症死亡率的9.1%，肝癌的发病率和死亡率在逐年增加[27]。肝癌病人往往确诊时已晚期，同时术后易复发、转移、耐药，故死亡率较高。因此深入了解肝癌的发病机制，寻找新的诊断及预后标志物对临床诊断和治疗工作有重要的指导意义。本次研究通过对肝癌中m6A异常甲基修饰调控因子进行分析，证实了其在肝癌发生发展中的重要性。从中筛选出具有影响肝癌病人预后的6个调控因子，并建立风险预测模型，其在临床诊疗中具有很大程度的参考意义。

既往已有大量研究报道m6A RNA甲基化调控因子在多种癌症中发挥重要作用[8-9, 25]。在乳腺癌中，METTL3可通过m6A修饰促进HBXIP的表达，从而促进肿瘤增殖[28]。对METTL3进行的m6A甲基化修饰是急性髓细胞白血病发生所必须的方式[29]。METTL3、YTHDF1和YTHDF2甲基化调控因子已被证明影响肝癌病人预后，并且可能作为肝癌新的诊断治疗生物标志[11-12]。Cheng等[30]发现，在肝癌中负责“书写”的KIAA1429调控因子通过m6A修饰抑制癌基因ID2表达，从而导致肝癌侵袭力增加。在本次研究中，通过总结目前常见的18种m6A RNA甲基化调控因子，从TCGA和ICGC数据库中获取基因表达图谱及临床信息，并分析其在肿瘤组织及正常组织表达的差异性，发现其中负责“读写”的METTL3、KIAA1429和RBM15B调控因子表达显著高于其他调控因子。

既往研究显示，乙型肝炎病毒感染是肝癌的主要致病因素之一，乙型肝炎病毒感染通过下调肝癌病人中RIG-I-like受体信号从而抑制Ⅰ型干扰素的表达[31]。多种恶性肿瘤与ERBB信号失调相关，比如肺癌、卵巢癌、乳腺癌，ERBB信号失调往往导致更差的生存预后[32-33]。过去的研究认为ERBB信号是肝脏急性损伤时的保护机制，随着研究的深入，发现肝脏在慢性条件刺激下，ERBB信号参与肝癌的形成[34]。Wnt信号主要参与结直肠癌生物学过程，在肝癌中，β-连环蛋白N末端区域突变，导致β-连环蛋白稳定性增加，在肝脏组织中累积，进而促进肿瘤进展。在本次研究中，GESA结果显示m6A RNA甲基化调控因子与上述介导肝癌发生发展的多种信号传导通路具有很大程度上的相关性。以上的实验结果提示m6A RNA异常甲基化调控因子可能参与肝癌发生、发展。在肝癌诊断和预后中起到至关重要的作用。

在TCGA数据库中，利用Cox回归分析筛选12个预后相关基因，进一步的LASSO 算法确定了6个影响肝癌病人生存预后的关键调节因子，包括KIAA1429、ZC3H13、METTL3、YTHDF1、YTHDF2、LRPPRC。由此建立风险预测模型，根据中位风险值将肝癌病人分为低风险组与高风险组。分析发现高风险组与低风险组之间病人总体生存率差异有统计学意义，单因素与多因素Cox回归分析结果均显示风险值可以作为独立预后的指标。然而，与既往的报道[34]不同的是，本研究的分析结果显示，在TCGA数据库中N分期未表现出独立的危险因素，这可能与本研究纳入的标本量较小有关。为了验证所建立模型的有效性，本研究使用ROC曲线和ICGC测试集数据进行了检验，最终确定了6个m6A甲基化调控因子构建的风险预测模型对判断肝癌预后的可靠性。

综上所述，本研究通过生物信息学分析证实m6A RNA甲基化调控因子在肝癌肿瘤组织中与正常组织中表达具有显著性差异，m6A RNA甲基化调控因子的异常表达与肝癌临床病理及病人预后密切相关，基于筛选的6个关键m6A RNA甲基化调控因子建立的风险预测模型对指导肝癌后期个性化治疗及改善预后具有重要的参考价值。