基于乳酸代谢相关基因的头颈部鳞状细胞癌分子亚型和临床特征的生物信息学分析

杨紫煦, 苏 畅, 王波元, 刘 冲, 李明贺,2

（1. 吉林大学口腔医院口腔颌面外科，吉林 长春 130021；2. 吉林大学第二医院口腔科，吉林 长春 130022）

头颈部鳞状细胞癌（head and neck squamous cell carcinoma，HNSCC） 是全球第六大常见癌症［1］，且发病率逐年上升，预计至2030 年将增加30%，即每年增加108 万新发病例［2-3］，患者5 年生存率仅为40%～50%［4］，晚期HNSCC 患者生存率为34.9%［5］。目前，HNSCC 患者的主要治疗方法为手术治疗，放疗和化疗可作为晚期和肿瘤转移患者的辅助疗法。HNSCC 患者的预后并不理想，因此研究HNSCC 的分子机制，确定新的治疗靶点，构建有价值的预后模型，可为预测HNSCC 患者预后和个体化用药提供参考，以改善其治疗效果。

肿 瘤 微 环 境 （tumor microenvironment，TME）中乳酸的增加有利于癌细胞增殖和免疫逃逸。乳酸不仅是糖酵解的产物，也是正常组织和恶性肿瘤组织之间的关键调节剂。肿瘤转移、复发和预后较差均与较高的乳酸水平有关联［6-7］。近年来，抑制乳酸代谢已被证实为针对恶性肿瘤的潜在治疗方法［8］。研究［9-12］表明：乳酸代谢相关基因（lactic acid metabolism-related genes，LRGs）在食管鳞状细胞癌、肺腺癌、乳腺癌和肾透明细胞癌等多种癌症的发生发展过程中发挥重要作用。但LRGs 在HNSCC 中的作用尚不明确。本研究对LRGs 表达谱进行综合评估，基于2 种LRGs 亚型鉴定差异表达基因（differentially expressed genes，DEGs），构建HNSCC 患者预后模型，以独立预测HNSCC 患者的生存情况，同时探讨 LRGs 与HNSCC 患者免疫细胞浸润和肿瘤相关评分的相关性，为HNSCC的临床诊断和治疗提供理论依据和思路。

1 资料与方法

1.1 数据收集和预处理由癌症基因图谱（The Cancer Genome Atlas， TCGA） 数据库下载HNSCC 患者基因表达数据和临床数据作为训练队列，其中HNSCC 组织样本522 份，正常组织样本44 份，并下载肿瘤突变谱。由基因表达综合（Gene Expression Omnibus，GEO） 数据库下载GSE27020、GSE41613 和GSE65858 数据集作为验证队列。由GeneCards［12］数据库提取5 034 个LRGs。

1.2 LRGs 差异表达分析采用limma 数据包对TCGA 数据库中HNSCC 组织和癌正常组织中LRGs mRNA 表达水平进行分析（|log2FC|＞1，P＜0.05），通过火山图和热图显示差异表达的LRGs。

1.3 聚类分析和TME 分析采用单因素Cox 回归分析筛选出与预后相关的差异表达LRGs（differentially expressed-LRGs，DE-LRGs）。对筛选出的DE-LRGs 进行Pearson 相关分析。基于DE-LRGs mRNA 表 达 水 平， 采 用ConsensusClusterPlus 数据包对TCGA 数据库中HNSCC 患者采用K-means 算法进行聚类分析。聚类分析是通过迭代寻找K 个分组的划分方案，在不同分组中找到数据均值，通过最小化损失函数确定分组的个数及对应的成员。损失函数：‖xi-μci‖2，其中xi为第i个患者对应的基因表达水平，ci为xi所属的组，μci为分组对应的均值，n为患者总数。采用Kaplan-Meier （K-M）生存曲线分析比较不同分组患者的预后。 采用CIBERSORT 算法进行TME 分析，预测K 个分组患者肿瘤样本中浸润免疫细胞的组成。

