識別糖酵解相關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和列線圖來證明風險模型預測的有效性。