喉鳞状细胞癌m6A相关基因预后预测模型的构建及分析

韩弈垣 王雪梅 贺 晴 曹晓林

喉癌是常见的头颈部恶性肿瘤，喉鳞状细胞癌(laryngeal squamous cell carcinoma，LSCC)是其最常见的病理类型，约占总发病数的95%[1]。我国喉癌发病数和死亡数占癌症总发病和总死亡的比例分别是 0.67%和0.57%，标化发生率和病死率分别为1.3/10万和0.7/10万，且男性发生率和病死率远高于女性[2,3]。在过去的10年中，尽管LSCC的诊断技术及治疗方案不断改进，其5年总生存率仍从66%下降至63%[4]。因此，寻找影响LSCC患者预后的生物靶点和分子机制越来越受到重视。

N6-甲基腺苷(N6-methyladenosine，m6A)是在真核生物中mRNA和lncRNA甲基化修饰的常见形式，受甲基转移酶、脱甲基酶以及结合蛋白的动态调节，并参与许多重要的生理过程，例如干细胞分化和多能性，昼夜节律，胚胎发生和DNA损伤反应等[5]。目前已知的m6A相关基因包括YTHDC1、YTHDF2、YTHDC2、YTHDF1、FTO、ALKBH5、HNRNPC、METTL14、ZC3H13、KIAA1429、RBM15、WTAP和METTL3[6]。近年来有研究报道，异常的m6A修饰与头颈部鳞状细胞癌的发生及预后密切相关[7,8]。因此，本研究旨在探寻m6A相关基因在LSCC组织中的差异表达情况，并在此基础上构建出LSCC的预后预测模型。

材料与方法

1.数据下载：从TCGA数据库下载RNA-seq数据及相关临床数据，筛选标准如下：①肿瘤原发部位：喉；②病理类型：鳞状细胞癌；③肿瘤数据来源：TCGA-头颈部鳞状细胞癌；④数据类型：RNA-seq；⑤具备完整的临床信息，包括生存时间、生存状态、年龄、性别、肿瘤分级、临床分期、肿瘤范围、淋巴结转移、远处转移。经筛选共纳入123例RNA-seq数据样本(包括111例肿瘤样本及12例正常样本)及92例具有完整临床资料的LSCC患者临床数据样本。

2.m6A相关基因差异表达分析：利用Perl10.0.0软件及ensembl数据库将样本中的ensembl id转换为基因名，并整合各样本RNA-seq数据，提取所有样本中m6A相关基因的表达量。利用R4.0.3软件中的Wilcoxon秩和检验计算肿瘤组织及正常组织中m6A相关基因(包括YTHDC1、YTHDF2、YTHDC2、YTHDF1、FTO、ALKBH5、HNRNPC、METTL14、ZC3H13、KIAA1429、RBM15、WTAP和METTL3)差异表达情况，并分别用pheatmap包、vioplot包进行热图及小提琴图的绘制。以P<0.05为条件筛选LSCC样本中显著差异表达的m6A相关基因。

3.预后预测模型的构建及评估：为进一步进行变量选择，避免m6A相关基因过度拟合，通过R4.0.3软件中的survival包及glmnet包对显著差异表达的m6A相关基因进行Lasso回归分析及交叉验证，并构建风险评分模型。根据该模型计算临床样本中各组患者的风险评分，并以风险评分中位值将患者分为高风险组和低风险组。利用survival包绘制Kaplan-Meier生存曲线，比较两组间LSCC患者5年生存差异并进行Log-Rank检验。然后用survivalROC包绘制时间依赖的受试者工作特征曲线(receiver operating characteristic，ROC)评估此模型的预测性能。

4.独立预后分析：通过R4.0.3软件中的survival包及forestplot包对该预后预测模型进行单因素和多因素COX回归分析，并绘制森林图，从而评估各临床变量(包括年龄、性别、肿瘤分级、临床分期、T分期、N分期和M分期)及该模型的风险评分与患者预后的关系，进而判断该模型能否作为独立的预后指标。

结 果

1.m6A相关基因差异表达分析：LSCC组织样本和正常组织比较后，共筛选出9个显著差异表达的m6A相关基因，包括METTL3、YTHDF1、YTHDC2、ALKBH5、HNRNPC、RBM15、WTAP、FTO和KIAA1429。除了YTHDC2在LSCC组织中表达下调以外，其余8个基因在LSCC组织中表达均上调(图1、图2)。

图1 m6A相关基因差异表达的热图与正常组织比较，*P<0.01，**P=0.000

图2 m6A相关基因差异表达的小提琴图

2.预后预测模型的构建：对上述9个显著差异的m6A相关基因进行Lasso回归分析，并使用交叉验证建立模型，避免m6A相关基因过度拟合，结果显示，当变量数为9时，模型的偏差值最小(图3、图4)。该预后预测模型的风险评分公式如下：risk score=(KIAA1429表达值×0.0757)+(HNRNPC表达值×0.0486)-(FTO表达值×0.1308)-(WTAP表达值×0.0055)-(RBM15表达值×0.2427)-(ALKBH5表达值×0.0015)-(YTHDC2表达值×0.1660)-(YTHDF1表达值×0.0444)-(METTL3表达值×0.2598)。