1.4 预后模型的构建和验证采用多因素Cox 回归分析和LASSO 回归分析对HNSCC 的独立预后基因进行分析，构建预后模型，根据风险评分的中位数，将HNSCC 患者分为高风险组和低风险组。采用K-M 生存曲线和受试者工作特征曲线（receiver operating characteristic curve，ROC）判断预后模型的预测价值，采用GSE27020、GSE41613和GSE65858 数据集验证预后模型。

1.5 肿瘤免疫细胞浸润分析采用单样本基因集富集分析（single-sample gene set enrichment analysis，ssGSEA）算法分析各样本的免疫细胞活性和免疫功能。采用CIBERSORT 算法预测低风险组和高风险组患者肿瘤样本中浸润免疫细胞的组成，比较不同分组患者主要组织相容性复合体（major histocompatibility complex，MHC）的表达情况。

1.6 肿瘤相关评分与风险评分相关性分析计算低风险和高风险组HNSCC 患者血管生成活性评分、EMT 评分、致瘤细胞因子评分和干性评分，采用Pearson 相关分析方法分析上述4 种肿瘤相关评分与风险评分的相关性。

1.7 统计学分析采用R统计软件（4.2.0版本）进行统计学分析。采用“Survival”软件包进行K-M生存分析、单因素Cox 分析和多因素 Cox 分析，采用“glmnet”软件包进行LASSO 回归分析，采用timeROC 软件包绘制ROC 曲线，根据ROC 曲线评估患者预后模型的准确性，采用GSVA 软件包进行ssGSEA分析。以P＜0.05为差异有统计学意义。

2 结 果

2.1 LRGs 差异表达分析结果基于GeneCards 数据库获得5 034 个LRGs。分析TCGA 数据库中HNSCC 组织和癌旁正常组织样本的表达数据，获得1 196 个DE-LRGs。见图 1。

图1 HNSCC 组织和癌旁正常组织中LRGs mRNA 表达水平Fig. 1 Expression levels of LRGs mRNA in HNSCC tissue and paracancerous normal tissues

2.2 聚类分型的鉴定及其与TME 的相关性分析单因素Cox 回归分析得到27 个HNSCC 预后相关基因，其中16 个为保护因素， 11 个为危险因素（图 2A）。相关分析结果显示：大多数预后相关基因之间存在明显相关性（图 2B）。进一步利用27 个预后相关基因进行聚类分型，当HNSCC 患者分为2 个亚组时，聚类效果最好，亚组内部一致性和稳定性较好（图3A～C）。K-M 生存曲线分析结果显示：分组2 患者的预后优于分组1（图3D）。CIBERSORT 算法结果显示：与分组1比较，分组2患者免疫细胞浸润水平更高，CD4+T 淋巴细胞、CD8+T 淋巴细胞、中性粒细胞、巨噬细胞和树突状细胞浸润水平差异有统计学意义（P＜0.05）。见图4。

图2 HNSCC 预后相关基因单因素Cox 回归分析的森林图（A）和相关性分析（B）Fig. 2 Univariate Cox regression analysis forest plot (A) and correlation analysis (B) on HNSCC prognosis-related genes

图3 聚类分析图（A～C）和聚类后HNSCC 患者K-M生存曲线（D）Fig. 3 Cluster analysis diagrams (A—C) and K-M survival curve (D) of HNSCC patients after cluster

图4 聚类分型后免疫相关分析Fig. 4 Immuno-correlation analysis after cluster typing

