趙旭輝,黃小敏,達(dá)德轉(zhuǎn),許 焱,崔曉東,李紅玲

(1.甘肅中醫(yī)藥大學(xué)第一臨床醫(yī)學(xué)院,甘肅 蘭州 730000;2.甘肅省人民醫(yī)院 腫瘤內(nèi)科,甘肅 蘭州 730000)

胃癌(gasrtic cancer, GC)是消化系統(tǒng)最常見的惡性腫瘤之一,為世界第四大常見癌癥,居癌癥相關(guān)死亡原因的第2位[1]。GC發(fā)病隱匿導(dǎo)致診斷時間較晚且預(yù)后不良,平均5年生存率低于20%[2]。目前手術(shù)、放療、化療臨床綜合治療效果不理想,含鉑類藥物的一線化療治療晚期GC已達(dá)到療效瓶頸[3-4]。雖然分子靶向治療以及PD-1/PD-L1免疫檢查點(diǎn)抑制劑對晚期GC有一定的療效,開辟了新的治療途徑。然而,一些GC患者的療效仍然較差[5-6]。20世紀(jì)初德國著名的生物化學(xué)家Otto Warburg首次從細(xì)胞代謝角度提出了關(guān)于腫瘤的“Warburg效應(yīng)”,即在腫瘤細(xì)胞中當(dāng)線粒體功能失調(diào)后,主要通過增強(qiáng)有氧酵解的方式來獲取充足的能量[7]。越來越多的證據(jù)表明,異常的能量代謝方式已成為腫瘤細(xì)胞生長的特殊模式,葡萄糖的異常消耗與乳酸的大量堆積使得腫瘤細(xì)胞更具快速生長、轉(zhuǎn)移的優(yōu)勢[8]。同時此過程也會使細(xì)胞外pH值明顯下降造成微環(huán)境的酸化,微環(huán)境的酸化可破壞鄰近的正常細(xì)胞且造成細(xì)胞外基質(zhì)降解,從而促進(jìn)腫瘤細(xì)胞的入侵[9]。有學(xué)者指出這種特殊的產(chǎn)能方式以及微環(huán)境的變化特征有利于腫瘤細(xì)胞逃避機(jī)體的凋亡機(jī)制,導(dǎo)致患者不良預(yù)后[10]。有文獻(xiàn)報道,糖酵解相關(guān)基因可預(yù)測癌癥患者的預(yù)后,包括透明細(xì)胞腎細(xì)胞癌[11]、肺腺癌[12]、肝細(xì)胞癌[13]、乳腺癌[14-15]、卵巢癌[16]和結(jié)直腸癌[17]。研究表明,糖酵解相關(guān)基因可用于有效評估GC患者預(yù)后[18-20]。然而,目前仍然缺乏關(guān)于糖酵解相關(guān)基因與GC患者預(yù)后關(guān)系的系統(tǒng)研究?；诖?本文通過利用生物信息學(xué)開發(fā)糖酵解相關(guān)基因以預(yù)測GC患者預(yù)后,為臨床提供參考。

1 資料與方法