图3 Lasso回归算法筛选变量动态过程图

图4 交叉验证参数λ选择过程图

3.预后预测模型的评估：整理从TCGA数据库下载的临床样本数据，筛除信息不全的临床样本，共得到92例具有完整临床资料的LSCC患者临床数据样本(表1)。根据风险评分公式计算92例LSSC患者的风险评分，并以风险评分中位值将患者分为高风险组和低风险组，其中高风险组36例，高风险组56例。Kaplan-Meier生存曲线结果显示，高风险组的总体生存率显著低于低风险组，且两组的5年生存率之间比较差异有统计学意义(P<0.01，图5)。时间依赖的ROC曲线分析结果显示该模型在预测LSCC患者预后生存期方面具有良好的效果(AUC=0.728，图6)。

表1 92例LSCC患者临床资料信息

图5 高、低风险组患者Kaplan-Meier生存曲线

图6 预后预测模型的ROC曲线及AUC值

4.独立预后分析：单因素及多因素COX回归分析表明：m6A相关基因预后预测模型的风险评分与LSCC患者的总体生存率显著相关(P=0.000)，且可作为独立的预后因素(P=0.000)。此外，性别和淋巴结转移情况同样与LSCC患者预后显著相关(P<0.05)，其中性别也可作为独立的预后因素(P=0.000，图7、图8)。

图7 各临床特征及风险评分的单因素COX回归分析森林图

图8 各临床特征及风险评分的多因素COX回归分析森林图

讨 论

LSCC的发生、发展涉及复杂的分子生物学改变，包括多种原癌基因的激活及抑癌基因的抑制，其具体的分子机制尚未完全阐明。近年来越来越多的研究证实异常的m6A修饰与头颈部鳞状细胞癌的发生与发展有关，并可依此来预测患者的生存预后，指导后续个体化治疗方案的制订[8～10]。

为了探寻m6A相关基因在LSCC组织中的差异表达情况，笔者从TCGA数据库下载了LSCC组织的基因表达情况及相关临床数据，并用Perl10.0.0软件及R4.0.3软件对数据进行了预处理及差异表达分析。笔者发现大多数m6A相关基因在LSCC组织中异常表达，包括8个显著上调的基因(METTL3、YTHDF1、ALKBH5、HNRNPC、RBM15、WTAP、FTO和KIAA1429)和1个显著下调的基因(YTHDC2)。有研究表明，METTL3高表达可导致肿瘤抑制因子SOCS2的m6A修饰增加，诱导SOCS2降解，促进肝细胞肝癌的发生；而KIAA1429也同样在肝癌组织中高表达，且与肝癌患者的预后相关[11,12]。YTHDF1可通过Wnt/β-catenin途径参与肠道恶性肿瘤的发展；同时也有研究证实YTHDF1与头颈部鳞状细胞癌患者的铁代谢相关，并可能成为靶向治疗的生物学标志物[10,13]。ALKBH5在胶质母细胞瘤和乳腺癌中表达上调，且与预后密切相关[14,15]。HNRNPC可通过EMT途径促进口腔鳞状细胞癌的恶性表达[16]。RBM15可作用于c-Myc抑制巨核细胞增殖，参与急性巨核细胞白血病的发展；而FTO、WTAP的异常表达可促进急性髓系白血病的发生[9,17,18]。YTHDC2作为抑癌基因，在头颈部鳞状细胞癌中显著下调，且与多种免疫细胞的免疫浸润水平相关[19]。

近年来，基于m6A相关基因的预后预测模型被用来评估和预测各种肿瘤患者的预后，指导后续个体化治疗方案的制定[20～22]。本研究通过对9个在LSCC组织中显著差异表达的m6A相关基因进行Lasso回归分析，构建了LSCC的m6A相关基因的预后预测模型。Kaplan-Meier生存曲线显示高风险组5年生存率明显低于低风险组，且时间依赖的ROC曲线结果显示该模型在预测LSCC患者预后生存期方面具有良好的效果(AUC值=0.728)。后续研究通过单因素及多因素COX回归分析证明该模型可作为独立的预后因子(P=0.000)。此外，本研究还发现性别也可作为独立的预后因子(P=0.000)，且男性患者的预后欠佳。这可能与男性更多地暴露于烟酒等危险因素及性激素的复杂作用有关[23]。

综上所述，本研究构建了基于9个m6A相关基因预测LSCC患者预后的风险评分模型。该模型具有较好的敏感度和特异性，且其可以作为独立的预后因子评估LSCC患者预后，对制定合理及有效的个体化治疗方案具有重要的参考价值。但该模型的应用仍需要大量临床研究数据支持和大规模、多中心的循证医学证据加以验证。