LncRNA m6A修饰在肺鳞癌的预后及免疫疗效预测中的鉴定和价值分析*

顾 睿，莫春芬，林 苹，张 洁，李 凯△

1.四川大学华西医院 肿瘤中心实验肿瘤研究室(成都 610041); 2.成都医学院 基础医学院免疫学教研室(成都 610500)

据统计，2018年全球新增肺癌病例209万人，且肺癌的病死率逐年升高[1]。肺鳞癌作为肺癌主要的组织学亚型，发病率占所有肺癌的25%～30%[2]。肺鳞癌临床预后较差，常规放化疗及靶向治疗疗效欠佳[3]。免疫检查点抑制剂(immune checkpoint inhibitors，ICI)的应用对于无免疫治疗禁忌证的患者，无论是单一免疫疗法还是联合免疫疗法，都已成为晚期肺鳞癌的一线治疗标准[4]。因此，迫切需要可靠的生物标志物用以评估肺鳞癌患者的预后及免疫疗效，从而针对患者病情制定个性化的诊疗方案。LncRNA是一类转录本长度>200 nt的非编码RNA,通过表观遗传调控、转录水平调控等多种方式参与编码蛋白质基因的功能调控[5]。M6A修饰是一种可逆性的RNA表观遗传调控方式，即腺苷酸(A)在甲基转移酶的催化作用下，第6位N发生甲基化修饰，通过特定的阅读蛋白识别发生m6A修饰的碱基，激活下游的调控通路[6]。LncRNA m6A修饰参与肺鳞癌发生发展的过程[7]。然而，lncRNA m6A修饰在肺鳞癌预后及免疫功能中的作用多为单个lncRNA的功能研究，不具有独立的预测价值。因此，探究联合lcRNA m6A修饰在预测肺鳞癌预后价值及免疫疗效价值中的意义重大。

1 资料与方法

1.1 数据提取与整理

从TCGA肺鳞癌数据库中下载转录组数据、体细胞突变及临床信息，得到501个肺鳞癌组织和49个非肿瘤组织的样本数据，剔除6例患者生存时间不详的肿瘤组织样本。利用Perl编程语言对肺鳞癌组织的原始数据进行整理和注释，得到14 086个LncRNA的表达矩阵。从TCGA数据库中提取已被广泛报道的23个m6A调控因子(甲基转移酶METTL3、METTL14、METTL16、WTAP、VIRMA、ZC3H13、RBM15、RBM15B；阅读蛋白YTHDC1、YTHDC2、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、LRPPRC、HNRNPA2B1、RBMX、IGFBP1、IGFBP2、IGFBP3；去甲基转移酶FTO、ALKBH5)的表达矩阵，通过R语言Rx64 4.1.2 “limma”软件包对lncRNA和m6A调控因子的表达矩阵整合鉴定出1 136个lncRNA m6A修饰，在此基础上提取肺鳞癌患者样本的生存时间和生存状态，最终得到27个具有肺鳞癌预后意义的lncRNA m6A修饰。

1.2 预后风险模型的构建及评估

基于上述27个lncRNA m6A修饰，利用R语言Rx64 4.1.2 “Survival”“caret”“glmnet”“survminer”软件包构建Lasso回归分析的内置函数cvfit和lambda，根据该函数中的“λ”值对构建预后风险模型的lncRNA m6A修饰个数进行取值，“λ”值代表该模型的过拟合程度，当“λ”值为13时模型的拟合程度最佳，对应地鉴定出13个lncRNA m6A修饰以构建预后风险模型，即LINC02332、AL591686.1、DHRS4-AS1、AC008734.1、AC138035.1、AL390195.3、AP001189.3、WT1-AS、AL391095.4、AC130651.1、AC015922.2、HORMAD2-AS1、AC243919.2。随后，利用风险评分公式对上述13个lncRNA m6A修饰的相关回归系数进行加权计算。为验证该预后风险模型对肺鳞癌预后的预测能力，首先根据5∶5留出法对该预后风险模型中的肺鳞癌样本进行随机抽样得到训练集和测试集，训练集和测试集的患者样本在年龄、性别、TMN分期等变量上基本无差异。之后将训练集的患者样本根据风险评分的中位值划分为高、低风险组，并以训练集的风险评分中位值划分测试集的高、低风险组。本研究利用“tidyverse” “ggplot2” “ggExtra”软件包分析13个LncRNA m6A修饰与纳入的23个m6A调控因子之间的相关性；“survival”“timeROC”“survminer”软件包比较高、低风险组间的总生存期(OS)、ROC曲线下面积(AUC)差异及独立风险因素；“org.Hs.eg.db”“clusterProfiler”“enrichplot”等软件包对高、低风险组间差异基因进行GO富集分析；“GSVA”“GSEABase”“reshape2”“pheatmap”软件包进行高低风险组差异基因免疫功能分析。

