基于生物信息学筛选影响胃癌患者预后的糖酵解相关基因

赵旭辉,黄小敏,达德转,许 焱,崔晓东,李红玲

(1.甘肃中医药大学第一临床医学院,甘肃 兰州 730000;2.甘肃省人民医院 肿瘤内科,甘肃 兰州 730000)

胃癌(gasrtic cancer, GC)是消化系统最常见的恶性肿瘤之一,为世界第四大常见癌症,居癌症相关死亡原因的第2位[1]。GC发病隐匿导致诊断时间较晚且预后不良,平均5年生存率低于20%[2]。目前手术、放疗、化疗临床综合治疗效果不理想,含铂类药物的一线化疗治疗晚期GC已达到疗效瓶颈[3-4]。虽然分子靶向治疗以及PD-1/PD-L1免疫检查点抑制剂对晚期GC有一定的疗效,开辟了新的治疗途径。然而,一些GC患者的疗效仍然较差[5-6]。20世纪初德国著名的生物化学家Otto Warburg首次从细胞代谢角度提出了关于肿瘤的“Warburg效应”,即在肿瘤细胞中当线粒体功能失调后,主要通过增强有氧酵解的方式来获取充足的能量[7]。越来越多的证据表明,异常的能量代谢方式已成为肿瘤细胞生长的特殊模式,葡萄糖的异常消耗与乳酸的大量堆积使得肿瘤细胞更具快速生长、转移的优势[8]。同时此过程也会使细胞外pH值明显下降造成微环境的酸化,微环境的酸化可破坏邻近的正常细胞且造成细胞外基质降解,从而促进肿瘤细胞的入侵[9]。有学者指出这种特殊的产能方式以及微环境的变化特征有利于肿瘤细胞逃避机体的凋亡机制,导致患者不良预后[10]。有文献报道,糖酵解相关基因可预测癌症患者的预后,包括透明细胞肾细胞癌[11]、肺腺癌[12]、肝细胞癌[13]、乳腺癌[14-15]、卵巢癌[16]和结直肠癌[17]。研究表明,糖酵解相关基因可用于有效评估GC患者预后[18-20]。然而,目前仍然缺乏关于糖酵解相关基因与GC患者预后关系的系统研究。基于此,本文通过利用生物信息学开发糖酵解相关基因以预测GC患者预后,为临床提供参考。

1 资料与方法

