m6A相关lncRNA 预测肺鳞癌患者预后和免疫反应分析①

毛宇泽,杨成鹏,姜腾蛟,骆 鹏,朱晓峰

(佳木斯大学附属第一医院心胸外科,黑龙江佳木斯 154003)

肺癌是世界上最常见的恶性肿瘤,非小细胞肺癌(NSCLC)约占所有肺癌的85%,包括腺癌、鳞癌和大细胞肺癌等[1]。肺鳞状细胞癌(LUSC)是一种起源于支气管鳞状上皮且与吸烟密切相关的恶性肿瘤,约占NSCLC病例的20%～30% 。与DNA和蛋白质相类似,RNA经历转录后修饰,其中包括N1-甲基化(m1A)、5-甲基化(m5C)、N6-甲基化(m6A)和7-甲基化(m7G)。m6A修饰是真核生物中调控基因表达的最丰富的RNA修饰。这些m6A修饰包括甲基转移酶、甲基结合蛋白和去甲基化酶,分别被命名为“编写器”、“阅读器”和“橡皮擦”[2～4]。研究表明,m6A调节因子通过调节肿瘤的生长和进展以及治疗反应,在NSCLC中发挥潜在作用[5]。m6A阅读器YTH n6-甲基化RNA结合蛋白1(YTHDF1)抑制NSCLC中的肿瘤生长和集落形成,而YTHDF1的低表达水平与临床结果差相关,并通过上调抗氧化系统促进癌细胞对顺铂的耐药性[6]。因此,m6A的修饰在NSCLC中起着至关重要的作用。

本研究在LUSC中筛选了m6A相关的lncRNA建立并验证了预后特征。进一步区分了高、低危组,并分析了免疫细胞在LUSC中的分布及免疫疗效。有助于预测LUSC的预后并区分LUSC患者是否受益于免疫治疗。

1 资料与方法

1.1 获取LUSC RNA-seq转录组及临床数据

