m 6ARNA甲基化调节基因预测甲状腺癌的临床预后的作用

吴洪渊，王媛媛，傅 友，戴 萌，吴依芬，王 刚

南方医科大学南方医院1放疗科，2健康管理科，广东 广州510515；南方医科大学附属东莞市人民医院3肿瘤内科，4放射科，广东 东莞523059

甲状腺癌是一种最常见的内分泌相关的恶性肿瘤，超过90%内分泌系统肿瘤都是甲状腺癌，其发病率在逐年上升[1-5]。甲状腺癌虽然不是一种高发肿瘤，但仍然是女性第六大肿瘤，仅次于乳腺癌、肺癌、直肠癌、宫颈癌和皮肤恶性肿瘤，女性的发病率是男性的2～3倍[6-7]。甲状腺癌死亡率较低，其10年生存率90%～95%[8]。虽然手术、放疗、左甲状腺素治疗方法在不断进步[9-10]，仍有超过10%甲状腺癌患者会复发[11-14]，其余治疗手段对于局部进展和远处转移的甲状腺癌效果有限，因此进一步明确甲状腺癌的分子机制有助于肿瘤预后的预测及肿瘤的精准治疗。

目前THCA的发病机制未完全阐述清楚。在基因水平对甲状腺癌进行研究发现，RET、BRAF等基因在甲状腺癌的发生、发展以及转移中发挥重要作用[15]。m6A作为mRNA中最常见的甲基化修饰，无论是在mRNA还是lncRNA，都存在大量的m6A修饰。m6A对癌症的作用反应在对癌症相关基因表达的调节中，可促进肿瘤细胞的生长及侵袭，亦可产生抑制作用。m6ARNA调节基因METTL3与乳腺癌、卵巢癌和肝癌的发生发展关系紧密[16-18]。此外，FTO和ALKBH5在乳腺癌中同样发挥重要作用[19-20]。虽然越来越多的研究证明m6ARNA甲基化在众多肿瘤的发生发展中的作用，但是其与THCA的关系仍不清楚，也没有可依靠m6A构建的预后模型。本研究使用TCGA数据分析了THCA与正常组织中，20个m6ARNA甲基化调节基因的表达差异，用Lasso回归分析选出4个基因构建风险模型。然后验证其在THCA中的预测作用。本研究发现m6ARNA调节基因在THCA的发生发展中起着重要的作用，并且可作为独立预测基因来预测THCA病人的预后，为进一步个体化治疗提供预测条件。

1 资料与方法

1.1 数据下载

从TCGA数据库下载RNA-seq表达数据（https://portal.gdc.cancer.gov/），纳入58个正常组织，497个肿瘤组织。患者临床信息（表1）。

1.2 m6A甲基化调节基因选择

选择20个被公认的m6A调节基因，从TCGA数据库中得到表达数据，在正常组织和肿瘤组织中用limma包用wilcox检验方法进行差异基因筛选，然后将临床病理资料进行比较。

1.3 生物信息学分析

所有生物信息分析都在R语言软件中进行。热图绘制采用pheatmap包，用corrplot包用Pearson相关性分析对20个m6A调节基因进行相关性分析。将20个m6A调节基因用glmnet包进行lasso风险回归分析，选取4个与预后相关基因，并根据风险值将病人分为高低风险组，用survival包进行高低风险组生存分析及survivalROC包进行ROC分析检验模型的有效性，Cox回归分析用来分析风险值与年龄、性别、stage分期、T分期、N分期之间的关系。森林图用forestplot包进行绘制。为了进行蛋白水平验证，进入The Human Protein Atlas（https://www.proteinatlas.org/）进 行 IGF2BP1、IGF2BP2、FTO与RBM15在正常组织和肿瘤组织中的免疫组化表达情况探索。

1.4 预后模型验证

随机选取原数据集中的一半患者，根据风险值将其分为高风险组和低风险组，在这两组中进行OS及ROC分析，并进行Cox分析评估风险值与年龄、性别、stage分期、T分期、N分期之间的关系。

1.5 统计学分析

用wilcoxon test统计方法分析肿瘤样品和正常样品中m6ARNA甲基化调节基因的差异表达。区分病人为高低风险的cutoff值为上面得到的风险值。采用Kaplan-Meier方法分析OS，利用卡方检验来分析风险值和临床特征变量。最后构建出的风险模型在临床亚组中进行生存分析验证。

2 结果

2.1 THCA中m6A调节基因的表达