1.1患者的临床信息和信使核糖核酸(messenger RNA,mRNA)表达数据集 我们从癌症基因组图谱(the cancer genome atlas, TCGA)数据库(https://cancergenome.nih.gov/)中提取了GC患者的临床资料和mRNA表达谱。纳入445例患者的临床信息,并统计了匹配的年龄,分期,分级,将样本的基因转录数据统一转化为log2(FPKM 值+1)后进行数据集序列整合。

1.2方法

1.2.1数据处理和风险评分计算 在分子特征数据库v4.3.2检索所有糖酵解相关基因集,包括5个不同的基因集(BIOCARTA_GLYCOLYSIS_PATHWAY、GO_GLYCOLYTIC_PROCESS、HALLMARK_GLYCOLYSIS、KEGG_ GLYCOLYSIS_GLUCONEOGENESIS、 REACTOME _GLYCOLYSIS),查阅相关文献后共筛选出371个糖酵解相关基因。接下来,我们采用单变量Cox回归分析来确定与总生存期(overall survival,OS)相关基因。我们应用最小绝对收缩和选择算子(least absolute shrinkage and selection operator, LASSO)回归分析来确定具有非零系数的潜在预测因子。然后,进行多变量Cox回归分析以确定最终建模中涉及的枢纽基因并获得相关数。

1.2.2预测模型的构建与验证 以中位风险评分为临界值,将445例GC患者分为两个亚组(高风险组和低风险组),并绘制Kaplan-Meier生存分析图以及无进展生存期(progression-free survival,PFS)曲线,以比较患者生存率与log-rank检验的差异。此外,在R软件(版本 4.0.2)中生成1年、3年和5年的受试者工作特征(receiver operating characteristic, ROC)曲线,以及临床相关特征的ROC 曲线,以计算每个预测因子模型的曲线下面积(area under the curve, AUC),以进一步评估模型的效率和准确性。P<0.05为差异有统计学意义。

1.2.3列线图的构建和评估 通过R软件(版本 4.0.2)构建列线图,并采用列线图模型比较GC患者的生存概率。生成校准图和C指数作为列线图性能的评估。

1.2.4高风险组和低风险组间的基因集变异分析(gene set variation analysis,GSVA) 为了进一步探索模型中高风险组和低风险组间生物途径激活状态的差异,使用“GSVA”R包进行GSVA。从MSigDB数据库下载“c2.cp.kegg.v7.4.symbols”基因集用于GSVA。GSVA通常用于以非参数和无监督方法估计表达数据集样品中途径和生物过程活性的变化。

1.2.5免疫细胞浸润分析 通过R语言对两组样本中免疫细胞的突变、免疫细胞浸润、免疫评分等差异进行了分析。使用“maftools”R包评估GC样本中22个免疫细胞的突变。GC样本中22个免疫细胞的浸润水平基于CIBERSORT(http://cibersort.stanford.edu/)进行评估。使用R包评估GC样本中22个免疫细胞的免疫和基质评分。

1.2.6统计学方法 使用Excel软件和R软件进行统计分析。通过单变量和多变量Cox回归分析评估单个指标的预后意义。采用Kaplan-Meier 生存曲线和log-rank检验评估预后结局。采用非配对t检验比较两组不同表达水平。R软件用于绘制热图,ROC曲线,富集,森林和校准图。P<0.05表示差异具有统计学意义。

2 结 果

2.1糖酵解相关基因预后特征的构建和验证 使用从TCGA获得的445个GC样本来构建模型并采用log-rank检验方法和利用“survival”R 包对样本基因进行单变量Cox回归分析探究糖酵解相关差异基因与患者OS的关系。单因素Cox回归分析显示,19个糖酵解相关差异基因与患者的预后相关。接下来,选择上述的19个糖酵解相关基因进行LASSO回归分析,得到15个糖酵解相关基因(见图1及表1)。基于LASSO回归分析的结果,使用15个糖酵解相关基因(PFKFB2、UHRF1、ACYP1、CLDN9、STC1、EFNA3、NUP50、、ADH4、ANGPTL4、PKP2、VCAN、HIF1A、LHX9、ANKZF1、ALDH3A2)建立和验证预测GC患者结局的风险模型。

图1 糖酵解相关差异基因与患者OS的单变量Cox回归分析

表1 通过LASSO回归分析得到与GC患者预后相关的15个糖酵解相关基因

采用log-rank检验方法,根据风险评分的中值将TCGA中的GC患者分为高风险组和低风险组,两组OS和FPS差异有统计学意义(均P<0.01)。ROC曲线表明,1年、3年、5年OS的AUC为分别为0.654、0.729、0.799,表明上述15个糖酵解相关基因标记在评估GC患者预后结局方面具有良好的特异度和敏感度。见图2～4。

图2 用于构建最终预测模型的 LASSO 回归分析 a.选择LASSO模型中的最优参数(λ);使用最小标准在最佳值绘制垂直虚线; b. LASSO 分布的部分似然偏差

图3 高风险组和低风险组的OS和PFS

图4 1年、3年和5年OS的ROC曲线

2.2整合糖酵解相关基因标记的列线图模型构建 构建基于基因特征的风险评分与临床病理特征(分期、性别、分级、风险和年龄)相结合的列线图模型,用以评估GC患者的生存可能性(图5a)。列线图证明了这15个基因标记在个体化生存估计中的价值。在校准、辨别和临床可用性方面评估列线图的性能,C指数为0.8557(图5b)。ROC曲线中风险(AUC=0.775)和列线图(AUC=0.753)具有较好的预测价值(图5c)。基于列线图单因素和多因素分析,发现在单因素分析中年龄[HR=1.024,95%置信区间(CI):1.006～1.042,P<0.05]、分级[HR=1.549,95%CI:1.244～1.929,P<0.01]和列线图[HR=1.049,95%CI:1.036～1.061,P<0.01]与预后结局相关。多因素Cox分析表明,列线图[HR=1.039,95%CI:1.020～1.059,P<0.01]可作为独立的预后因子。见图5。

图5 基于糖酵解和临床因素的列线图 a.包含年龄和风险评分的列线图,以预测1年、3年和5年的OS; b.列线图的校准图; c.风险、列线图和其他临床特征的ROC曲线

2.3基于15个基因特征的风险评分作为独立的预后因素 采用单因素和多因素分析,分析各临床病理特征的预后价值。单因素分析结果表明,年龄[HR=1.026,95%置信区间(CI):1.008～1.044,P<0.01],分级[HR=1.534,95%CI:1.241～1.896,P<0.01]和风险系数[HR=3.690,95%CI:2.468～5.518,P<0.01]与预后结局相关(见图6a)。随后,多因素Cox分析表明,年龄[HR=1.038,95%置信区间(CI):1.019～1.057,P<0.01],分级[HR=1.549,95%CI:1.238～1.937,P<0.01]和风险系数[HR=3.686,95%CI:2.390～5.686,P<0.01]是独立的预后指标(见图6b)。同样,通过绘制ROC曲线发现年龄(AUC=0.606),分级(AUC=0.606), 风险系数(0.799)具有很好的预测价值(见图6c)。

图6 临床特征的单因素、多因素分析及ROC曲线 a.单变量分析; b.多变量分析; c.ROC曲线

2.4基于糖酵解相关基因的临床相关性分析 临床相关性分析发现,高低风险组年龄和性别差异无统计学意义。G2与G3患者分级差异有统计学意义。在TNM分期中,T1期患者与T2、T3、T4期患者差异有统计学意义(P<0.01),T1期之后的患者风险明显上升。 M0与M1期患者风险差异具有统计学意义,N0期与N1、N3期患者风险差异具有统计学意义,而N0与N2期患者风险差异无统计学意义。见图7。

图7 基于糖酵解相关基因的临床相关性分析 a.年龄; b.性别; c.分级; d.分期; e-g.TNM分期

2.5糖酵解相关基因的表达水平与免疫细胞浸润的分析 不同水平的免疫浸润可能在一定程度上解释了为什么具有相同组织学类型的癌症患者可能具有不同的临床结果[20]。我们分析糖酵解相关基因与免疫浸润水平的相关性和免疫相关功能的相关性,旨在揭示15个糖酵解相关基因标记影响GC患者预后的潜在机制。箱线图分析显示,高风险组细胞浸润水平(巨噬细胞 M1)高于低风险组(P<0.05),而高风险组激活的CD4+T细胞、辅助性滤泡T细胞、单核细胞(P<0.05)水平低于低风险组(见图8a)。关于免疫相关功能的分析显示,高风险组CCR、Parainflammation、Type_I_IFN_Reponse、Type_II_ IFN_Reponse功能高于低风险组(均P<0.05),而高风险组MHC_class_I功能低于低风险组(见图8b)。

图8 箱线图分析 a.两组免疫细胞浸润分布; b.免疫相关功能

2.6GSVA分析确定的高低风险亚组差异生物学行为 进一步对上述糖酵解相关基因进行GSVA分析,以探索高低风险组间生物学通路激活状态的差异。利用热图对相关生物化过程进行可视化分析并根据P值大小进行排列。结果表明,高风险组致癌激活途径显著富集,ECM_RECEPTOR_INTERACTION、FOCAL_ADHESION、CHEMOKINE_SIGNALING_PATHWAY、TGF_BETA_SIGNALING_PATHWAY等。这些通路可以调控肿瘤发生发展,因而可能为高危人群预后不良提供一定的理论依据。见图9。

图9 热图 注:红色表示激活途径,蓝色表示抑制途径

3 讨 论

近年来,关于代谢组学的研究越来越受到关注,人类恶性肿瘤与能量代谢有着密不可分的联系[21]。在氧气充足的条件下,肿瘤细胞也会通过无氧糖酵解的方式给肿瘤细胞提供更多的能量,使肿瘤细胞得以快速生长(即“Warburg效应”)[22]。关于肿瘤细胞与糖酵解的研究近年来层出不穷。Hu等[23]发现,在胰腺癌中,UHRF1可通过抑制SIRT4促进无氧糖酵和肿瘤增殖。Zhao等[24]发现LncRNA MIR17HG可以通过介导糖酵解相关正反馈回路促进结直肠癌肝转移。还有研究发现,EB病毒编码的miRNAs也可以调节细胞的能量代谢,增加细胞摄取葡萄糖的能力并促进细胞内糖酵解速率。Cai等[25]还发现EBV-miR-BART1可依赖糖酵解并通过PTEN影响 HIF1α从而诱导鼻咽癌细胞的转移与扩散。然而,年龄和转移诊断等临床病理学特征不足以准确预测癌症患者预后。因此,越来越多的mRNA被确定为肿瘤进展或预后的生物标志物,并且已经评估了生物标志物的临床意义[26-27]。本研究中,在TCGA中获得了445例GC患者的临床资料和mRNA表达谱,然后在分子特征数据库v4.0上检索了5个不同糖酵解相关基因集,查阅相关文献后共筛选出371个糖酵解相关基因。接下来,我们采用单变量Cox回归分析来得到与预后相关的19个糖酵解相关的差异基因,接着对这些基因进行LASSO回归分析,基于LASSO回归分析的结果,使用15个糖酵解相关基因建立和验证预测GC患者结局和系数的风险模型。

其次,采用log-rank检验方法,根据风险评分的中值将TCGA-STAD中的GC患者分为高风险组和低风险组,进一步通过 Kaplan-Meier 曲线、ROC 曲线、单因素及多因素 Cox 回归和通过绘制列线图分析考察风险评分模型预测能力,结果显示高风险组病死率高于低风险组,ROC曲线提示该模型具有较好的短期预测能力。通过GSVA分析发现,高危人群的致癌激活通路(ECM_RECEPTOR_INTERACTION、TGF_BETA_SIGNALING_PATHWAY、FOCAL_ADHESION等)显著富集,为高危人群预后不良提供了更全面、更有说服力的解释。最后,将风险和年龄确定为独立的OS相关变量,并将其纳入列线图,结果表明列线图可被视为GC患者临床诊断和治疗的有效工具。这些结果表明,基因标记具有较强预测GC患者预后的能力,对临床治疗决策具有一定的指导意义。此外通过CIBERSORT探讨了糖酵解相关基因与免疫浸润水平的相关性和免疫相关功能的相关性,揭示了15个糖酵解相关基因标记影响GC预后的潜在机制。我们发现高风险组中细胞浸润水平(巨噬细胞 M1)高于低风险组,而高风险组激活的CD4+T细胞、辅助性滤泡T细胞、单核细胞(均P<0.05)水平低于低风险组。显然,高风险组反映了一个免疫抑制的肿瘤微环境,相较于低风险组缺乏免疫细胞,与高风险组的预后不良相吻合。关于免疫相关功能的分析显示,高风险组CCR、Parainflammation、Type_I_IFN_Reponse、Type_II_ IFN_Reponse功能高于低风险组,而高风险组MHC_class_I低于低风险组。以上说明这15个糖酵解相关基因标记与GC患者的免疫细胞浸润与免疫相关功能密切相关。

综上,建立的15个糖酵解相关基因标记(PFKFB2、UHRF1、ACYP1、CLDN9、STC1、EFNA3、NUP50、ADH4、ANGPTL4、PKP2、VCAN、HIF1A、LHX9、、ANKZF1和ALDH3A2)可作为预测GC患者预后的可靠工具,可能为GC的潜在发病机制提供新的见解。