2.3 HNSCC 预后模型的构建采用多因素Cox回归分析和LASSO 回归分析对HNSCC 患者的独立预后基因进行分析，构建预后模型，选择9 个基因加入模型，包括次黄嘌呤磷酸核糖基转移酶 1（hypoxanthine phosphoribosyltransferase 1，HPRT1）、淀粉样蛋白前体蛋白 （amyloid precursor protein， APP）、 糖原磷酸化酶 L （glycogen phosphorylase L，PYGL）、尿激酶型纤溶酶原激活物（urokinase-type plasminogen activator，PLAU）、大麻素受体2（cannabinoid receptor 2，CNR2）、斯钙素2 （stanniocalcin 2，STC2）、核苷酸结合寡聚化结构域样受体蛋白 1 （nucleotide-binding oligomerization domain-like receptor protein 1，NLRP1）、整合素连接激酶（integrin-linked kinase，ILK）和叉头框蛋白B1（forkhead-box B1，FOXB1）。风险评分公式如下：风险评分=0.015 195×HPRT1 表达水平+0.002 176×APP 表达水平+0.005 053 ×PYGL 表达水平+0.002 458×PLAU 表达水平-0.983 8×CNR2 表达水平+0.023 846×STC2 表达水平-0.090 1×NLRP1表达水平-3.414 72×ILK 表达水平+0.641 958×FOXB1 表达水平。 根据风险评分中位数将HNSCC 患者分为低风险组和高风险组。

TCGA 数据集中，与低风险组比较，高风险组HNSCC 患者预后较差，1、3 和5 年的曲线下面积（area under curve，AUC）分别为0.699、0.724和0.660 （图5）。采用GSE41613、GSE65858 和GSE27020 数据集验证，预后模型具有良好的效率（图6），在2 个临床亚组中，与低风险组比较，高风险组HNSCC 患者预后较差。

图5 TCGA 队列的生存分析、风险评分、生存状态和ROC 分析Fig. 5 Survival analysis, risk scores, survival status， and ROC analysis in TCGA cohort

图6 GSE41613、GSE65858 和GSE27020 数据集中HNSCC 患者的生存分析、风险评分、生存状态和ROC 分析Fig. 6 Survival analysis, risk scores, survival status, and ROC analysis of HNSCC patients in GSE41613 ，GSE65858,and GSE27020 datasets

2.4 HNSCC 患者肿瘤免疫细胞浸润分析ssGSEA 结果显示：低风险组HNSCC 患者较高风险组患者有更高的免疫细胞浸润和更多的免疫相关通路（图7）。低风险组HNSCC 患者浸润免疫细胞主要包括 B 淋巴细胞、浆细胞、CD8+T 淋巴细胞、肥大细胞和树突状细胞（图8A），低风险组HNSCC 患者MHC 表达水平较高（图8B）。

图7 ssGSEA 算法分析高风险组和低风险组HNSCC 患者肿瘤免疫细胞浸润情况Fig. 7 Infiltration of tumor immune cells of HNSCC patients in high and low risk groups analyzed by ssGSEA algorithm

图8 CIBERSORT 算法分析高风险组和低风险组HNSCC 患者肿瘤免疫细胞浸润情况（A)和MHC 表达情况（B)Fig. 8 Infiltration of tumor immune cells (A) and expressions of MHC(B) of HNSCC patients in high and low risk groups analyzed by CIBERSORT algorithm

2.5 HNSCC 患者肿瘤相关评分与风险评分的相关性高风险组HNSCC 患者血管生成活性和致瘤细胞因子评分较高，低风险组HNSCC 患者干性评分较高（图9），血管生成活性和致瘤细胞因子评分与风险评分呈正相关关系（r=0.85，P=8.8e-16；r=0.24，P=5.3e-08），干性评分与风险评分呈负相关关系（r=-0.12，P=0.006 1）。

图9 高风险组和低风险组HNSCC 患者肿瘤相关评分Fig. 9 Tumor-related scores of HNSCC patients in high and low risk groups

3 讨 论

肿瘤组织中乳酸的积累是恶性肿瘤发展的关键和早期事件［13］。近年来，乳酸代谢受到越来越多的关注，但HNSCC 与乳酸代谢相关性分析的研究较少。本文作者分析了乳酸代谢与HNSCC 发展的相关性，以确认其在癌症生物学中的作用。

