孫一馳， 曾憲琳， 張世超

(1.貴州醫(yī)科大學(xué) 醫(yī)學(xué)影像學(xué)院， 貴州 貴陽 550025； 2.貴州省細(xì)胞免疫治療工程研究中心， 貴州 貴陽 550025； 3.貴州省普通高等學(xué)校感染免疫與抗體工程特色重點(diǎn)實(shí)驗(yàn)室， 貴州 貴陽 550025)

頭頸鱗狀細(xì)胞癌(squamous cell carcinoma of head and neck, SCCHN)是全球最常見的惡性腫瘤之一，每年奪去約350 000人的生命[1]，SCCHN局部復(fù)發(fā)、頸部淋巴結(jié)轉(zhuǎn)移和耐藥是晚期患者死亡的主要原因[2]。雖然SCCHN的治療方法已有改進(jìn)，但SCCHN患者預(yù)后較差，5年生存率僅為50%[3]。因此，如何對SCCHN患者進(jìn)行早期診斷和有效治療是亟待解決的問題。無論是否存在氧氣，腫瘤細(xì)胞都主要通過糖酵解進(jìn)行代謝。大量葡萄糖隨著乳酸的產(chǎn)生而被消耗，這種現(xiàn)象稱為有氧糖酵解或Warburg效應(yīng)[4]。長鏈非編碼 RNA (long noncoding RNA, lncRNA)被定義為一類不編碼蛋白質(zhì)的調(diào)節(jié)性RNA，其分子長度超過200 個(gè)堿基，在腫瘤發(fā)生發(fā)展中起著關(guān)鍵作用[5-6]。近年來的研究表明lncRNA在腫瘤代謝中起關(guān)鍵調(diào)控作用，并參與糖代謝通路[7]。例如，lncRNA中CDKN2B反義RNA 1 (CDKN2B antisense RNA 1，ANRIL)上調(diào)葡萄糖轉(zhuǎn)運(yùn)蛋白1 (glucose transporter 1, GLUT1)和乳酸脫氫酶A(Lactate dehydrogenase A, LDHA)的表達(dá)，從而增加葡萄糖攝取并促進(jìn)鼻咽癌的發(fā)展[8]，LINC00092直接與6-磷酸果糖-2-激酶/果糖-2，6-二磷酸酶2(6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 2，PFKFB2)結(jié)合能增強(qiáng)糖酵解并促進(jìn)腫瘤的發(fā)生和發(fā)展[9]。在膀胱癌中，lncRNA中過表達(dá)的尿路上皮癌胚抗原1 (urothelial cancer associated 1，UCA1)通過上調(diào)己糖激酶2 (Hexokinase 2, HK2) 促進(jìn)糖酵解[10]。然而，參與SCCHN糖酵解重編程的lncRNA對SCCHN發(fā)生發(fā)展的作用仍不清楚。因此，本研究從TCGA數(shù)據(jù)集中獲得了具有顯著預(yù)后價(jià)值的lncRNA。并基于這些GRLs構(gòu)建了預(yù)后風(fēng)險(xiǎn)模型，并確定SCCHN高/低風(fēng)險(xiǎn)組在信號(hào)通路、免疫微環(huán)境、免疫檢查點(diǎn)、m6A相關(guān)基因方面的差異。

1 材料與方法

1.1 材料