1.1患者的臨床信息和信使核糖核酸(messenger RNA,mRNA)表達(dá)數(shù)據(jù)集 我們從癌癥基因組圖譜(the cancer genome atlas, TCGA)數(shù)據(jù)庫(https://cancergenome.nih.gov/)中提取了GC患者的臨床資料和mRNA表達(dá)譜。納入445例患者的臨床信息,并統(tǒng)計(jì)了匹配的年齡,分期,分級,將樣本的基因轉(zhuǎn)錄數(shù)據(jù)統(tǒng)一轉(zhuǎn)化為log2(FPKM 值+1)后進(jìn)行數(shù)據(jù)集序列整合。

1.2方法

1.2.1數(shù)據(jù)處理和風(fēng)險評分計(jì)算 在分子特征數(shù)據(jù)庫v4.3.2檢索所有糖酵解相關(guān)基因集,包括5個不同的基因集(BIOCARTA_GLYCOLYSIS_PATHWAY、GO_GLYCOLYTIC_PROCESS、HALLMARK_GLYCOLYSIS、KEGG_ GLYCOLYSIS_GLUCONEOGENESIS、 REACTOME _GLYCOLYSIS),查閱相關(guān)文獻(xiàn)后共篩選出371個糖酵解相關(guān)基因。接下來,我們采用單變量Cox回歸分析來確定與總生存期(overall survival,OS)相關(guān)基因。我們應(yīng)用最小絕對收縮和選擇算子(least absolute shrinkage and selection operator, LASSO)回歸分析來確定具有非零系數(shù)的潛在預(yù)測因子。然后,進(jìn)行多變量Cox回歸分析以確定最終建模中涉及的樞紐基因并獲得相關(guān)數(shù)。

1.2.2預(yù)測模型的構(gòu)建與驗(yàn)證 以中位風(fēng)險評分為臨界值,將445例GC患者分為兩個亞組(高風(fēng)險組和低風(fēng)險組),并繪制Kaplan-Meier生存分析圖以及無進(jìn)展生存期(progression-free survival,PFS)曲線,以比較患者生存率與log-rank檢驗(yàn)的差異。此外,在R軟件(版本 4.0.2)中生成1年、3年和5年的受試者工作特征(receiver operating characteristic, ROC)曲線,以及臨床相關(guān)特征的ROC 曲線,以計(jì)算每個預(yù)測因子模型的曲線下面積(area under the curve, AUC),以進(jìn)一步評估模型的效率和準(zhǔn)確性。P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

1.2.3列線圖的構(gòu)建和評估 通過R軟件(版本 4.0.2)構(gòu)建列線圖,并采用列線圖模型比較GC患者的生存概率。生成校準(zhǔn)圖和C指數(shù)作為列線圖性能的評估。

1.2.4高風(fēng)險組和低風(fēng)險組間的基因集變異分析(gene set variation analysis,GSVA) 為了進(jìn)一步探索模型中高風(fēng)險組和低風(fēng)險組間生物途徑激活狀態(tài)的差異,使用“GSVA”R包進(jìn)行GSVA。從MSigDB數(shù)據(jù)庫下載“c2.cp.kegg.v7.4.symbols”基因集用于GSVA。GSVA通常用于以非參數(shù)和無監(jiān)督方法估計(jì)表達(dá)數(shù)據(jù)集樣品中途徑和生物過程活性的變化。

1.2.5免疫細(xì)胞浸潤分析 通過R語言對兩組樣本中免疫細(xì)胞的突變、免疫細(xì)胞浸潤、免疫評分等差異進(jìn)行了分析。使用“maftools”R包評估GC樣本中22個免疫細(xì)胞的突變。GC樣本中22個免疫細(xì)胞的浸潤水平基于CIBERSORT(http://cibersort.stanford.edu/)進(jìn)行評估。使用R包評估GC樣本中22個免疫細(xì)胞的免疫和基質(zhì)評分。

1.2.6統(tǒng)計(jì)學(xué)方法 使用Excel軟件和R軟件進(jìn)行統(tǒng)計(jì)分析。通過單變量和多變量Cox回歸分析評估單個指標(biāo)的預(yù)后意義。采用Kaplan-Meier 生存曲線和log-rank檢驗(yàn)評估預(yù)后結(jié)局。采用非配對t檢驗(yàn)比較兩組不同表達(dá)水平。R軟件用于繪制熱圖,ROC曲線,富集,森林和校準(zhǔn)圖。P<0.05表示差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1糖酵解相關(guān)基因預(yù)后特征的構(gòu)建和驗(yàn)證 使用從TCGA獲得的445個GC樣本來構(gòu)建模型并采用log-rank檢驗(yàn)方法和利用“survival”R 包對樣本基因進(jìn)行單變量Cox回歸分析探究糖酵解相關(guān)差異基因與患者OS的關(guān)系。單因素Cox回歸分析顯示,19個糖酵解相關(guān)差異基因與患者的預(yù)后相關(guān)。接下來,選擇上述的19個糖酵解相關(guān)基因進(jìn)行LASSO回歸分析,得到15個糖酵解相關(guān)基因(見圖1及表1)。基于LASSO回歸分析的結(jié)果,使用15個糖酵解相關(guān)基因(PFKFB2、UHRF1、ACYP1、CLDN9、STC1、EFNA3、NUP50、、ADH4、ANGPTL4、PKP2、VCAN、HIF1A、LHX9、ANKZF1、ALDH3A2)建立和驗(yàn)證預(yù)測GC患者結(jié)局的風(fēng)險模型。

圖1 糖酵解相關(guān)差異基因與患者OS的單變量Cox回歸分析

表1 通過LASSO回歸分析得到與GC患者預(yù)后相關(guān)的15個糖酵解相關(guān)基因

采用log-rank檢驗(yàn)方法,根據(jù)風(fēng)險評分的中值將TCGA中的GC患者分為高風(fēng)險組和低風(fēng)險組,兩組OS和FPS差異有統(tǒng)計(jì)學(xué)意義(均P<0.01)。ROC曲線表明,1年、3年、5年OS的AUC為分別為0.654、0.729、0.799,表明上述15個糖酵解相關(guān)基因標(biāo)記在評估GC患者預(yù)后結(jié)局方面具有良好的特異度和敏感度。見圖2～4。

圖2 用于構(gòu)建最終預(yù)測模型的 LASSO 回歸分析 a.選擇LASSO模型中的最優(yōu)參數(shù)(λ);使用最小標(biāo)準(zhǔn)在最佳值繪制垂直虛線; b. LASSO 分布的部分似然偏差

圖3 高風(fēng)險組和低風(fēng)險組的OS和PFS

圖4 1年、3年和5年OS的ROC曲線

2.2整合糖酵解相關(guān)基因標(biāo)記的列線圖模型構(gòu)建 構(gòu)建基于基因特征的風(fēng)險評分與臨床病理特征(分期、性別、分級、風(fēng)險和年齡)相結(jié)合的列線圖模型,用以評估GC患者的生存可能性(圖5a)。列線圖證明了這15個基因標(biāo)記在個體化生存估計(jì)中的價值。在校準(zhǔn)、辨別和臨床可用性方面評估列線圖的性能,C指數(shù)為0.8557(圖5b)。ROC曲線中風(fēng)險(AUC=0.775)和列線圖(AUC=0.753)具有較好的預(yù)測價值(圖5c)?；诹芯€圖單因素和多因素分析,發(fā)現(xiàn)在單因素分析中年齡[HR=1.024,95%置信區(qū)間(CI):1.006～1.042,P<0.05]、分級[HR=1.549,95%CI:1.244～1.929,P<0.01]和列線圖[HR=1.049,95%CI:1.036～1.061,P<0.01]與預(yù)后結(jié)局相關(guān)。多因素Cox分析表明,列線圖[HR=1.039,95%CI:1.020～1.059,P<0.01]可作為獨(dú)立的預(yù)后因子。見圖5。

圖5 基于糖酵解和臨床因素的列線圖 a.包含年齡和風(fēng)險評分的列線圖,以預(yù)測1年、3年和5年的OS; b.列線圖的校準(zhǔn)圖; c.風(fēng)險、列線圖和其他臨床特征的ROC曲線

2.3基于15個基因特征的風(fēng)險評分作為獨(dú)立的預(yù)后因素 采用單因素和多因素分析,分析各臨床病理特征的預(yù)后價值。單因素分析結(jié)果表明,年齡[HR=1.026,95%置信區(qū)間(CI):1.008～1.044,P<0.01],分級[HR=1.534,95%CI:1.241～1.896,P<0.01]和風(fēng)險系數(shù)[HR=3.690,95%CI:2.468～5.518,P<0.01]與預(yù)后結(jié)局相關(guān)(見圖6a)。隨后,多因素Cox分析表明,年齡[HR=1.038,95%置信區(qū)間(CI):1.019～1.057,P<0.01],分級[HR=1.549,95%CI:1.238～1.937,P<0.01]和風(fēng)險系數(shù)[HR=3.686,95%CI:2.390～5.686,P<0.01]是獨(dú)立的預(yù)后指標(biāo)(見圖6b)。同樣,通過繪制ROC曲線發(fā)現(xiàn)年齡(AUC=0.606),分級(AUC=0.606), 風(fēng)險系數(shù)(0.799)具有很好的預(yù)測價值(見圖6c)。

圖6 臨床特征的單因素、多因素分析及ROC曲線 a.單變量分析; b.多變量分析; c.ROC曲線

2.4基于糖酵解相關(guān)基因的臨床相關(guān)性分析 臨床相關(guān)性分析發(fā)現(xiàn),高低風(fēng)險組年齡和性別差異無統(tǒng)計(jì)學(xué)意義。G2與G3患者分級差異有統(tǒng)計(jì)學(xué)意義。在TNM分期中,T1期患者與T2、T3、T4期患者差異有統(tǒng)計(jì)學(xué)意義(P<0.01),T1期之后的患者風(fēng)險明顯上升。 M0與M1期患者風(fēng)險差異具有統(tǒng)計(jì)學(xué)意義,N0期與N1、N3期患者風(fēng)險差異具有統(tǒng)計(jì)學(xué)意義,而N0與N2期患者風(fēng)險差異無統(tǒng)計(jì)學(xué)意義。見圖7。

圖7 基于糖酵解相關(guān)基因的臨床相關(guān)性分析 a.年齡; b.性別; c.分級; d.分期; e-g.TNM分期

2.5糖酵解相關(guān)基因的表達(dá)水平與免疫細(xì)胞浸潤的分析 不同水平的免疫浸潤可能在一定程度上解釋了為什么具有相同組織學(xué)類型的癌癥患者可能具有不同的臨床結(jié)果[20]。我們分析糖酵解相關(guān)基因與免疫浸潤水平的相關(guān)性和免疫相關(guān)功能的相關(guān)性,旨在揭示15個糖酵解相關(guān)基因標(biāo)記影響GC患者預(yù)后的潛在機(jī)制。箱線圖分析顯示,高風(fēng)險組細(xì)胞浸潤水平(巨噬細(xì)胞 M1)高于低風(fēng)險組(P<0.05),而高風(fēng)險組激活的CD4+T細(xì)胞、輔助性濾泡T細(xì)胞、單核細(xì)胞(P<0.05)水平低于低風(fēng)險組(見圖8a)。關(guān)于免疫相關(guān)功能的分析顯示,高風(fēng)險組CCR、Parainflammation、Type_I_IFN_Reponse、Type_II_ IFN_Reponse功能高于低風(fēng)險組(均P<0.05),而高風(fēng)險組MHC_class_I功能低于低風(fēng)險組(見圖8b)。

圖8 箱線圖分析 a.兩組免疫細(xì)胞浸潤分布; b.免疫相關(guān)功能

2.6GSVA分析確定的高低風(fēng)險亞組差異生物學(xué)行為 進(jìn)一步對上述糖酵解相關(guān)基因進(jìn)行GSVA分析,以探索高低風(fēng)險組間生物學(xué)通路激活狀態(tài)的差異。利用熱圖對相關(guān)生物化過程進(jìn)行可視化分析并根據(jù)P值大小進(jìn)行排列。結(jié)果表明,高風(fēng)險組致癌激活途徑顯著富集,ECM_RECEPTOR_INTERACTION、FOCAL_ADHESION、CHEMOKINE_SIGNALING_PATHWAY、TGF_BETA_SIGNALING_PATHWAY等。這些通路可以調(diào)控腫瘤發(fā)生發(fā)展,因而可能為高危人群預(yù)后不良提供一定的理論依據(jù)。見圖9。

圖9 熱圖 注:紅色表示激活途徑,藍(lán)色表示抑制途徑

3 討 論

近年來,關(guān)于代謝組學(xué)的研究越來越受到關(guān)注,人類惡性腫瘤與能量代謝有著密不可分的聯(lián)系[21]。在氧氣充足的條件下,腫瘤細(xì)胞也會通過無氧糖酵解的方式給腫瘤細(xì)胞提供更多的能量,使腫瘤細(xì)胞得以快速生長(即“Warburg效應(yīng)”)[22]。關(guān)于腫瘤細(xì)胞與糖酵解的研究近年來層出不窮。Hu等[23]發(fā)現(xiàn),在胰腺癌中,UHRF1可通過抑制SIRT4促進(jìn)無氧糖酵和腫瘤增殖。Zhao等[24]發(fā)現(xiàn)LncRNA MIR17HG可以通過介導(dǎo)糖酵解相關(guān)正反饋回路促進(jìn)結(jié)直腸癌肝轉(zhuǎn)移。還有研究發(fā)現(xiàn),EB病毒編碼的miRNAs也可以調(diào)節(jié)細(xì)胞的能量代謝,增加細(xì)胞攝取葡萄糖的能力并促進(jìn)細(xì)胞內(nèi)糖酵解速率。Cai等[25]還發(fā)現(xiàn)EBV-miR-BART1可依賴糖酵解并通過PTEN影響 HIF1α從而誘導(dǎo)鼻咽癌細(xì)胞的轉(zhuǎn)移與擴(kuò)散。然而,年齡和轉(zhuǎn)移診斷等臨床病理學(xué)特征不足以準(zhǔn)確預(yù)測癌癥患者預(yù)后。因此,越來越多的mRNA被確定為腫瘤進(jìn)展或預(yù)后的生物標(biāo)志物,并且已經(jīng)評估了生物標(biāo)志物的臨床意義[26-27]。本研究中,在TCGA中獲得了445例GC患者的臨床資料和mRNA表達(dá)譜,然后在分子特征數(shù)據(jù)庫v4.0上檢索了5個不同糖酵解相關(guān)基因集,查閱相關(guān)文獻(xiàn)后共篩選出371個糖酵解相關(guān)基因。接下來,我們采用單變量Cox回歸分析來得到與預(yù)后相關(guān)的19個糖酵解相關(guān)的差異基因,接著對這些基因進(jìn)行LASSO回歸分析,基于LASSO回歸分析的結(jié)果,使用15個糖酵解相關(guān)基因建立和驗(yàn)證預(yù)測GC患者結(jié)局和系數(shù)的風(fēng)險模型。

其次,采用log-rank檢驗(yàn)方法,根據(jù)風(fēng)險評分的中值將TCGA-STAD中的GC患者分為高風(fēng)險組和低風(fēng)險組,進(jìn)一步通過 Kaplan-Meier 曲線、ROC 曲線、單因素及多因素 Cox 回歸和通過繪制列線圖分析考察風(fēng)險評分模型預(yù)測能力,結(jié)果顯示高風(fēng)險組病死率高于低風(fēng)險組,ROC曲線提示該模型具有較好的短期預(yù)測能力。通過GSVA分析發(fā)現(xiàn),高危人群的致癌激活通路(ECM_RECEPTOR_INTERACTION、TGF_BETA_SIGNALING_PATHWAY、FOCAL_ADHESION等)顯著富集,為高危人群預(yù)后不良提供了更全面、更有說服力的解釋。最后,將風(fēng)險和年齡確定為獨(dú)立的OS相關(guān)變量,并將其納入列線圖,結(jié)果表明列線圖可被視為GC患者臨床診斷和治療的有效工具。這些結(jié)果表明,基因標(biāo)記具有較強(qiáng)預(yù)測GC患者預(yù)后的能力,對臨床治療決策具有一定的指導(dǎo)意義。此外通過CIBERSORT探討了糖酵解相關(guān)基因與免疫浸潤水平的相關(guān)性和免疫相關(guān)功能的相關(guān)性,揭示了15個糖酵解相關(guān)基因標(biāo)記影響GC預(yù)后的潛在機(jī)制。我們發(fā)現(xiàn)高風(fēng)險組中細(xì)胞浸潤水平(巨噬細(xì)胞 M1)高于低風(fēng)險組,而高風(fēng)險組激活的CD4+T細(xì)胞、輔助性濾泡T細(xì)胞、單核細(xì)胞(均P<0.05)水平低于低風(fēng)險組。顯然,高風(fēng)險組反映了一個免疫抑制的腫瘤微環(huán)境,相較于低風(fēng)險組缺乏免疫細(xì)胞,與高風(fēng)險組的預(yù)后不良相吻合。關(guān)于免疫相關(guān)功能的分析顯示,高風(fēng)險組CCR、Parainflammation、Type_I_IFN_Reponse、Type_II_ IFN_Reponse功能高于低風(fēng)險組,而高風(fēng)險組MHC_class_I低于低風(fēng)險組。以上說明這15個糖酵解相關(guān)基因標(biāo)記與GC患者的免疫細(xì)胞浸潤與免疫相關(guān)功能密切相關(guān)。

綜上,建立的15個糖酵解相關(guān)基因標(biāo)記(PFKFB2、UHRF1、ACYP1、CLDN9、STC1、EFNA3、NUP50、ADH4、ANGPTL4、PKP2、VCAN、HIF1A、LHX9、、ANKZF1和ALDH3A2)可作為預(yù)測GC患者預(yù)后的可靠工具,可能為GC的潛在發(fā)病機(jī)制提供新的見解。