1.3 免疫疗效评估

从TCGA数据库中下载肺鳞癌的体细胞突变数据，并应用perl编程语言和“maftools”软件包进行数据提取和可视化整理，得到高、低风险组每个样本的肿瘤突变负荷(tumor mutation burden,TMB)，即每100万个碱基中发生突变的碱基数，通过“limma” “ggplot2”软件包分析高、低风险组间的TMB差异。从肿瘤免疫功能障碍和排斥(tumor immune dysfunction and exclusion,TIDE)数据库中下载NSCLC患者的TIDE评分，利用“ggpubr” “limma”软件包分析高、低风险组间的TIDE差异。通过CIBERSORT算法计算高、低风险组样本的免疫细胞评分、基质细胞评分及TME总评分，分析高、低风险组间的TME差异。

1.4 统计学方法

采用Perl编程语言和R语言Rx64 4.1.2进行数据处理，采用SPASS 22.0软件进行数据统计分析。单因素Cox回归分析取95%置信区间(CI)。采用双侧Log-rank检验分析Kaplan-Meier(K-M)曲线，Wilcoxon检验分析两组变量之间的差异。检验水准α除Pearson分析设定为0.001外，其余统计学分析均设定为0.05。

2 结果

2.1 构建预后风险模型

基于TCGA肺鳞癌数据库，利用相关性分析得到1 136个lncRNA m6A修饰(图1A)。利用单因素Cox回归分析得到具有预后意义的27个lncRNA m6A修饰(图1B)，其中AL591686.1、AC008734.1、AP001189.3、WT1-AS、AL391095.4、AC015922.2、HORMAD2-AS1、DSCR9与肺鳞癌预后呈负相关，其他19个lncRNA m6A修饰呈正相关。利用Lasso回归分析cvfit和lambda函数，鉴定出13个联合lncRNA m6A修饰预后标志物用以构建预后风险模型，即LINC02332、AL591686.1、DHRS4-AS1、AC008734.1、AC138035.1、AL390195.3、AP001189.3、WT1-AS、AL391095.4、AC130651.1、AC015922.2、HORMAD2-AS1、AC243919.2(图1C～D)。上述13个LncRNA m6A修饰的表达与不同的m6A调控因子的表达相关，其中AC138035.1与HNRNPA2B1、HORMAD2-AS1与IGFBP1、AC243919.2与YTHDC1的表达之间存在较强的正相关，AP001189.3与HNRNPC的表达之间存在较强的负相关(图1E)。

图1 预后风险模型的构建

2.2 预后风险模型的评估

训练集和测试集基本互斥且不存在分组差异(表1)。随着样本风险评分增加，患者的生存率下降(图2A～B)。K-M曲线分析表明,在训练集和测试集中均显示高风险组患者的OS低于低风险组(图2C～D)。同时，单因素Cox回归分析显示,相较于年龄、性别、临床分期等临床特征，风险评分是肺鳞癌相关的OS预后指标(图3A～B);多因素Cox回归分析表明,风险评分作为独立的预后因子，预测价值高(P<0.001)(图3C～D)。

表1 训练集与测试集TCGA数据库肺鳞癌患者样本的临床信息及分组情况[n(%)]

图2 预后风险模型的风险趋势与K-M曲线

图3 单因素、多因素Cox回归分析预后风险因素

2.3 GO通路富集分析及免疫功能差异分析

GO通路富集分析结果显示，上述67个差异基因主要执行信号受体激活因子活性、受体配体活性等功能，并参与ERK1/2级联调控等生物学过程(图4)，高、低风险组的差异基因与趋化因子受体(CCR)、人类白细胞抗原(HLA)、干扰素(IFN)Ⅰ/Ⅱ型反应、MHC Ⅰ类分子、细胞活性、炎症刺激、T细胞共抑制/共刺激分子、免疫检查点等13种免疫分子及功能密切相关(图5)。

图4 高、低风险组差异基因的GO通路富集分析 图5 高、低风险组免疫相关功能差异分析

2.4 免疫治疗疗效分析

目前对于免疫治疗疗效的分析主要包括TMB、TIDE及TME评分等方法。TMB能够衡量体细胞的非同义突变，即产生新抗原或肽段被自身免疫系统识别为非自身抗原的概率[8]。TMB瀑布图显示,TP53在高、低风险组中突变频率最高(图6A～B)，高、低风险组的TMB并无明显差异(图6C)。TIDE评分表明,高风险组的TIDE评分高于低风险组(图6D)。此外，TME评分表明,高风险组的免疫细胞评分(图6E)、基质细胞评分(图6F)以及TME总评分(图6G)均高于低风险组。

