基于糖酵解风险基因老年结直肠癌患者预后风险模型构建*

何耀明， 李尚仁， 李晓平， 王成兴， 谭群英， 梁伟俊， 周朝荣， 赵景林△

(江门市中心医院 a.胃肠外科；b.乳腺外科，广东 江门 529000)

结肠癌是常见恶性肿瘤，根据2014年中国恶性肿瘤发病和死亡分析报告[1]，结肠癌发病率位于恶性肿瘤的第3位，每年新发病例数约37万，每年因结直肠癌死亡患者18万，位居第5位。随着我国肠镜检查的逐渐普及，结直肠癌的病死率不但无显著下降，而且高龄(≥65岁)患者病死率逐年上升,加上社会老龄化问题越来越严重，降低老年人结直肠癌病死率具有重要意义。肿瘤组织具有特殊的能量代谢特点，尽管在有氧环境下，肿瘤细胞始终以糖酵解这种快速产生ATP的糖代谢方式为细胞生长、增殖提供能量，即著名的Warburg效应[2]。目前多认为，糖酵解相关基因[如缺氧诱导因子1α(HIF-1α)[3]、葡萄糖转运蛋白-1(GLUT-1)[4]、己糖激酶2(HK-2)[5]等]在结直肠癌组织中表达增加，基因[如GLUT-1、乳酸脱氢酶(LDH-A)、丙酮酸激酶同工酶2(PKM2)等]在结肠癌肝转移瘤中比原发灶有更高表达[6]，均提示糖酵解相关基因与结肠癌发生、进展有重要作用。早发现、早治疗、早干预可降低结肠癌病死率，若每个结直肠癌术后患者能缩短检查周期，提早发现隐匿病灶并及时处理，对延长患者的生存率有重要意义，但过度检查会造成医疗资源浪费，所以一个预测模型的应用[7]，可将高危结肠癌患者区分开，从而有针对性地提早干预时间，是降低结直肠癌患者死亡的有效途径之一。然而，目前针对老年结直肠癌患者预后预测模型的报道较少，因此，本研究根据糖酵解相关基因结肠癌发生、进展的重要作用，基于TCGA结直肠癌数据库，确定与老年结直肠癌进展的相关风险基因及其相关系数，联合各项临床指标，建立老年结直肠癌的预后预测模型，从而为高危老年结直肠癌患者提早干预提供依据。

1 资料与方法