從TCGA數(shù)據(jù)庫(https://cancergenome.nih.gov/)下載SCCHN患者的轉(zhuǎn)錄組數(shù)據(jù)和臨床數(shù)據(jù)[11]，轉(zhuǎn)錄組數(shù)據(jù)包含44例正常組織樣本和502例頭頸癌組織樣，臨床數(shù)據(jù)包含501例頭頸癌患者的生存時(shí)間和生存狀態(tài)。

1.2 方法

1.2.1數(shù)據(jù)預(yù)處理 利用Perl語言腳本整合SCCHN患者的轉(zhuǎn)錄組數(shù)據(jù)，獲得Ensemble表達(dá)矩陣，同時(shí)，從Ensembl數(shù)據(jù)庫(https://asia.ensembl.org/index.html)下載人類基因注釋文件,將Ensemble表達(dá)矩陣轉(zhuǎn)化為Gene Symbol表達(dá)矩陣，根據(jù)Gene Symbol表達(dá)矩陣與GENCODE數(shù)據(jù)庫中l(wèi)ncRNA染色體的位置進(jìn)行比較識(shí)別lncRNAs[12]，獲得SCCHN的lncRNA表達(dá)矩陣。從MSigDB數(shù)據(jù)庫中獲得了總共274個(gè)糖酵解基因,利用Perl腳本從Gene Symbol中提取糖酵解基因表達(dá)矩陣。

1.2.2差異表達(dá)糖酵解基因和lncRNA篩選 使用R軟件包limma篩選差異的糖酵解基因和lncRNA,篩選標(biāo)準(zhǔn)：|log2(fold-change)|>1以及矯正后的P(false discovery rate,fdr)<0.05，對于糖酵解基因表達(dá)矩陣和lncRNA表達(dá)矩陣篩選出的差異基因和差異lncRNA分別，利用R包ggplots繪制正常樣本和腫瘤樣本的火山圖。

1.2.3差異糖酵解基因的富集分析 為了探索潛在糖酵解基因與SCCHN生物通路之間的關(guān)系和相關(guān)程度，使用 DAVID對差異的糖酵解基因進(jìn)行GO-BP富集和KEGG通路分析，設(shè)定P<0.05為顯著性基因富集[13]。

1.2.4糖酵解相關(guān)lncRNA 在R軟件中，本研究使用糖酵解基因lncRNA共表達(dá)方法分析糖酵解基因和lncRNA的相關(guān)性，相關(guān)系數(shù)Cor>0.4且P<0.05的lncRNA被認(rèn)為是糖酵解相關(guān)lncRNA。

1.2.5糖酵解相關(guān)lncRNA模型構(gòu)建 通過單因素Cox回歸(P<0.05)篩選預(yù)后糖酵解相關(guān)lncRNA后，通過多因素Cox回歸確定關(guān)鍵的lncRNA[14]?；谝韵鹿綐?gòu)建SCCHN患者的風(fēng)險(xiǎn)評分模型，計(jì)算每個(gè)SCCHN患者的風(fēng)險(xiǎn)得分。風(fēng)險(xiǎn)評分=∑βlncRNA×ExplncRNA(βlncRNA表示多因素Cox回歸分析中計(jì)算出的每個(gè)lncRNA的回歸系數(shù)，Exp lncRNA代表樣本中每個(gè)lncRNA的表達(dá)值。)

1.2.6生存分析 以SCCHN患者的中位風(fēng)險(xiǎn)得分作為臨界值，進(jìn)一步將SCCHN患者分為高危組和低危組(低危組得分<中位，高危組得分≥中位)。采用R包中的“Survival”繪制Kaplan-Meier生存曲線；使用R包“timeROC”包繪制ROC曲線。此外，通過校準(zhǔn)曲線構(gòu)建和列線圖，以確定作為獨(dú)立預(yù)后因素風(fēng)險(xiǎn)模型的準(zhǔn)確性。

1.2.7模型的特點(diǎn)和應(yīng)用 免疫微環(huán)境與SCCHN的發(fā)生發(fā)展密切相關(guān)。根據(jù)SCCHN患者的表達(dá)數(shù)據(jù)，通過 ESTIMATE 估計(jì)免疫和基質(zhì)評分，以代表基質(zhì)和免疫細(xì)胞的存在[21]；此外，使用兩種算法 CIBERSORT、MCPcounter來比較不同風(fēng)險(xiǎn)組中各種免疫細(xì)胞比例的差異[15]；同時(shí)提取免疫檢查點(diǎn)基因、m6A相關(guān)基因的表達(dá)水平，通過Wilcox檢驗(yàn)分析它們在不同風(fēng)險(xiǎn)組中的表達(dá)差異。

1.3 統(tǒng)計(jì)學(xué)分析

R軟件(v3.6.3)用于統(tǒng)計(jì)分析，Wilcox檢驗(yàn)用于組間比較，Pearson用于分析糖酵解基因與lncRNA之間的相關(guān)性，單因素和多因素Cox回歸分析影響SCCHN患者總生存期的相關(guān)因素；P<0.05為差異有統(tǒng)計(jì)意義。

2 結(jié)果

2.1 差異和富集分析

根據(jù)所設(shè)定的顯著性閾值，從274個(gè)糖酵解基因中篩選得到76個(gè)糖酵解差異表達(dá)基因(圖1A)，其中上調(diào)的有52個(gè)，下調(diào)的有24個(gè)。從lncRNA表達(dá)矩陣中篩選得到965個(gè)差異的lncRNA(圖1B)，其中上調(diào)的有864個(gè)，下調(diào)的有101個(gè)。此外， 76 個(gè)糖酵解差異表達(dá)基因富集到GO和KEGG 通路上(圖1C～D)。結(jié)果顯示，大多數(shù)差異的糖酵解基因都富集于代謝路徑。鑒定出的糖酵解基因與腫瘤發(fā)生發(fā)展中的幾個(gè)重要生物學(xué)過程有關(guān)，如缺氧通路、代謝通路、AMPK信號(hào)通路、HIF-1信號(hào)通路，進(jìn)一步證明糖酵解與腫瘤缺氧微環(huán)境密切相關(guān)。

注：A為差異表達(dá)的糖酵解相關(guān)mRNA的火山圖，B為差異表達(dá)的糖酵解相關(guān)lncRNA的火山圖；紅色圓為上調(diào)，綠色圓為下調(diào)；C為差異的糖酵解基因GO分析結(jié)果、D為差異的糖酵解相關(guān)lncRNA的KEGG分析結(jié)果，色標(biāo)表示P值，直方圖大小表示計(jì)數(shù)。圖1 糖酵解基因的差異和富集分析Fig.1 The differences and enrichment analyses of glycolysis-related genes

2.2 基于CRLs的風(fēng)險(xiǎn)模型構(gòu)建與驗(yàn)證

對Pearson相關(guān)分析 (Cor>0.4,P<0.05)獲得的298個(gè)GRLs進(jìn)行單因素Cox回歸分析，結(jié)果顯示有60個(gè)GRLs與SCCHN患者的預(yù)后顯著相關(guān)(圖2A)。本研究進(jìn)一步對這60個(gè)lncRNA進(jìn)行多因素Cox回歸分析，確定了25個(gè)關(guān)鍵的lncRNA用于構(gòu)建SCCHN患者的預(yù)后風(fēng)險(xiǎn)模型。生存曲線展示了高風(fēng)險(xiǎn)的SCCHN患者預(yù)后較差(圖2B)，高低風(fēng)險(xiǎn)組患者的存活狀態(tài)、總生存期及l(fā)ncRNA的表達(dá)均具有顯著差異(圖2C～E)。

注：A為預(yù)后相關(guān)lncRNA森林圖，B為生存曲線，C為高/低風(fēng)險(xiǎn)SCCHN患者的RS分布圖，D為SCCHN患者的生存時(shí)間散點(diǎn)圖，E為SCCHN患者的生存熱圖。圖2 lncRNA與SCCHN患者的生存和風(fēng)險(xiǎn)分析Fig.2 Survival and risk analysis of lncRNA in SCCHN patients

由臨床特征和風(fēng)險(xiǎn)模型的單因素和多因素Cox 回歸分析結(jié)果可知，風(fēng)險(xiǎn)評分P<0.001。因此，風(fēng)險(xiǎn)評分可作為SCCHN患者的獨(dú)立預(yù)后因素(圖3A～B)。同時(shí)本研究通過繪制ROC曲線來評估風(fēng)險(xiǎn)模型的可靠性，結(jié)果表明風(fēng)險(xiǎn)評分的曲線下面積(AUC)值高于年齡、性別、臨床分期和組織學(xué)分級(jí)，表明風(fēng)險(xiǎn)評分對SCCHN患者預(yù)后的預(yù)測能力要優(yōu)于傳統(tǒng)的臨床參數(shù)(圖3C～D)。

注：A為風(fēng)險(xiǎn)模型的單因素分析，B為風(fēng)險(xiǎn)模型的多因素分析，C為 SCCHN患者1、2、3年生存期的ROC曲線，D為各臨床指標(biāo)與RS的ROC曲線比較。圖3 臨床特征和風(fēng)險(xiǎn)模型的Cox和ROC分析Fig.3 Cox analysis of clinical characteristics and ROC analysis of the risk model

基于臨床特征和風(fēng)險(xiǎn)模型進(jìn)一步構(gòu)建決策曲線DCA(圖4A)和列線圖(圖4B)。DCA曲線和列線圖的結(jié)果都證明了構(gòu)建模型的準(zhǔn)確性。糖酵解相關(guān)lncRNAs表達(dá)與臨床特征的熱圖(圖4C)表明，糖酵解相關(guān)基因與預(yù)后具有密切關(guān)系，結(jié)果與構(gòu)建模型預(yù)測的結(jié)果吻合。

注：A為不同臨床因素繪制的決策曲線圖，B為不同臨床因素預(yù)后預(yù)測列線圖，C為25種糖酵解表達(dá)水平與臨床特征之間關(guān)聯(lián)的熱圖。圖4 風(fēng)險(xiǎn)模型的驗(yàn)證Fig.4 Validating the risk model

2.3 免疫細(xì)胞和免疫檢查點(diǎn)分析

免疫檢查點(diǎn)抑制劑可阻斷表達(dá)免疫檢查點(diǎn)腫瘤細(xì)胞對免疫細(xì)胞的抑制作用。鑒于此，為了進(jìn)一步探索風(fēng)險(xiǎn)模型的臨床應(yīng)用，對相關(guān)的免疫細(xì)胞(圖5A)和免疫功能(圖5B)進(jìn)行了分析，結(jié)果表明高/低風(fēng)險(xiǎn)組的免疫細(xì)胞浸潤類型有明顯差異，這為SCCHN的免疫治療提供了理論依據(jù)。本研究比較了兩個(gè)風(fēng)險(xiǎn)組之間免疫檢查點(diǎn)基因的差異(圖5C)可知，低風(fēng)險(xiǎn)組中免疫檢查點(diǎn)基因表達(dá)量較高，說明RS較低的SCCHN患者對免疫檢查點(diǎn)阻斷(ICB)治療更敏感。還比較了兩個(gè)風(fēng)險(xiǎn)組之間m6A的表達(dá)水平，共匹配了12個(gè)m6A 調(diào)節(jié)因子，除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個(gè) m6A 調(diào)節(jié)因子的表達(dá)水平在低風(fēng)險(xiǎn)組中較高 (圖5D)。

注：A為高低風(fēng)險(xiǎn)組免疫細(xì)胞熱圖，B為高低風(fēng)險(xiǎn)組免疫功能差異分析，C為高低風(fēng)險(xiǎn)組免疫檢查點(diǎn)差異分析，D為高低風(fēng)險(xiǎn)組m6A差異分析。圖5 高低風(fēng)險(xiǎn)組患者m6A免疫和基因分析Fig.5 Immunity and m6A analyses in high risk group

3 討論

糖酵解是腫瘤細(xì)胞代謝的主要方式，可以為腫瘤細(xì)胞增殖提供所需的能量[16-17]，其最終產(chǎn)物乳酸被釋放到細(xì)胞外，可維持酸性腫瘤微環(huán)境并促進(jìn)毛細(xì)血管的形成，最終導(dǎo)致腫瘤細(xì)胞的增殖、耐藥、侵襲和轉(zhuǎn)移[18]。已有研究表明，糖酵解基因的上調(diào)與SCCHN的增殖、轉(zhuǎn)移以及免疫逃逸有著密切關(guān)系[19-20]。然而，靶向SCCHN糖酵解途徑的藥物的作用機(jī)制尚不清楚，因此對其分子機(jī)制的深入研究仍具有重要意義。

已有研究表明，lncRNA在腫瘤的發(fā)生和發(fā)展中發(fā)揮著重要作用，且在肝細(xì)胞癌、結(jié)直腸癌和胃癌等多種癌癥中，lncRNA可通過重編程葡萄糖代謝和糖酵解影響腫瘤的發(fā)生發(fā)展[21]。例如，lncRNA AGPG(actin gamma 1 pseudogene 25) 通過抑制 Lys302的泛素化和隨后的蛋白酶體依賴性PFKFB3(6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3) 降解來增加 PFKFB3 的穩(wěn)定性,并激活糖酵解途徑，導(dǎo)致食管癌細(xì)胞的代謝重編程，促進(jìn)食管癌細(xì)胞的遷移[22]。然而，在SCCHN中對糖酵解相關(guān)的lncRNAs的研究仍然很少。

因此，本研究通過Pearson分析、單因素Cox回歸發(fā)現(xiàn)了GRLs與SCCHN患者預(yù)后顯著相關(guān)，并建立了由25個(gè)關(guān)鍵的GRLs組成的風(fēng)險(xiǎn)模型。根據(jù)功能富集分析，本研究發(fā)現(xiàn)鑒定出的GRLs與HIF-1信號(hào)通路、氨基酸代謝和核苷酸代謝等有關(guān)。因此，推測上述lncRNA可能通過糖酵解誘導(dǎo)腫瘤微環(huán)境的缺氧、腫瘤細(xì)胞增殖和轉(zhuǎn)移等。此外，本研究還通過ROC曲線對風(fēng)險(xiǎn)模型的可靠性進(jìn)行了評估，發(fā)現(xiàn)AUC值分別為0.740、0.730及0.744，說明本研究的模型具有較好的預(yù)測性。之后本研究進(jìn)一步探索了風(fēng)險(xiǎn)模型的臨床意義，發(fā)現(xiàn)兩個(gè)風(fēng)險(xiǎn)組表現(xiàn)出不同的免疫細(xì)胞浸潤、免疫功能、免疫檢查點(diǎn)、m6A調(diào)節(jié)因子水平，且高風(fēng)險(xiǎn)組與不良的臨床表現(xiàn)和免疫功能低下顯著相關(guān)。同時(shí)，Shen等[23]發(fā)現(xiàn)，腫瘤糖酵解誘導(dǎo)的缺氧和酸性微環(huán)境可導(dǎo)致代謝介導(dǎo)的T細(xì)胞功能障礙，Jung等[24]發(fā)現(xiàn)腫瘤的糖酵解可誘導(dǎo)腫瘤免疫抑制和免疫逃逸，這與本研究通過富集分析以及免疫功能差異分析得到的結(jié)果一致。因此，靶向代謝的腫瘤免疫治療策略有望提高免疫治療的療效，本研究也希望通過進(jìn)一步的研究探索糖酵解(通過m6A修飾或免疫檢查點(diǎn))的分子機(jī)制，以提高免疫治療的療效。

綜上所述，本研究中篩選出的糖酵解相關(guān)的lncRNAs在 SCCHN的發(fā)生和進(jìn)展中起著重要作用，基于糖酵解相關(guān)lncRNAs的風(fēng)險(xiǎn)模型是SCCHN患者的獨(dú)立預(yù)后因素。所構(gòu)建的基于糖酵解相關(guān)lncRNAs的風(fēng)險(xiǎn)模型是SCCHN患者的獨(dú)立預(yù)后因素。通過綜合分析，糖酵解相關(guān)的lncRNAs的模型為SCCHN的臨床治療提供了一定的指導(dǎo)建議。但本研究仍存在一些局限性。首先，由于可以注釋lncRNA表達(dá)的SCCHN樣本數(shù)量有限，需要更多具有同源信息的患者納入研究并證明本研究的可信度。其次，本研究僅通過生物信息學(xué)分析探討了lncRNA的功能，需要實(shí)驗(yàn)數(shù)據(jù)來支持這些結(jié)論。盡管有這些限制，本研究還是使用了ROC和列線圖來證明風(fēng)險(xiǎn)模型預(yù)測的有效性。