本研究首先通过TCGA 和GeneCards 数据库获取1 196 个在HNSCC 中差异表达的LRGs。 进一步对HNSCC 患者的总生存时间（overall survival，OS）进行单因素Cox 回归分析，得到了27 个预后相关基因，对预后基因进行进一步优化，纳入9 个关键基因建立HNSCC 预后模型，包括HPRT1、APP、PYGL、PLAU、CNR2、STC2、NLRP1、ILK 和FOXB1，并在3 个外部队列中进行验证。此外，通过聚类分析将HNSCC 患者分为2 组，与分组1 比较，分组2 患者的预后更好，其免疫细胞浸润水平更高，可能与分组2 患者OS 更长有关。

本研究中，根据风险评分中位数将患者分为高风险组和低风险组，2 组患者在预后和TME 等方面均存在明显差异，表明该预后模型可用于指导HNSCC 患者的治疗，具有预测HNSCC 患者预后的价值。

近年来，免疫疗法在治疗皮肤癌、膀胱癌、肺癌、肾癌和错配修复缺陷肿瘤方面显示出明显效果［14］。因此，寻找可靠的生物标记物对HNSCC 患者进行免疫治疗非常重要［15］。

HPRT1 在多种组织中高表达，与HNSCC 患者预后不良有明显关联［16］，可能是HNSCC 潜在的治疗靶点。STC2 是一种分泌糖蛋白，在调节钙稳态、葡萄糖稳态和磷转移中起重要作用，是一种潜在的泛癌预后标志物和新的免疫治疗靶点［17］。NLRP1 是一种先天性免疫传感器，与肺腺癌和胰腺癌［18］等疾病发生有关。

APP 对癌症免疫疗法具有耐药性，可能与使MHC-I APP 失活的基因发生突变有关［19］。PYGL是结肠癌发生的风险基因［20］，还与食管鳞癌［21］、胰腺癌［22］和胶质瘤［23］的发生有关。PLAU 可能是葛根芩连汤治疗溃疡性结肠炎的靶基因［24］，还与胰腺癌的预后［25］及甲状腺癌［26］、胃癌［27］和乳腺癌［28］的发生有关联。CNR2 可通过内源性大麻素损害特异性T 淋巴细胞的功能，导致肿瘤患者的OS 缩短［29］。ILK 与整合素和生长因子受体信号通路有关，由于2 种信号通路均有助于癌细胞的抵抗，ILK 可能是潜在的新的肿瘤治疗靶点［30］。FOXB1 是叉头框（forkhead-box ，FOX） 基因家族的成员， FOX 家族基因的失序会导致癌症的发生，全面研究FOX 家族基因的表达谱、基因改变和表观遗传变化，有助于开发新的人类疾病治疗方法和预防手段［31］。基于上述9 个基因建立HNSCC预后模型， 通过TCGA 数据库和GSE65858、GSE27020 及GSE41613 数据集验证，高风险组HNSCC 患者OS 较短，从预后模型相关基因的回归分析来看，9 个基因均具有较好的诊断性能。因此该预后模型可以独立预测HNSCC 患者的总体生存率。

在本研究中，低风险组HNSCC 患者免疫细胞表达更丰富，提示其对免疫治疗更敏感。本研究结果表明：乳酸代谢对TME 有一定的调节作用，可能有助于发现肿瘤免疫的调控机制，为TME 的研究提供新的思路。

综上所述，本研究中由9 个LRGs（HPRT1、APP、PYGL、PLAU、CNR2、STC2、NLRP1、ILK 和FOXB1） 构建的HNSCC 的LRGs 预后模型，可以预测HNSCC 患者的生存情况和治疗反应。

利益冲突声明：

所有作者声明不存在利益冲突。

作者贡献声明：

杨紫煦参与研究设计、数据收集整理和论文撰写，苏畅参与研究设计，王波元参与数据收集和分析，刘冲参与论文校对和图片整合，李明贺参与论文审校。