识别糖酵解相关LncRNA对头颈鳞状细胞癌的预后价值*

孙一驰， 曾宪琳， 张世超

(1.贵州医科大学 医学影像学院， 贵州 贵阳 550025； 2.贵州省细胞免疫治疗工程研究中心， 贵州 贵阳 550025； 3.贵州省普通高等学校感染免疫与抗体工程特色重点实验室， 贵州 贵阳 550025)

头颈鳞状细胞癌(squamous cell carcinoma of head and neck, SCCHN)是全球最常见的恶性肿瘤之一，每年夺去约350 000人的生命[1]，SCCHN局部复发、颈部淋巴结转移和耐药是晚期患者死亡的主要原因[2]。虽然SCCHN的治疗方法已有改进，但SCCHN患者预后较差，5年生存率仅为50%[3]。因此，如何对SCCHN患者进行早期诊断和有效治疗是亟待解决的问题。无论是否存在氧气，肿瘤细胞都主要通过糖酵解进行代谢。大量葡萄糖随着乳酸的产生而被消耗，这种现象称为有氧糖酵解或Warburg效应[4]。长链非编码 RNA (long noncoding RNA, lncRNA)被定义为一类不编码蛋白质的调节性RNA，其分子长度超过200 个碱基，在肿瘤发生发展中起着关键作用[5-6]。近年来的研究表明lncRNA在肿瘤代谢中起关键调控作用，并参与糖代谢通路[7]。例如，lncRNA中CDKN2B反义RNA 1 (CDKN2B antisense RNA 1，ANRIL)上调葡萄糖转运蛋白1 (glucose transporter 1, GLUT1)和乳酸脱氢酶A(Lactate dehydrogenase A, LDHA)的表达，从而增加葡萄糖摄取并促进鼻咽癌的发展[8]，LINC00092直接与6-磷酸果糖-2-激酶/果糖-2，6-二磷酸酶2(6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 2，PFKFB2)结合能增强糖酵解并促进肿瘤的发生和发展[9]。在膀胱癌中，lncRNA中过表达的尿路上皮癌胚抗原1 (urothelial cancer associated 1，UCA1)通过上调己糖激酶2 (Hexokinase 2, HK2) 促进糖酵解[10]。然而，参与SCCHN糖酵解重编程的lncRNA对SCCHN发生发展的作用仍不清楚。因此，本研究从TCGA数据集中获得了具有显著预后价值的lncRNA。并基于这些GRLs构建了预后风险模型，并确定SCCHN高/低风险组在信号通路、免疫微环境、免疫检查点、m6A相关基因方面的差异。

1 材料与方法

1.1 材料