图6 免疫疗效分析

3 讨论

相较于其他类型的非小细胞肺癌，肺鳞癌体细胞突变频率更高、免疫原性更强，患者更适合应用ICI免疫疗法改善预后[9]，因此探索能够评估肺鳞癌患者的预后及ICI疗效的生物标志物不可或缺。研究[10]发现,AL049840.3、AC008770.3、AL355312.3等11个lncRNA m6A修饰与胃癌的不良预后密切相关，通过聚类分析发现,C1和C2亚组的胃癌患者免疫治疗的应用前景更好。在膀胱癌患者样本的分析中，研究[11]发现,PTOV1-AS2、EHMT2-AS1、KCNQ1OT1等lncRNA m6A修饰可能是参与膀胱癌免疫微环境的重要组成。本研究涉及的LncRNA m6A修饰中，DHRS4-AS1被报道通过拮抗TP53和TET1调节肿瘤干细胞在NSCLC中的干性，从而发挥抑癌作用[12]。研究[13]发现，AP001189.3与铁蛋白代谢相关，并能够影响结肠腺癌患者的结局。WT1-AS广泛表达于实体肿瘤和血液肿瘤中，在胶质瘤中通过调节miR-494-3p抑制恶化[14]，在非霍奇金淋巴瘤的患者中检测到高表达WT1的特异性T细胞[15]，同时研究[16]发现，WT1能够作为胸腺上皮肿瘤的免疫治疗靶点，且能作为免疫抑制细胞因子IL-10(IL-10)新的转录激活因子，调控IL-10基因表达，从而发挥免疫调节作用[17]。因此，本研究基于上述lncRNA m6A修饰所建立的预后风险模型可信度高，以预后模型的回归系数为基础得到风险评分不仅能够作为独立的预后因子，同时能够较准确地预测患者的1、3、5年生存情况。上述结果表明，本研究所构建的预后风险模型在肺鳞癌预后的预测方面是高度可信且稳定的。

ICI免疫治疗作为肺鳞癌患者改善预后的关键突破口，其应答效率可通过TMB、TIDE和TME评分进行评估[18]。在一项晚期NSCLC的研究[19]中发现,高TMB患者ICI的有效作用时间为64周，比低TMB患者更长。而另一项研究[20]中，在阿特珠单抗(PD-L1抑制剂)单药治疗时，评估晚期NSCLC患者治疗前的TMB发现，治疗后高TMB患者的无进展生存期(PFS)得到明显改善。本研究中，低风险组和高风险组整体突变频率无明显差异。有研究[21]通过对TGFB1和SOX10等影响TIDE的基因进行综合分析，在转移性黑色素瘤中发现，TIDE评分较高的患者OS更高。在肝癌的疗效评估中，研究[22]发现，TIDE评分较高的高危组(113/186)比低危组(89/186)免疫治疗应答患者人数更多。类似地，本研究中高风险组的TIDE得分明显高于低风险组，初步说明ICI在高风险组中疗效更佳。此外，本研究发现，高风险组的患者TME评分更高，即高风险组患者的ICI应用更有效。上述结果表明，本研究所构建的预后风险模型预测肺鳞癌的免疫疗效可靠性强。

然而，本研究仍存在一些不足。首先，用于构建肺鳞癌患者预后风险模型的样本量相对不足，因其仅包含了TCGA数据库。其次，研究结果仅依赖于生物信息学分析，后期仍需通过大量基础实验和临床试验验证13个lncRNA m6A修饰在肺鳞癌中的作用，同时对接受免疫治疗和化疗肺鳞癌患者进行多中心大规模的验证性临床试验。

本研究鉴定出LINC02332、AL591686.1、DHRS4-AS1、AC008734.1、AC138035.1、AL390195.3、AP001189.3、WT1-AS、AL391095.4、AC130651.1、AC015922.2、HORMAD2-AS1、AC243919.2共13个lncRNA m6A修饰，并以此构建肺鳞癌的预后风险模型，K-M曲线显示,高风险组非鳞癌患者OS明显低于低风险组；AUC显示,此预后风险模型用于预测肺鳞癌预后情况可信度高。随后,通过分析TIDE、TME证实了此预后风险模型能够有效预测ICI疗效。

综上所述，本研究鉴定的13个lncRNA m6A修饰联合应用对评估肺鳞癌预后及ICI疗效具有重要价值，为肺鳞癌ICI疗法的应用及敏感性化疗药物的开发提供了重要的理论依据。