1.1 资料 从TCGA数据库获得结直肠癌转录组数据和相关临床资料(https://www.cancer.gov/) (Tomczak, Czerwinska & Wiznerowicz, 2015)。再进一步从中筛选出年龄≥65岁的结直肠癌患者相关数据用作分析，其中包含癌旁或正常组织32例，肿瘤组织340例，用于生存分析包含完整随访信息的患者321例，糖酵解相关基因集通过MsigDB数据库获得(http://www.hmdb.ca)。包括BIOCARTA_GLYCOLYSIS_PATHWAY，GOBP_FRUCTOSE_1_6_BISPHOSPHATE_METABOLIC_PROCESS，HALLMARK_GLYCOLYSIS，GOMF_LACTATE_TRANSMEMBRANE_TRANSPORTER_ACTIVITY，GOBP_LACTATE_TRANSMEMBRANE_TRANSPORT，KEGG_GLYCOLYSIS_GLUCONEOGENESIS，MODULE_306，REACTOME_GLYCOLYSIS，WP_GLYCOLYSIS_AND_GLUCONEOGENESIS 9个基因集。

1.2 方法

1.2.1 基因功能富集分析和差异分析获得糖酵解差异表达基因 上述糖酵解基因集通过GSEA 4.1.0分析软件进行结直肠癌肿瘤组织和正常组织的基因富集分析，分析结果以标准化富集分数(NES)>1.5和P值(P<0.05)设定临界值，获得相关糖酵解基因集。 从相关糖酵解基因集中提取基因，在TCGA数据集提取患者糖酵解相关基因的表达矩阵，表达矩阵使用“edgeR”包归一化处理，使用Log2FC>1和错误发现率(fales discovery rate,FDR)<0.05作为临界值，对肿瘤组织和正常组织相关糖酵解基因的表达量进行差异分析，获得糖酵解差异表达基因。

1.2.3 风险模型的效能验证 由TCGA数据获得的风险模型，通过单因素COX回归和多因素COX回归明确独立预后影响因子，包括：风险分数、性别、肿瘤分期、淋巴管浸润、脉管浸润。风险模型的验证分别在内部临床不同的亚组和外部GEO (https://www.ncbi.nlm.nih.gov/geo/)数据集GSE17536两方面进行，高低风险组结果对比使用Kaplan-Meier生存分析，ROC曲线计算AUCs值，评估风险模型的预测效能。

1.2.4 绘制列线图 根据多因素COX回归分析的结果，绘制分期，肿瘤位置(右侧：升结肠、横结肠；左侧：结肠脾曲、降结肠、乙状结肠、直肠)、分数的列线图，预测1、3、5 a的生存率，列线图通过C-index和校准图，临床决策曲线分析法，ROC曲线进行预测效能验证效能。

2 结果

2.1 老年结肠癌患者肿瘤组织和正常组织糖酵解相关基因集富集分析结果 通过GSEA富集分析9个糖酵解相关基因集发现，HALLMARK_GLYCOLYSIS和REACTOME_GLYCOLYSIS两个基因集在老年结肠癌患者的肿瘤组织对比正常组织，糖酵解相关基因富集差异显著。从这两个显著基因集中获得糖酵解相关基因254个，再与老年结直肠癌转录数据合并，得到老年结直肠癌糖酵解相关基因表达矩阵，进行下一步基因差异表达分析。见表1，图1。

表1 糖酵解基因集GSEA富集分析

注：A：HALLMARK_GLYCOLYSIS；B：REACTOME_GLYCOLYSIS图1 肿瘤组织与正常组织GSEA富集分析结果

2.2 糖酵解差异表达基因结果 从老年结直肠癌肿瘤基因表达数据和临床数据(340个肿瘤组织样品和32个癌旁正常组织样品)，结合获得的254个糖酵解相关基因，以Log2FC>1和FDR<0.05作为临界值进行差异分析，获得糖酵解差异基因46个，其中上调基因37个，下调基因9个。见图2。

注：A：热图展示254个糖酵解相关基因在肿瘤组织和正常组织中的表达情况；B：火山图展示老年结直肠癌肿瘤组织对比正常组织，37个基因显著上调，9个基因显著下调图2 老年结直肠癌糖酵解相关基因在340个肿瘤组织和32个癌旁正常组织的表达差异分析

2.3 糖酵解基因风险模型构建

2.3.1 筛选预后相关糖酵解基因 结合老年结直肠癌患者生存资料及差异表达糖酵解基因的表达矩阵，使用COX单因素回归方程，得到CAPN5、STC2、VEGFA、FKBP4、PKP2、NDC1、DEPDC1、COL5A1、TGFA、ANKZF1、GOT2 11个风险相关基因，再根据LASSO回归模型进一步优化，最终确认模型风险基因CAPN5、STC2、VEGFA、FKBP4、PKP2、NDC1、COL5A1、TGFA、ANKZF1、GOT2及其分子系数coefi。见图3。

注：A：单因素COX回归获得11个预后风险相关糖酵解基因；B-C：LASSO回归进一步优化确定10个风险相关基因及其分子系数图3 构建预后风险模型

2.3.2 构建风险模型 根据糖酵解相关基因的表达量和对应的分子系数，计算得出各个患者的风险分数，根据风险分数中位值将患者分为高风险组及低风险组，Kaplan-Meier生存分析提示高风险组患者生存期显著短于低风险组，P=0.000 1，ROC曲线分析，1、3、5 a生存率的AUC值分别为0.705、0.730、0.674。见图4。

注：A：根据患者风险分数由低到高排序分为高风险组和低风险组；B：高低风险组患者生存状态；C：高低风险组糖酵解预后基因的表达；D：高低风险组的Kaplan-Meier生存分析；E：ROC曲线图4 预后风险模型的构建

进行单因素和多因素COX回归，评价老年结直肠癌的影响预后的独立因素，结果发现，在单因素COX回归分析中，分期、脉管浸润、淋巴管浸润和风险分数是老年结直肠癌预后的独立因素，而在多因素COX回归分析中，只有分期和风险分数是影响老年结直肠癌预后的独立因素。见图5。

注：A：单因素COX回归；B：多因素COX回归图5 影响老年结直肠癌预后的独立因素

2.4 风险分数模型效能外验证 GEO数据GSE17536中年龄≥65岁患者进行风险模型外验证，外验证数据集中包括老年结直肠癌患者99例，根据TCGA老年结直肠癌风险分数中位数0.148 9作为临界值，将外验证数据集分为高风险组46例，低风险组53例。进行Kaplan-Meier生存分析提示高风险组患者生存期短于低风险组，差异有统计学意义(P=0.023)。见图6。

图6 外验证数据集GSE17536老年结直肠癌患者Kaplan-Meier生存分析

2.5 老年结直肠癌预测模型的构建 根据COX回归分析结果，基于患者风险得分，淋巴管浸润、脉管浸润、T分期、N分期、M分期构建列阵图，预测老年结直肠癌患者1、3、5 a生存率。使用C-index，校准曲线，ROC曲线，临床决策曲线(DCA)评估列阵图对老年结直肠癌患者的预测效能。预测模型C-index=0.730 48，校准曲线与理想曲线相近，ROC曲线评价提示模型预测1、3、5 a生存率的AUC值分别为0.766、0.742、0.783。均提示模型具有良好的预测效能。见图7。

注：A：构建列阵图预测老年结直肠癌患者生存率；B：预测模型的ROC曲线，预测1、3、5 a生存的AUC值分别是0.766、0.742、0.783；C：预测1、3、5 a生存率的校准曲线；D-F：1、3、5 a生存预测模型的临床决策曲线图7 老年结直肠癌预测模型的构建

3 讨论

结直肠癌是我国常见恶性肿瘤，虽然其患者的病死率基本与以往相仿，但是老年(≥65岁)结直肠癌患者病死率却逐年上升。肿瘤细胞具有特殊的代谢方式——Warburg效应，尽管在有氧条件下，肿瘤细胞的能量代谢仍以产生能量更快的糖酵解方式为主，糖酵解是肿瘤细胞侵袭转移的重要能量来源，而且目前多认为，糖酵解相关基因与结肠癌患者的预后有重要联系。因此，基于肿瘤的糖酵解代谢特征，本研究建立了老年结直肠癌糖酵解相关基因预后模型，有利于高危老年结直肠癌患者临床决策确立，及早干预，延长老年结直肠癌患者生存期。

糖酵解相关基因的风险分数与肿瘤分期一样，是老年结直肠癌预后不良的独立影响因素，高风险组患者，其生存率显著低于低风险组(P<0.001)，ROC曲线分析1、3、5 a生存率的AUC值分别是0.705、0.730、0.674，提示风险分数对老年1、3 a生存率有较好的预测效能，而对5 a生存率的预测效能一般，但当结合相关临床信息时，5 a生存率的AUC值提升至0.783，提示模型具有良好的预测效能。

目前对结直肠癌预后模型的使用报道不多，糖酵解相关基因预后模型更少见，曾有糖酵解相关基因结直肠癌预后模型的报道[8]，该模型主要建立在糖酵解预后不良上调基因的基础上，前期利用相同方法建立老年人结直肠癌预后模型发现，该方法无法在外数据集中得到验证，其中一个重要原因是TCGA数据库主要使用二代测序方法，与GEO探针检测方法存在本质区别，基因表达量差距较大，仅以TCGA数据库计算得出的由上调基因构成的高低风险临界值难以应用于外验证集中。本研究根据COX单因素回归方程所得风险基因再通过LASSO回归模型优化最终获得10个风险相关基因，包括提示不良预后的STC2、VEGFA、FKBP4、COL5A1、TGFA、ANKZF1和良好预后的CAPN5、PKP2、NDC1、GOT2，由预后良好基因和预后不良基因及其分子系数coefi共同计算得出风险分数可应用于GEO数据集，具有良好的预测效能。

VEGFA编码血管内皮生长因子α，调节肿瘤微环境中血管生成[9]，与结肠癌肿瘤进展、转移相关[10]，而STC2编码一种分泌型同型二聚体糖蛋白，参与调节多种蛋白活化相关，目前认为与结肠癌肿瘤大小、血管浸润、淋巴转移、肿瘤分化等多种病理临床指标具有相关性[11]，也有研究[12]认为，与转移结直肠癌VEGFR靶向治疗耐药相关。FKBP4是热休克蛋白90的共同伴侣分子， ANKZF1编码为锌指肽基TRNA水解酶，通常在细胞应激的状态下，FKBP4、ANKZF1均表达增加，维持线粒体的结构稳定性和完整性[13-14]。FKBP4在结肠癌肿瘤组织中表达增加[15]，有研究[16-17]表明，FKBP4表达的增强与乳腺癌、肺癌进展，预后不良相关，ANKZF1高表达提示结直肠癌患者预后不良[18]。TGFA，编码变形生长因子alpha，是EGFR的重要配体之一，能激活EGFR，调节上皮—间质迁移，促进结肠癌细胞转移和侵袭[19]。COL5A胶原蛋白成员，参与细胞外基质构成，在乳腺癌、卵巢癌、肺癌及肾癌及中高表达，COL5A表达增加与预后不良相关[20]。PKP2属于斑菲素蛋白Plakophilins (PKPs)家族，能与蛋白激酶C结合相关，PKP2是Wnt/β-catenin转录目标，可负反馈Wnt/β-catenin的转录活性，抑制肿瘤细胞中Wnt/β-catenin通路的过度激活[21]。而关于CAPN5、NDC1、GOT2在肿瘤中的作用报道很少，作用机制还有待进一步阐明。

本研究也存在一定的局限性，首先，文中的预测模型缺少患者的化疗、放疗、免疫治疗等治疗信息，可能会对预测效能造成一定影响；其次，本研究是基于生物信息学的分析，风险基因的相关功能与机制研究还有待于进一步明确。总之，本研究基于糖酵解相关风险基因建立的预测模型对老年结直肠癌患者预后有良好预测效能。