对20个m6A调节基因进行分析，结果显示，17个m6A RNA调节基因都有表达差异。IGF2BP2、IGF2BP3、HNRNPC表达上调，而YTHDC2、YTHDC1、ALKBH5、HNRNPA2B1、ZC3H13、FTO、METTL14、WTAP、YTHDF1、RBM15、KIAA1429、METTL3、YTHDF3、IGF2BP1表达下调（图1）。

2.2 m6A甲基化调节基因相关性

20个调节基因中除了IGF2BP1和IGF2BP3，剩下18个基因都具有相关性，其中相关性最大的是YTHDF3和KIAA1429（图2）。

2.3 m6A调节基因在THCA中的预后作用

20个基因进行lasso回归分析后筛选得到4个基因，分别是IGF2BP2（-9.95183573118394e-05）；

FTO（0.000109707439996991）；

IGF2BP1（0.0033484012761759）；

RBM15（0.00063213881736861）（图3A、B）。根据风险值的中位值，将病人分为高危组和低危组，高风险组OS与低风险组相比明显较低（图3C，P=4.278e-02）。ROC曲线用来评估预测模型的敏感性和特异性，结果显示AUC值为0.783（图3D）。

为验证所筛选出的4个基因在正常组织和肿瘤组织中的蛋白水平的表达，从The Human ProteinAtlas官网进行四个基因检索，结果显示：IGF2BP1在甲状腺正常组织和肿瘤组织中均为弱表达；IGF2BP2在甲状腺正常组织中表达高于肿瘤组织；FTO在正常组织和肿瘤组织中均高表达；RBM15未找到相关表达情况（图3E～J）。

2.4 THCA病人中风险值的作用

高低风险组病人临床信息（表1）。热图显示选择出来的4个基因在高低风险组中的表达量及与临床特征的相关性（图4A）。结果显示，在stage T分期（P＜0.05）和N分期（P＜0.01）中差异有统计学意义。单因素分析显示，年龄（以60岁划分为两组）、stage、T、riskScore与OS相关（P＜0.001）；多因素分析显示，年龄和riskScore与OS相关（P＜0.001，图4B、C）。

为进一步研究riskScore在临床亚组中的预测价值，根据病人临床特征进行分组（stage、T、N、M、age），高风险和低风险组的OS在N0（P=4.911e-02），stage3/4（P=4.322e-02）亚组中有差异，且高风险组OS更低。其他分组未见明显异常。

2.5 预测模型的验证

随机选中251例患者纳入为验证组，根据风险值将患者分为高风险组（n=121）、低风险组（n=130）。生存分析显示，与低风险组相比，高风险组的OS更短（P=8.756e-03，图6A），ROC显示在验证组中AUC值为0.778（图6B）。在验证组中，单因素分析表明年龄、stage分期、T分期和风险值与OS明显相关（图6C）。

图1 m6ARNA甲基化调节基因在THCA中的表达Fig.1 The expression of m6A RNA methylation regulators in THCA The heat map(A)and The violin plot(B)of 20 m6A RNA methylation regulators between THCA tumor samples and normalsamples.*P＜0.05,**P＜0.01,***P＜0.001.

图2 THCA中20个m6ARNA甲基化调节基因的相关性Fig.2 The correlation between 20 m6A RNA methylation regulators in THCA.

3 讨论

多组学研究表明THCA是一种分子水平高度异质性的疾病[21-22]，分子表达可作为一种THCA诊断、预后及治疗反应的潜在预测因素，可克服临床和组织病理学特征的局限性。mRNA、miRNA和lncRNA已被广泛探索，很多mRNA、miRNA和lncRNA相关的标志已已有研究证明可提高THCA诊断和预后[23-28]。多种基因总是相互作用来调节肿瘤的发生发展[29]。最近几年，m6ARNA甲基化正成为不同肿瘤研究领域中的热点，其与mRNA、miRNA和lncRNA之间的作用在很多肿瘤中也被广泛研究。m6ARNA甲基化调节基因在肿瘤的发生发展中起着重要的作用[30]，但其在THCA中的作用机制不明，因此很有必要探讨m6ARNA调节基因对THCA的作用。