从TCGA数据库(https://cancergenome.nih.gov/)下载SCCHN患者的转录组数据和临床数据[11]，转录组数据包含44例正常组织样本和502例头颈癌组织样，临床数据包含501例头颈癌患者的生存时间和生存状态。

1.2 方法

1.2.1数据预处理 利用Perl语言脚本整合SCCHN患者的转录组数据，获得Ensemble表达矩阵，同时，从Ensembl数据库(https://asia.ensembl.org/index.html)下载人类基因注释文件,将Ensemble表达矩阵转化为Gene Symbol表达矩阵，根据Gene Symbol表达矩阵与GENCODE数据库中lncRNA染色体的位置进行比较识别lncRNAs[12]，获得SCCHN的lncRNA表达矩阵。从MSigDB数据库中获得了总共274个糖酵解基因,利用Perl脚本从Gene Symbol中提取糖酵解基因表达矩阵。

1.2.2差异表达糖酵解基因和lncRNA筛选 使用R软件包limma筛选差异的糖酵解基因和lncRNA,筛选标准：|log2(fold-change)|>1以及矫正后的P(false discovery rate,fdr)<0.05，对于糖酵解基因表达矩阵和lncRNA表达矩阵筛选出的差异基因和差异lncRNA分别，利用R包ggplots绘制正常样本和肿瘤样本的火山图。

1.2.3差异糖酵解基因的富集分析 为了探索潜在糖酵解基因与SCCHN生物通路之间的关系和相关程度，使用 DAVID对差异的糖酵解基因进行GO-BP富集和KEGG通路分析，设定P<0.05为显著性基因富集[13]。

1.2.4糖酵解相关lncRNA 在R软件中，本研究使用糖酵解基因lncRNA共表达方法分析糖酵解基因和lncRNA的相关性，相关系数Cor>0.4且P<0.05的lncRNA被认为是糖酵解相关lncRNA。

1.2.5糖酵解相关lncRNA模型构建 通过单因素Cox回归(P<0.05)筛选预后糖酵解相关lncRNA后，通过多因素Cox回归确定关键的lncRNA[14]。基于以下公式构建SCCHN患者的风险评分模型，计算每个SCCHN患者的风险得分。风险评分=∑βlncRNA×ExplncRNA(βlncRNA表示多因素Cox回归分析中计算出的每个lncRNA的回归系数，Exp lncRNA代表样本中每个lncRNA的表达值。)

1.2.6生存分析 以SCCHN患者的中位风险得分作为临界值，进一步将SCCHN患者分为高危组和低危组(低危组得分<中位，高危组得分≥中位)。采用R包中的“Survival”绘制Kaplan-Meier生存曲线；使用R包“timeROC”包绘制ROC曲线。此外，通过校准曲线构建和列线图，以确定作为独立预后因素风险模型的准确性。

1.2.7模型的特点和应用 免疫微环境与SCCHN的发生发展密切相关。根据SCCHN患者的表达数据，通过 ESTIMATE 估计免疫和基质评分，以代表基质和免疫细胞的存在[21]；此外，使用两种算法 CIBERSORT、MCPcounter来比较不同风险组中各种免疫细胞比例的差异[15]；同时提取免疫检查点基因、m6A相关基因的表达水平，通过Wilcox检验分析它们在不同风险组中的表达差异。

1.3 统计学分析

R软件(v3.6.3)用于统计分析，Wilcox检验用于组间比较，Pearson用于分析糖酵解基因与lncRNA之间的相关性，单因素和多因素Cox回归分析影响SCCHN患者总生存期的相关因素；P<0.05为差异有统计意义。

2 结果

2.1 差异和富集分析

根据所设定的显著性阈值，从274个糖酵解基因中筛选得到76个糖酵解差异表达基因(图1A)，其中上调的有52个，下调的有24个。从lncRNA表达矩阵中筛选得到965个差异的lncRNA(图1B)，其中上调的有864个，下调的有101个。此外， 76 个糖酵解差异表达基因富集到GO和KEGG 通路上(图1C～D)。结果显示，大多数差异的糖酵解基因都富集于代谢路径。鉴定出的糖酵解基因与肿瘤发生发展中的几个重要生物学过程有关，如缺氧通路、代谢通路、AMPK信号通路、HIF-1信号通路，进一步证明糖酵解与肿瘤缺氧微环境密切相关。

注：A为差异表达的糖酵解相关mRNA的火山图，B为差异表达的糖酵解相关lncRNA的火山图；红色圆为上调，绿色圆为下调；C为差异的糖酵解基因GO分析结果、D为差异的糖酵解相关lncRNA的KEGG分析结果，色标表示P值，直方图大小表示计数。图1 糖酵解基因的差异和富集分析Fig.1 The differences and enrichment analyses of glycolysis-related genes

2.2 基于CRLs的风险模型构建与验证

对Pearson相关分析 (Cor>0.4,P<0.05)获得的298个GRLs进行单因素Cox回归分析，结果显示有60个GRLs与SCCHN患者的预后显著相关(图2A)。本研究进一步对这60个lncRNA进行多因素Cox回归分析，确定了25个关键的lncRNA用于构建SCCHN患者的预后风险模型。生存曲线展示了高风险的SCCHN患者预后较差(图2B)，高低风险组患者的存活状态、总生存期及lncRNA的表达均具有显著差异(图2C～E)。

注：A为预后相关lncRNA森林图，B为生存曲线，C为高/低风险SCCHN患者的RS分布图，D为SCCHN患者的生存时间散点图，E为SCCHN患者的生存热图。图2 lncRNA与SCCHN患者的生存和风险分析Fig.2 Survival and risk analysis of lncRNA in SCCHN patients

由临床特征和风险模型的单因素和多因素Cox 回归分析结果可知，风险评分P<0.001。因此，风险评分可作为SCCHN患者的独立预后因素(图3A～B)。同时本研究通过绘制ROC曲线来评估风险模型的可靠性，结果表明风险评分的曲线下面积(AUC)值高于年龄、性别、临床分期和组织学分级，表明风险评分对SCCHN患者预后的预测能力要优于传统的临床参数(图3C～D)。

注：A为风险模型的单因素分析，B为风险模型的多因素分析，C为 SCCHN患者1、2、3年生存期的ROC曲线，D为各临床指标与RS的ROC曲线比较。图3 临床特征和风险模型的Cox和ROC分析Fig.3 Cox analysis of clinical characteristics and ROC analysis of the risk model

基于临床特征和风险模型进一步构建决策曲线DCA(图4A)和列线图(图4B)。DCA曲线和列线图的结果都证明了构建模型的准确性。糖酵解相关lncRNAs表达与临床特征的热图(图4C)表明，糖酵解相关基因与预后具有密切关系，结果与构建模型预测的结果吻合。

注：A为不同临床因素绘制的决策曲线图，B为不同临床因素预后预测列线图，C为25种糖酵解表达水平与临床特征之间关联的热图。图4 风险模型的验证Fig.4 Validating the risk model

2.3 免疫细胞和免疫检查点分析

免疫检查点抑制剂可阻断表达免疫检查点肿瘤细胞对免疫细胞的抑制作用。鉴于此，为了进一步探索风险模型的临床应用，对相关的免疫细胞(图5A)和免疫功能(图5B)进行了分析，结果表明高/低风险组的免疫细胞浸润类型有明显差异，这为SCCHN的免疫治疗提供了理论依据。本研究比较了两个风险组之间免疫检查点基因的差异(图5C)可知，低风险组中免疫检查点基因表达量较高，说明RS较低的SCCHN患者对免疫检查点阻断(ICB)治疗更敏感。还比较了两个风险组之间m6A的表达水平，共匹配了12个m6A 调节因子，除YTHDF1(YTH N6-methyladenosine RNA binding protein 1)、HNRNPC(heterogeneous nuclear ribonucleoprotein C)、METTL3(methyltransferase 3, N6-adenosine-methyltransferase complex catalytic subunit)、YTHDF2(YTH N6-methyladenosine RNA binding protein 2) 和ZC3H13(zinc finger CCCH-type containing 13) 外，其余7个 m6A 调节因子的表达水平在低风险组中较高 (图5D)。

注：A为高低风险组免疫细胞热图，B为高低风险组免疫功能差异分析，C为高低风险组免疫检查点差异分析，D为高低风险组m6A差异分析。图5 高低风险组患者m6A免疫和基因分析Fig.5 Immunity and m6A analyses in high risk group

3 讨论

糖酵解是肿瘤细胞代谢的主要方式，可以为肿瘤细胞增殖提供所需的能量[16-17]，其最终产物乳酸被释放到细胞外，可维持酸性肿瘤微环境并促进毛细血管的形成，最终导致肿瘤细胞的增殖、耐药、侵袭和转移[18]。已有研究表明，糖酵解基因的上调与SCCHN的增殖、转移以及免疫逃逸有着密切关系[19-20]。然而，靶向SCCHN糖酵解途径的药物的作用机制尚不清楚，因此对其分子机制的深入研究仍具有重要意义。

已有研究表明，lncRNA在肿瘤的发生和发展中发挥着重要作用，且在肝细胞癌、结直肠癌和胃癌等多种癌症中，lncRNA可通过重编程葡萄糖代谢和糖酵解影响肿瘤的发生发展[21]。例如，lncRNA AGPG(actin gamma 1 pseudogene 25) 通过抑制 Lys302的泛素化和随后的蛋白酶体依赖性PFKFB3(6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3) 降解来增加 PFKFB3 的稳定性,并激活糖酵解途径，导致食管癌细胞的代谢重编程，促进食管癌细胞的迁移[22]。然而，在SCCHN中对糖酵解相关的lncRNAs的研究仍然很少。

因此，本研究通过Pearson分析、单因素Cox回归发现了GRLs与SCCHN患者预后显著相关，并建立了由25个关键的GRLs组成的风险模型。根据功能富集分析，本研究发现鉴定出的GRLs与HIF-1信号通路、氨基酸代谢和核苷酸代谢等有关。因此，推测上述lncRNA可能通过糖酵解诱导肿瘤微环境的缺氧、肿瘤细胞增殖和转移等。此外，本研究还通过ROC曲线对风险模型的可靠性进行了评估，发现AUC值分别为0.740、0.730及0.744，说明本研究的模型具有较好的预测性。之后本研究进一步探索了风险模型的临床意义，发现两个风险组表现出不同的免疫细胞浸润、免疫功能、免疫检查点、m6A调节因子水平，且高风险组与不良的临床表现和免疫功能低下显著相关。同时，Shen等[23]发现，肿瘤糖酵解诱导的缺氧和酸性微环境可导致代谢介导的T细胞功能障碍，Jung等[24]发现肿瘤的糖酵解可诱导肿瘤免疫抑制和免疫逃逸，这与本研究通过富集分析以及免疫功能差异分析得到的结果一致。因此，靶向代谢的肿瘤免疫治疗策略有望提高免疫治疗的疗效，本研究也希望通过进一步的研究探索糖酵解(通过m6A修饰或免疫检查点)的分子机制，以提高免疫治疗的疗效。

综上所述，本研究中筛选出的糖酵解相关的lncRNAs在 SCCHN的发生和进展中起着重要作用，基于糖酵解相关lncRNAs的风险模型是SCCHN患者的独立预后因素。所构建的基于糖酵解相关lncRNAs的风险模型是SCCHN患者的独立预后因素。通过综合分析，糖酵解相关的lncRNAs的模型为SCCHN的临床治疗提供了一定的指导建议。但本研究仍存在一些局限性。首先，由于可以注释lncRNA表达的SCCHN样本数量有限，需要更多具有同源信息的患者纳入研究并证明本研究的可信度。其次，本研究仅通过生物信息学分析探讨了lncRNA的功能，需要实验数据来支持这些结论。尽管有这些限制，本研究还是使用了ROC和列线图来证明风险模型预测的有效性。