从TCGA数据库(https://cancergenome.nih.gov/)和TCGA-LUSC(n=551)下载LUSC的标准化RNA-seq数据和临床信息(年龄、性别、TNM分期、生存时间和状态),其中包括501个肿瘤样本和49个正常样本。我们在分析中排除了缺失生存信息的样本。

1.2 m6A相关lncRNA的鉴定

从TCGA-LUSC中提取30个m6A调控因子的RNA表达矩阵,包括“编写器”(METTL3,METTL14,METTL16,METTL5,WTAP,VIRMA,RBM15,RBM15B,ZC3H13,CBLL1,和ZCCHC4),“阅读器”(YTHDF1,YTHDF2,YTHDF3,YTHDC1,YTHDC2,HNRNPA2B1,HNRNPC,FMF31,IGF2BP1,IGF2BP2,IGF2BP3EPLAVL1,G3BP1G3BP2,PRRC2A,RBMX),和“橡皮擦”(FTO和ALKBH5)。此外,采用Pearson相关分析对m6A相关的lncRNA进行筛选(|R|>0.4,P<0.001)。

1.3 预后m6A相关lncRNA特征的建立和验证

采用单因素Cox回归分析筛选预后后与m6A相关的lncRNA(P<0.05)。采用R包“glmnet”进行10倍交叉验证,以进一步缩小范围。最后,利用λmin(λmin = -4.5)鉴定了7个lncRNA并建立了与m6A相关的lncRNA预后特征,使用以下公式计算了风险评分:Risk score = coef (lncRNA1) × expr (lncRNA1) + coef (lncRNA2) ×expr (lncRNA2) +……+ coef (lncRNA) × expr (lncRNAn), 其中coef为系数,expr(lncRNAn)为lncRNAs的表达量。

将数据集以7:3的比例被划分为训练集和测试集。根据每个患者的风险评分,使用R软件进行生存分析、风险曲线和绘制热图。使用ROC评估了模型对预测预后的敏感性。AUC使用R包“cesivivalROC”计算。采用单因素和多因素Cox回归分析来评估该模型作为一个独立的预后因素的潜力。采用分层分析来评估该模型对不同临床病理特征的预测能力。

1.4 功能分析

使用GSEA软件揭示高危组和低危组之间差异表达的功能通路。从Molecular Signatures Database数据库下载并分析两个广泛应用的基因集(h.all.v7.2.symbols.gmt [cancer hallmarks]和c7.all. v7.4.symbols.gmt [immunologic signatures])。为了获得标准化富集分数,进行1000次基因集排列。其中具有|NES|>1, nominalP<0.05,and q<0.25的通路被认为是显著富集的。

1.5 肿瘤突变负荷(TMB)和肿瘤浸润性免疫细胞(TICs)分析

采用CIBERSORT算法(http://cibersort.stanford.edu)在LUSC样本中获得了22个不同的TICs。采用TIDE(http://tide.dfci.harvard.edu/)来预测各亚组中免疫治疗的疗效。此外,从TCGA数据库中下载了LUSC患者资料,并使用“maftools”软件包进行分析和可视化。最后,计算每个患者的TMB。

1.6 统计学方法

采用R版本4.1.0和Perl编程语言(版本5.34.0)进行。基于预后m6A相关lncRNA的表达情况,采用K-M生存曲线分析和log-rank检验来比较不同亚组间的OS。此外,进行单因素和多因素Cox回归分析来评估特征在预测OS方面的独立预后价值。采用Wilcoxon检验来确定两组间变量的差异。P<0.05为差异有统计学意义。

2 结果

2.1 LUSC中m6A相关lncRNA的鉴定

为识别LUSC中m6A相关的lncRNA,从癌症基因组图谱(TCGA)-LUSC数据集中提取30个m6A调控因子和14086个lncRNA表达矩阵进行皮尔逊相关分析,筛选出373个m6A相关的lncRNA(|R|>0.4,P<0.001)。与m6A相关的lncRNA共表达网络见图1a。为了筛选LUSC中重要的lncRNA,将临床预后和lncRNA的表达整合到模型中。使用单变量Cox回归分析筛选了9个与预后相关的lncRNA,见图1b。其中 AP001189.3和HORMAD2-AS1(HR>1、P<0.05)与预后不良相关,而AL122125.1、AP001469.3、AC138035.1、AC002、AC002467.1、AC020658.4、AC243919.2、PRC1-AS1与预后良好相关(HR <1、P<0.05)。

图1 m6A相关LncRNAs和m6A相关预后LncRNAs的鉴定a.共表达网络; b.森林图。

2.2 m6A相关lncRNA特征的建立和验证

基于上述与m6A相关的lncRNA,完成LUSC的lncRNA特征的建立,对9个预测的m6A相关的lncRNA进行了Lasso回归分析。确认了AL122125.1、AP001469.1、AC138035.1、AC002467.1、AC243919.1、AC243919.2、AP001189.3、HORMAD1-As1构建风险模型。为了探究lncRNA签名的作用,将整个集合(n=495)分为训练集(n=347)和测试集(N = 148),并对这些集合进行了Kaplan-Meier生存分析。结果显示,高危组的预后比低危组差(训练组和测试组的P=0.003和P=0.002;图2a,b)。此外,两组患者的风险曲线显示,随着风险评分的增加,LUSC的预后进一步恶化(图2c,d)。

构建一个受试者工作特征曲线(ROC)来评估我们的模型的敏感性。训练组1年、3年和5年总生存率(OS)的曲线下面积(AUC)分别为0.581、0.581和0.608。在测试集中,1-、3-和5年OS的AUC分别为0.655、0.655和0.555。这些结果表明,已建立的LUSC的预后m6A相关特征具有一般的敏感性(图2e,f)。K-M生存曲线显示,AP001189.3和HORMAD2-AS1表达水平较高与预后不良相关,AC002467.1、AC243919.2、AC138035.1、AL122125.1和AP001469.3表达水平较低,与LUSC良好预后相关(图3a～g)。

图2 m6A相关的lncRNA的预后特征a-b.训练集和测试集中的Kaplan-Meier生存曲线;c-d.训练集和测试集中,LUSC患者的风险评分和生存状态的分布; e-f .训练集和测试集中,与m6A相关的预后LncRNAs特征的ROC曲线。

图3 Kaplan-Meier生存曲线a.AP001189.3; b.HORMAD2-AS1; c.AC002467.1; d.AC243919.2; e.AC138035.1; f.AL122125.1; g.AP001469.3(P<0.05)。

2.3 m6A相关LncRNA特征的分层和独立预后分析

为了更好地评估已建立的m6A相关预后特征是否可以作为一个独立的预后因素,基于临床病理特征对训练集和测试集进行了单变量和多变量Cox回归分析。基于预后m6A相关lncRNA特征的风险评分可以预测训练和测试集中LUSC的预后(单变量: HR=3.437、P<0.001和HR=6.960、P=0.005;多变量: HR=3.475、P<0.001和HR=10.386、P=0.003;图4a～d)。K～M生存分析结果显示,低危组在年龄(≤65岁、>65岁)、性别(男、女)、TNM期(I～II期、III～IV期)、T期(T1～2)、N期(N0、N1～3)、M期(M0)(P<0.05;图4e-l)方面均有更好的预后。这些结果表明,高危评分的模型与预后不良相关,可作为LUSC的独立预后因素。

图4 m6A相关预后特征的独立预后和分层分析a-b.训练集中m6A相关预后特征的单变量和多变量回归分析;c-d.测试集中m6A相关预后特征的单变量和多变量回归分析;e-l .高危组和低危组的Kaplan-Meier生存曲线(P-value<0.05)。

图5 高低风险组的功能富集分析a.在低风险组中具有丰富的肿瘤特征:细胞周期;b-d.高危组中富集的肿瘤特征: MAPK信号通路、细胞凋亡、JAK/STAT信号通路;e-h.免疫相关信号通路富集(e:低危组,f-h:高危组)。

2.4 m6A预后相关lncRNA特征的功能分析

为了检验高风险组和低风险组之间的功能差异,我们进行了基因集富集分析(GSEA),发现有几种不同的肿瘤特征在两个风险组之间有显著差异富集。细胞周期通路在低危组中富集(图6a),而MAPK、凋亡和JAK/STAT信号通路在高危组中富集(图6b～d)。值得注意的是,两个危险组之间许多差异表达的免疫相关信号通路与T细胞、B细胞和单核细胞相关(图6e中的低风险组,图6f～h中的高危组)。因此,研究结果表明,预后后与m6A相关的lncRNA与肿瘤发生和免疫通路相关。

2.5 肿瘤突变负荷(TMB)和肿瘤浸润性免疫细胞(TICs)分析结果

我们分析了高低风险组中每个LUSC患者的突变谱。在总数中,显著突变的前20个基因为TP53、TTN、TTN、CSMD3、MUC16、MUC16、RYR2、LRP1B、USH2A、SYNE1、ZFHX4、FAM135B、KMT2D、XIRP2、SPTA1、CDH10、NAV3、PCDH15、PAPPA2、RYR3、DNAH5、PKHD1(图6a低危组,图6b高危组)。低危组的TMB显著高于高危组(图6c),且风险评分与TMB呈负相关(图6d;R=-0.13,P<0.05),提示低危组可能表现出更好的免疫应答。

图6 TMBa.低风险组突变基因;b.高危组突变基因;c.低危组的TMB高于高危组(P<0.05);d.风险评分与TMB呈负相关(R=-0.13,P=0.0032)。

比较两组之间的免疫细胞的差异,发现低风险组比高危组中幼稚B细胞,滤泡辅助T细胞,和激活NK细胞更丰富,而CD4记忆休息T细胞,M2巨噬细胞和中性粒细胞相对较少(图7)。然后我们进行了肿瘤免疫功能障碍和排斥(TIDE)分析,探讨两个亚组之间的差异,发现高危组的TIDE评分高于低危组,提示高危组的肿瘤免疫逃逸和免疫治疗反应较低(图7b)。此外,我们还进行了免疫细胞与风险评分之间的相关性分析。我们发现M2巨噬细胞、中性粒细胞和CD4记忆静息T细胞与风险评分呈正相关,而活化的NK细胞、滤泡辅助T细胞和幼稚B细胞与风险评分呈负相关(P<0.05)。这些结果表明,与m6A相关的预后lncRNA与肿瘤免疫相关,可能成为新的免疫治疗的生物标志物。

图7 高、低危组的免疫学特征a.TICs差异;b.箱线图。

3 讨论

肺癌是一种高度异质性的疾病,由于不同的组织学和分子特征,不同亚型的肺癌的治疗效果也有所不同[7]。晚期肺鳞癌尽管其治疗方法有限,但免疫治疗已成为LUSC的标准治疗方法。然而,只有一小部分表达高水平免疫检查点如PD-L1的LUSC患者从免疫治疗中获益[8]。因此,探索其他生物标志物来识别可能受益于免疫治疗的LUSC患者是至关重要的。

N6-甲基化(m6A)是一种在真核生物中调控基因表达RNA修饰[9]。有研究表明,m6A调控因子的异常表达通过在肿瘤生长、增殖、侵袭和转移中发挥作用,促进肿瘤的发展和进展[10]。m6A相关基因在肿瘤中起着关键作用,但在LUSC中的功能与作用知之甚少。

在本研究中,建立并验证了一个基于LUSC中与m6A相关的7个lncRNA的预后特征。在高低风险组中确定了特异性的癌症和免疫相关通路,并研究了预后m6A相关lncRNA与肿瘤免疫之间的关系。M1巨噬细胞具有促炎和吞噬肿瘤细胞的作用,而M2巨噬细胞具有抗炎作用并促进肿瘤增殖[11]。结果显示,高危组中M2巨噬细胞的数量高于低危组,提示高危组患者有更高的肿瘤进展风险。此外,TMB被广泛称为基于免疫检查点抑制剂(ICI)治疗应答的替代或补充生物标志物[12],高TMB患者可能比低TMB患者从ICI治疗中获益更多[13]。在我们的研究中,低危组的TMB评分高于高危组,TIDE评分较低,高危组相反,说明低危组对免疫治疗的反应可能高于高危组。

综上,本研究是对LUSC中预后m6A相关lncRNA的综合分析,建立并验证了一个与m6A相关的lncRNA预后特征,可用于区分受益于免疫治疗的LUSC患者。本研究也加深了对m6A修饰在肿瘤中重要作用的理解,并希望对未来研究有指导作用。