图3 THCA中20个m6ARNA甲基化调节基因的风险特征Fig.3 The risk signature of 20 m6ARNA methylation regulators in THCA.A,B:The coefficients calculated by Lasso Cox regression analysis;C:Kaplan-Meier overall survival curve for THCApatients between high risk and low risk group;D:ROC curve showed the predictive value of risk signature.E,G,I:The IH Cexpression of IGF2BP1、IGF2BP2、RBM15 innormaltissue;F,H,J:The IH Cexpression of IGF2BP1、IGF2BP2、RBM15intumortissue.

有研究表明，在甲状腺癌细胞，METTL3表达升高，可促进THCA的发生发展，但其未在组织水平中验证[31]。本研究分析了m6ARNA甲基化调节基因THCA和正常组织中的表达，发现IGF2BP2、IGF2BP3、HNRNPC表达上调，而RBM15、METTL14、YTHDC1、FTO、ZC3H13、IGF2BP1、YTHDF3、YTHDC2、WTAP、HNRNPA2B1、 KIAA1429、 METTL3、 ALKBH5、YTHDF1表达下调。METTL3的表达水平与上述细胞学研究水平有所区别。表达上调的基因中，IGF2BP2和IGF2BP3均有文献报道。IGF2BP2促进导致结直肠癌的发病机理和发展[32]。在胰腺癌中，IGF2BP3表达较高，且与患者的整体生存有关[33]。在卵巢透明细胞癌中，过表达的IGF2BP3是一个潜在的致癌基因[34]。肺腺癌中，HNRNPC表达上升[35]。表达下调的基因中，部分基因的文献报道与本文一致。ALKBH5在结肠癌中被下调[36]。YTHDC2在直肠癌中是保护基因[37]。METTL3和YTHDF1与肺腺癌更好的OS和RFS相关[35]。过量表达METTL14在体外抑制HepG2细胞的迁移和侵袭，并在体内抑制细胞增殖和转移[38]。IGF2BP1基因表达情况在不同文献报道自相矛盾。在HCC中，IGF2BP1抑制细胞增殖并促进细胞凋亡[39]。而在卵巢癌中，IGF2BP1升高且预后不良[40]。KIAA1429、RBM15 、HNRNPA2B1、FTO的文献报道结果与本文相反。KIAA1429、RBM15 、HNRNPA2B1在肺腺癌中高表达[35]。同样，FTO在乳腺癌中高水平表达，并与其不良预后相关[20]。究其原因，这可能与肿瘤异质性有关，同一基因对不同癌症的影响不同。

表1 风险分组临床特征(n,%)Tab.1 The clinical in formation in different risk group

图4 风险评分和临床特征对THCA患者预后的影响Fig.4 Effects of the risk signature and clinical information on the prognosis of CESC patients.A:The heat map showed the different expression of m6A RNA methylation regulators and the distribution of clinical information between high risk group and low risk group;Univariate Cox regression(B)and Multivariate Cox regression (C) analysis of the correlation between clinical information and overall survival.

图5 risk Score在THCA患者临床亚组中的预测价值Fig.5 Prognostic value of the risk signature in THCA patients stratified into different clinical subgroup.Kaplan-Meiersurvival for patientin(A)N0and(B)stage3/4.

图6 在验证集中验证预后模型Fig.6 Validation of the prognostic signature in tested group.A:Kaplan-Meier overall survival curve for THCA patients between high risk and low risk group in tested group;B:ROC curve showed the predictive value of risk signature in tested group;C:Univariate Cox regression analysis the correlation between clinical information and overall survival in tested group.*P＜0.05.

m6ARNA甲基化调节基因对癌症的早期诊断和治疗产生重要影响，与癌症的预后值也紧密联系[41]。经Lasso cox回归分析，有4个m6ARNA甲基化调节基因IGF2BP2、FTO、IGF2BP1、RBM15，被筛选出以构建风险模型。利用ROC分析预测风险评分，AUC为0.713，说明本研究的风险模型构建良好。随后按风险评分分为高危组和低风险组的CESC患者，经单变量和多变量cox回归分析显示，风险评分是一个独立的甲状腺癌预测模型因素。此外，我们进一步评价了临床特征因素中的风险评分，结果发现高风险和低风险组的OS在N0和stage3/4亚组中有显著性差异，且高风险组OS均明显更低。当然本研究有一定的局限性，即只分析了TCGA中的THCA数据，没有进行临床实验验证，这将是后期研究重点。

综上所述，m6ARNA甲基化调节基因与THCA临床特征密切相关。本研究选择了4个m6ARNA甲基化基因，构建了一个可以作为独立预后因子的风险模型，为THCA中m6ARNA甲基化的后续研究提供了重要证据。