1. 前言
头颈部鳞状细胞癌是世界上最常见的恶性肿瘤之一，是一组具有异质性的上皮源性肿瘤，主要发生在口腔、唇部、咽喉部、鼻窦和其他头颈部位 [1]。有发病率高、死亡率高、复发率高等特点。HNSCC最常见的危险因素是饮酒、吸烟和HPV病毒感染。目前的治疗方式是放疗、化疗、手术联合的综合治疗方案。远处转移率高(20%~30%)、局部复发率高(50%~60%)和放化疗抵抗等特点增加了HNSCC患者的治疗难度 [2] [3]。据统计，HNSCC的发病率在逐年上升，患者5年生存率 < 50%，局部进展期的复发或转移的患者生存期 < 1年。由于HNSCC的病因复杂性和异质性，并且缺乏早期诊断和预测的生物标志物导致HNSCC患者预后不良，因此寻找有助于预测HNSCC的预后和指导治疗的新型分子标志物具有重要的意义。
细胞代谢重编程是细胞为满足能量的需求，通过改变代谢模式促进肿瘤细胞的生长和增殖，被认为是肿瘤细胞独特的标志之一。Warburg Effect是指肿瘤细胞改变能量代谢方式使自身迅速适应低氧、酸性等不利于其生长的肿瘤微环境而快速增殖的现象，代谢重编程主要包括糖酵解、脂质代谢、氨基酸代谢等。本研究是探讨糖酵解基因与头颈部鳞状细胞癌患者预后之间的关系 [4]。研究发现，随着糖酵解程度的增加，不同的信号分子通路被激活，促进了肿瘤细胞的增殖、侵袭和转移，也增强了肿瘤细胞的耐药性 [5]，提示我们糖酵解对肿瘤的预后、诊断和治疗有重要的意义。
本研究通过下载TCGA数据库中头颈部鳞状细胞癌的转录组数据和临床信息，在GSEA数据库中下载了5个糖酵解相关基因集，通过基因集富集分析筛选出254个表达量有显著差异的糖酵解相关基因，进行后续分析。我们使用单因素、多因素COX回归分析、LASSO回归分析等方法构建了由16个糖酵解基因(STC2、STC1、GNPDA1、EXT2、PRPS1、HK1、GMPPB、TXN、HMMR、CHPF、PGK1、PYGL、SOD1、CTH、GALE、PLOD2)组成的HNSCC的预后风险模型。此外，我们证实该预后模型可以作为HNSCC的独立预测因子，高风险组患者的预后明显要差于低风险组患者的预后，并且预测性能要高于其他临床特征。
2. 方法
2.1. 数据下载
从TCGA数据库(https://portal.gdc.cancer.gov/)中下载头颈部鳞状细胞癌的转录组数据(n = 546)和临床信息数据(n = 528)。临床信息包括患者的年龄、性别、肿瘤的TNM分期、分级。
2.2. GSEA富集分析
GSEA官网(http://www.gsea-msigdb.org/gsea/index.jsp)中下载糖酵解基因集，通过生物信息学的方法研究HNSCC肿瘤组织和正常组织中糖酵解基因集是否有显著差异，筛选标准为P < 0.05。将基因集中有显著差异的基因用于后续的研究。
2.3. 糖酵解基因预后风险模型的构建
使用R软件进行单因素、多因素COX回归、LASSO回归分析识别与HNSCC预后相关的糖酵解基
因，并构建预后风险模型。风险评分的计算公式是风险得分=
，使用survival包进行Kaplan-Meier生存分析，比较高低风险组的总生存期的差异。采用受试者工作特征曲线评估该模型的区分能力。我们使用多因素COX回归分析评估该模型的独立预后价值，将风险评分与临床特征相结合构建动态列线图，通过校准曲线评估该列线图的预测性能。
2.4. 统计学分析
WilCoxon检验分析肿瘤组织和正常组织差异表达基因，过滤标准设定为|log2Fc| ≥ 1和FDR < 0.05。采用Mann-Whitney对高低风险组的免疫细胞、免疫成分和免疫功能的ssGSEA评分进行比较，采用Benjamini-Hochberg法调整P值。使用Kaplan-Meier方法比较不同组间的总生存期(OS)。使用Pearson相关分析评估糖酵解预后基因表达水平或预后模型风险评分和肿瘤干性评分、药物敏感性、免疫和基质评分之间的相关性。使用R软件中的ggplot2、pheatmap、venn、survminer、igraph、survival包进行绘图。
3. 结果
3.1. GSEA基因集筛选
使用GSEA4.2.1软件进行糖酵解基因集富集分析，明确糖酵解基因在HNSCC肿瘤组织和正常组织中是否有显著差异。结果显示在HALLMARK_GLYCOLYSIS和REACTOME_GLYCOLYSIS两个基因集中存在显著富集，P < 0.05 (图1(a)、图1(b))。从以上基因集中选择254个基因进行后续分析。
3.2. 糖酵解预后风险模型的构建与评估
我们使用单因素COX回归分析、多因素COX回归分析和LASSO回归分析构建糖酵解预后风险模型，结果显示16个糖酵解基因(STC2、STC1、GNPDA1、EXT2、PRPS1、HK1、GMPPB、TXN、HMMR、
(a) 
(b)
Figure 1. GSEA enrichment pathway analysis
图1. GSEA富集通路分析
CHPF、PGK1、PYGL、SOD1、CTH、GALE、PLOD2)与HNSCC患者的总生存期相关(P < 0.05)。风险评分的计算公式如下：风险评分 = (STC2的表达量*0.048) + (STC1的表达量*0.011) + (GNPDA1的表达量*0.175) + (EXT2的表达量*0.146) + (PRPS1的表达量*0.052) + (HK1的表达量*0.277) + (GMPPB的表达量*−0.179) + (TXN的表达量*0.061) + (HMMR的表达量*0.067) + (CHPF的表达量*0.100) + (PGK1的表达量*0.206) + (PYGL的表达量*0.134) + (SOD1的表达量*0.221) + (CTH的表达量*0.014) + (GALE的表达量*0.129) + (PLOD2的表达量*0.022)。根据中位截断值，将HNSCC患者分为高、低风险两组(图2(a))，每个患者的生存时间如图2(b)所示。我们使用PCA和t-SNE的方法进行降维，区分方式主要是两种(图2(c)、图2(d))。K-M生存分析结果显示低风险组的总生存期要明显高于高风险组(图2(e))，ROC曲线结果显示该模型1、3、5年的区分度分别是0.666、0.709和0.742 (图2(f))。为避免过度拟合，我们使用LASSO回归分析，最终得到了16个预后风险基因用于构建预后模型(图3(a)、图3(b))。
3.3. 糖酵解相关基因预后风险模型的独立预后价值
我们为了确定风险评分是否可以作为HNSCC的独立预测因子，我们采用单因素和多因素COX回归分析的方法进行验证。单因素COX回归分析结果显示，风险评分与总生存期呈明显的相关(图4(a))，HR = 1.935，95% CI = 1.513~2.476，P < 0.001；多因素COX回归分析的结果显示，风险评分依旧可以作为HNSCC的独立预测因子HR = 1.926，95% CI = 1.490~2.490，P < 0.001。
3.4. 列线图的构建与评估
为了增强该风险模型的预测性能，我们将风险评分与临床特征(患者性别、年龄、肿瘤分级、肿瘤分期、T分期、N分期)相结合，构建出一个动态列线图预测模型，并绘制校准曲线评估该模型的预测值和实际值的差异程度，结果如图5(a)所示。校准曲线结果显示该列线图在第1、3、5年的预测准确性较好(图5(b))。
(a) 
(b)
Figure 3. LASSO regression analysis. (a) Partial likelihood deviation plots were drawn by LASSO regression analysis; (b) LASSO regression coefficient analysis
图3. LASSO回归分析。(a) LASSO回归分析法绘制部分似然偏差图；(b) LASSO回归系数分析
(a) 
(b)
Figure 4. Construction of a glycolytic prognostic risk model for head and neck squamous cell carcinoma. (a) Univariate independent prognostic analysis; (b) Multivariate independent prognostic analysis
图4. 头颈部鳞状细胞癌糖酵解预后风险模型的构建。(a) 单因素独立预后分析；(b) 多因素独立预后分析
3.5. 预后模型基因突变和差异分析
使用cBioPortal数据库对参与构建模型的糖酵解基因进行突变情况分析，结果显示，GNPDA1、HMMR、PLOD2、PGK1、PRPS1、PYGL、SOD1、STC2、CTH、EXT2、CHPF、STC1在肿瘤组织中的表达量明显高于正常组织(图6(b)~(d)、图6(f)~(i)、图6(k)~(m)、图6(p)、图6(q))；而HK1、TXN、GALE、GMPPB在肿瘤组织中的表达量明显低于正常组织(图6(e)、图6(j)、图6(n)、图6(o))。突变情况如图6(a)所示。
(a)
(b)
Figure 5. Construction and validation of nomogram. (a) Line segment dynamic nomogram; (b) Calibration curve
图5. 诺莫图的构建和验证。(a) 线段式动态诺莫图；(b) 校准曲线
3.6. 数据分层分析
在TCGA数据库中，根据临床病理特征(年龄、性别、肿瘤分级、肿瘤分期、T分期、N分期、M分期)进行了亚组分析，根据K-M生存曲线，在按年龄、性别、肿瘤分级、肿瘤分期、T/N/M分期进行分层的HNSCC患者中，风险评分仍然可以作为一个预测预后的稳定因素(图7)。
4. 讨论
头颈部鳞状细胞癌是一种比较常见的恶性肿瘤，发病机制复杂、复发率高等特点使得患者的发病率逐年上升，5年生存率 < 50%。由于表型的复杂性，预测HNSCC患者的预后难度增加。预后模型的构建对临床决策很有帮助，更有利于制定个体化的治疗方案，进行精准治疗，以延长患者的生存时间。Warburg反应是指肿瘤细胞萄糖的利用方式从氧化磷酸化转变为糖酵解的转变，被认为是癌症的一大关键特征，肿瘤细胞糖酵解的增强是由于此反应中关键酶的表达或活性增强。所以基于糖酵解基因构建HNSCC的预后模型具有重要意义。
本研究首先通过GSEA富集分析筛选了在头颈部鳞状细胞癌患者中表达量具有显著差异的糖酵解基因集，最终确定了两个基因集，并纳入后续分析。通过LASSO回归分析、单因素和多因素COX回归分析，最终确定16个有价值的糖酵解基因(STC2、STC1、GNPDA1、EXT2、PRPS1、HK1、GMPPB、TXN、HMMR、CHPF、PGK1、PYGL、SOD1、CTH、GALE、PLOD2)用于构建预后风险模型。STC家族是一种糖蛋白激素，参与各种生物过程，包括钙磷平衡、细胞增殖、细胞凋亡、氧化应激和肿瘤发生发展等 [6]。STC2 (斯钙素2)经研究发现在结直肠癌中可以作为有意义的血清分子标记物用于疾病的诊断，STC1 (斯钙素1)在诊断和预测乳头状甲状腺癌的预后中的作用被证实 [7]。GNPDA1 (葡萄糖6-磷酸脱氨酶1)催化葡萄糖胺-6-磷酸转化为6-磷酸–果糖，将己胺系统与糖酵解途径联系起来，提供能量。研究证实，在肝细胞癌中，GNPDA1大量表达，可以作为肝细胞癌的诊断和治疗新的靶点 [8]。EXT2被发现在多发性外生性骨瘤中存在突变，可以作为遗传咨询和产前诊断重要的生物分子标志物 [9]。PRPS1可以控制嘌呤的生物合成、为DNA和RNA的生物合成提供前体物质，而且可以通过底物R-5P的反馈通路调节磷酸戊糖途径。研究发现PRPS1的磷酸化促进核苷酸的合成和肿瘤的发生 [10]。HK1作为糖酵解途径的第一个限速酶，在癌症中具有重要作用，在卵巢癌中，HK1的高表达与患者的预后不良相关 [11]。GMPPB经研究发现与先天性糖基化紊乱疾病的发病密切相关 [12]，TXN是一种典型的氧化还原蛋白，可以通过多种方式调节细胞的氧化还原状态，在结直肠癌中，对于放疗抵抗的患者，使用TXN抑制剂联合放疗可以作为一种有效的治疗措施，有助于提高放疗敏感性 [13]。HMMR在结直肠癌中表达明显上调，可以促进癌症的生长、侵袭和远处转移 [14] [15]。CHPF的转录表达和DNA甲基化与免疫浸润和免疫标记物相关，CHPF的上调促进了肿瘤细胞的恶性行为，并与乳腺癌患者预后不良相关 [16]。PGK1催化1,3-BPG和ADP转化为3-PG和ATP，是糖酵解的关键过程，并且在乳腺癌 [17]、胰腺导管腺癌 [18]、子宫内膜癌 [19]、结肠癌中明显上调 [20]。PYGL被发现在小鼠中PYGL的缺乏会导致糖原累积病以及Yap依赖性的肝脏肿大和肿瘤发生 [21]。SOD1可以调节KRAS突变的非小细胞肺癌的核糖体的生物合成 [22]。CTH可以保护乳腺癌细胞不受活性氧和化疗药物诱导的细胞凋亡的影响，被认为是乳腺癌诊断和治疗的新型生物分子标记物 [23]。PLOD2在食管鳞状细胞癌中高表达，与肿瘤的侵袭和转移密切相关 [24]。GALE是半乳糖代谢的关键酶，在胶质母细胞瘤的发病中发挥重要作用 [25]。
我们通过K-M生存分析发现高风险组患者的总生存期要明显低于低风险组。多因素COX回归分析结果显示，风险评分可以作为独立的预测因子预测HNSCC患者的预后。我们对风险评分和临床特征进行整合构建出动态列线图，以提高该预后模型的预测性能，并通过校准曲线进行验证，结果显示列线图在预测患者预后方面具有良好的准确性。通过数据分层分析发现与性别、年龄、肿瘤分级、分期等临床特征相比，风险评分依旧是很稳定的预测因子。我们还通过cBioPortal数据库分析16个基因的突变状态，不仅帮助我们了解HNSCC的遗传学特点，还为HNSCC精准治疗提供方向。
我们的研究存在一些局限性。首先，所有的数据都是在TCGA公共数据平台中获得的，属于回顾性的研究；其次，我们需要进行大规模的队列研究，用于验证此模型的预测性能以及临床适用性。我们还可以通过基础实验探究这16个基因在头颈部鳞状细胞癌预后中的作用。
综上所述，基于16个糖酵解基因构建的HNSCC的预后风险模型，可以为HNSCC患者的诊断和治疗提供个体化方案，为临床决策提供指导，为未来HNSCC发病机制的研究提供新的靶点和方向。