基于WGCNA分析及風(fēng)險(xiǎn)評(píng)分構(gòu)建頭頸部鱗狀細(xì)胞癌的預(yù)后模型

金祥

頭頸部鱗狀細(xì)胞癌(head and neck squamous cell carcinoma,HNSCC)多發(fā)于口腔、鼻、咽、喉的黏膜上皮細(xì)胞,是全世界第六大常見的惡性腫瘤[1]。危險(xiǎn)因素包括吸煙、酗酒、環(huán)境致癌物質(zhì)及HPV感染等。大多數(shù)HNSCC的患者發(fā)現(xiàn)確診時(shí)已處于晚期[2]。盡管手術(shù)技術(shù)、放化療及免疫治療的應(yīng)用顯示良好的效果,但五年生存率僅有50%[3]。其中一個(gè)重要的原因是缺乏針對(duì)早期HNSCC的有效生物標(biāo)志物的認(rèn)識(shí)與應(yīng)用以便患者及時(shí)醫(yī)治。近年來,靶向治療已多應(yīng)用于臨床治療HNSCC一線方案[4-5],因此,尋找HNSCC的預(yù)后標(biāo)志物對(duì)于為患者提供有效的診療方案及延長(zhǎng)患者的生存時(shí)間存在迫切而又深遠(yuǎn)的意義。在本研究中,使用R軟件整合了GEO(gene expression omnibus)[6]和ArrayExpress數(shù)據(jù)庫[7]的芯片數(shù)據(jù),用TCGA(the cancer genome atlas)[8]的樣本數(shù)據(jù)確定了11 個(gè)與HNSCC的生存有關(guān)的候選基因,使用LASSO回歸分析及列線圖構(gòu)建了預(yù)測(cè)患者生存率的預(yù)后模型。

1 資料與方法

研究方法及技術(shù)路線如圖1。

圖1 預(yù)后模型構(gòu)建的流程圖

1.1 數(shù)據(jù)準(zhǔn)備

68 個(gè)HNSCC樣本(GSE23036)[9]的芯片數(shù)據(jù)從GEO (https://www.ncbi.nlm.nih.gov/geo/) 數(shù)據(jù)庫下載,其中包括63 個(gè)癌癥樣本及5 個(gè)正常組織樣本。E-MTAB-8588隊(duì)列的芯片數(shù)據(jù)包括117 個(gè)癌癥樣本及98 個(gè)正常組織樣本從ArrayExpress (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-8588/) 數(shù)據(jù)庫中獲取。 545 個(gè)樣本的測(cè)序數(shù)據(jù)及臨床數(shù)據(jù)從TCGA (https://portal.gdc.cancer.gov/) 數(shù)據(jù)庫獲取。

1.2 獲取GEO和ArrayExpress數(shù)據(jù)集的交集基因

以校準(zhǔn)后的P值<0.05及|log2FC|>1.0為閾值,使用R語言的limma包對(duì)GSE23036和E-MTAB-8588校正后的數(shù)據(jù)進(jìn)行差異分析。利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis,WGCNA)包[10]進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析,得到兩個(gè)隊(duì)列P<0.05的模塊。選取兩個(gè)隊(duì)列的差異基因及P值最小的模塊基因取交集,利用韋恩圖[11]進(jìn)行可視化。得到的交集基因GO和KEGG通路富集分析,GO分為3 個(gè)部分:生物過程(BP)、細(xì)胞組成(CC)和分子功能(MF)。

1.3 在TCGA隊(duì)列中建立預(yù)后模型

首先利用單因素Cox分析確定交集基因中與HNSCC生存相關(guān)的基因。利用R包glmnet進(jìn)行LASSO回歸分析建立風(fēng)險(xiǎn)模型,每個(gè)患者的風(fēng)險(xiǎn)值運(yùn)用以下公式計(jì)算:

Coefi表示模型中基因的回歸系數(shù),Expi表示基因的表達(dá)量。以中位風(fēng)險(xiǎn)值為標(biāo)準(zhǔn)將患者分為高低風(fēng)險(xiǎn)兩組。

1.4 基因預(yù)后模型的驗(yàn)證

使用Kaplan-Meier (KM)生存分析比較高低風(fēng)險(xiǎn)兩組患者的生存率,繪制操作者工作特征(ROC)曲線,利用曲線下面積(AUC)驗(yàn)證模型的預(yù)測(cè)能力。隨后,利用單因素及多因素Cox回歸分析驗(yàn)證風(fēng)險(xiǎn)評(píng)分和臨床變量(性別、年齡、腫瘤分級(jí)、分期及TN分期)的預(yù)后獨(dú)立性,利用R包rms建立列線圖預(yù)測(cè)患者生存率。在高低風(fēng)險(xiǎn)組中使用GSEA軟件(4.1.0版本)進(jìn)行基因集富集分析。在HPA(Human Protein Atlas,https://www.proteinatlas.org/)數(shù)據(jù)庫中驗(yàn)證模型中的基因表達(dá)水平。

1.5 統(tǒng)計(jì)學(xué)方法

使用PERL軟件(ActivePerl- 5.28.msi,http://www.perl.org) 及R軟件(4.1.1版本)進(jìn)行所有的數(shù)據(jù)分析。

2 結(jié) 果

2.1 在GEO及ArrayExpress中篩選出交集基因

以|log2FC|>1.0和FDR<0.05為閾值在E-MTAB-8588和GSE23036隊(duì)列中分別篩選出1 163和570 個(gè)差異基因(圖2A～B)。應(yīng)用WGCNA篩選特異模塊,在E-MTAB-8588和GSE23036隊(duì)列中選取軟閾值β=7和5,分別得到14和12 個(gè)模塊(圖2C～D);選取最有意義的模塊,E-MTAB-8588中的turquoise及GSE23036的green模塊,分別包括2 307及458 個(gè)基因,與差異基因取交集,得到60 個(gè)交集基因(圖3A)。GO富集分析顯示60 個(gè)基因在脂肪酸代謝途徑、藥物代謝過程、環(huán)氧合酶P450途徑等生物過程中富集(圖3B),KEGG顯示基因在細(xì)胞色素P450代謝、化學(xué)致癌、花生四烯酸代謝等通路富集(圖3C)。

圖2 基因差異分析及共表達(dá)網(wǎng)絡(luò)分析

圖3 交集基因及富集分析

2.2 在TCGA中建立預(yù)后模型

利用單因素回歸分析確定60 個(gè)交集基因中有11 個(gè)基因(EXPH5,ELOVL6,ACOX3,NAGK,LPIN1,TJP3,NBEAL2,SASH1,MGST2,TSPAN6,ADH7)與HNSCC的生存有關(guān)(圖4A);其中相關(guān)性分析(圖4B)顯示NBEAL2與EXPH5及SASH1顯著正相關(guān)(R=0.66)。熱圖可視化11 個(gè)基因在TCGA正常及腫瘤樣本中的表達(dá)量(圖4C)。通過LASSO回歸建立11基因的預(yù)后模型(圖5A～B),利用基因回歸系數(shù)及表達(dá)量計(jì)算風(fēng)險(xiǎn)值:riskscore=(-0.233)*EXPH5+(0.104)*ELOVL6+(-0.026)*ACOX3+(-0.036)*NAGK+(-0.197)*LPIN1+(-0.054)*TJP3+(-0.005)*NBEAL2+(0.014)*SASH1+(0.018)*MGST2+(0.022)*TSPAN6+(0.003)*ADH7基于中位風(fēng)險(xiǎn)值,498 例患者被分為高低風(fēng)險(xiǎn)兩組。圖5C顯示風(fēng)險(xiǎn)值和生存狀態(tài)的分布及預(yù)后基因隨風(fēng)險(xiǎn)值的表達(dá)量。繪制熱圖評(píng)估高低風(fēng)險(xiǎn)組中預(yù)后基因表達(dá)量與臨床變量的相關(guān)性(圖6),T分期(P<0.001),臨床分期(P<0.01)及N分期(P<0.05)在高低風(fēng)險(xiǎn)兩組中有顯著差異。

圖4 預(yù)后相關(guān)基因分析 圖5 構(gòu)建LASSO回歸風(fēng)險(xiǎn)模型

圖6 臨床變量與預(yù)后基因相關(guān)性的熱圖

2.3 驗(yàn)證模型的有效性

Kaplan-Meier生存分析顯示高風(fēng)險(xiǎn)組患者預(yù)后不良(P=8.002e-07,圖7A)。同時(shí)繪制ROC曲線,1、3、 5 年的AUC分別為0.67、 0.70、 0.66(圖7B～7D),提示模型具有一定的預(yù)測(cè)效果。進(jìn)行單因素及多因素Cox回歸分析評(píng)估風(fēng)險(xiǎn)比(HR),95%可信區(qū)間(CI)及P值,單因素回歸分析顯示風(fēng)險(xiǎn)評(píng)分、臨床分期、年齡和TN分期與患者生存有關(guān)(圖7E);多因素分析結(jié)果顯示風(fēng)險(xiǎn)評(píng)分(HR=1.784,95%CI:1.461-2.180,P<0.001),年齡(HR=1.023,95%CI:1.005-1.040,P<0.001)及N分期(HR=1.332,95%CI:1.054-1.683,P=0.016)具有獨(dú)立預(yù)后性(圖7F)。隨后,選擇上述3 個(gè)具有顯著意義的預(yù)測(cè)因子進(jìn)行列線圖的繪制(圖7G),通過繪制分值的軸向垂線計(jì)算總分值,可以預(yù)測(cè)患者的1、 3、 5 年的生存率。進(jìn)行GSEA富集分析評(píng)估預(yù)后模型相關(guān)的通路和功能。圖8A列出了前7 個(gè)最有顯著意義的通路,蛋白酶體、核糖體、嘌呤代謝、RNA聚合酶嘧啶代謝、戊糖磷酸途徑和谷胱甘肽代謝途徑富集在高風(fēng)險(xiǎn)組中,T細(xì)胞、B細(xì)胞受體信號(hào)通路、乙醚脂質(zhì)代謝、醛固酮調(diào)節(jié)鈉重吸收、原發(fā)性免疫缺陷、FCεR I信號(hào)通路以及JAK/STAT信號(hào)通路富集在低風(fēng)險(xiǎn)組。最后,使用HPA數(shù)據(jù)庫評(píng)估基于免疫組化染色的預(yù)后基因表達(dá)水平(圖8B)。結(jié)果顯示大多數(shù)預(yù)后基因在正常組織中的表達(dá)量高于腫瘤組織,與TCGA的表達(dá)水平一致。

圖7 預(yù)后模型的驗(yàn)證與列線圖的構(gòu)建

圖8 GSEA富集分析及HPA數(shù)據(jù)庫驗(yàn)證

3 討 論

目前,高通量技術(shù)已廣泛應(yīng)用于HNSCC基因組、轉(zhuǎn)錄組和蛋白組水平的研究,但如今仍然缺少較為理想且穩(wěn)定的生物標(biāo)志物。因此,通過從公共數(shù)據(jù)庫挖掘并整合高通量數(shù)據(jù)進(jìn)而建立有效的預(yù)后標(biāo)志物是有必要的。本研究通過利用R包WGCNA和limma分析GSE23036和E-MTAB-8588的芯片數(shù)據(jù)并得到60 個(gè)交集基因,在TCGA隊(duì)列中進(jìn)行單因素分析篩選出11 個(gè)基因與HNSCC生存相關(guān)。隨后利用LASSO回歸分析建立11 個(gè)基因的預(yù)后模型,單因素及多因素Cox分析證明風(fēng)險(xiǎn)評(píng)分是獨(dú)立的預(yù)后因子,最后建立了一個(gè)列線圖以輔助臨床實(shí)踐。

本研究中的基因預(yù)后模型包括11 個(gè)基因:EXPH5、ELOVL6、ACOX3、NAGK、LPIN1、TJP3、NBEAL2、SASH1、MGST2、TSPAN6和ADH7。EXPH5是一種Rab27b的效應(yīng)蛋白,可以通過調(diào)控角質(zhì)形成細(xì)胞調(diào)節(jié)表皮內(nèi)穩(wěn)態(tài)[12],在結(jié)直腸癌中發(fā)現(xiàn)其與COL1A2融合[13]。ELOVL6是超長(zhǎng)鏈脂肪酸延伸酶,主要參與飽和與不飽和脂肪酸的延伸;其高表達(dá)與肝細(xì)胞癌的良好預(yù)后有關(guān)[14]。ACOX3是一種?；o酶A氧化酶,參與支鏈脂肪酸的代謝,被認(rèn)為是宮頸癌的潛在生物標(biāo)志物[15]。LIPIN1調(diào)節(jié)脂質(zhì)和能量代謝,其過表達(dá)可促進(jìn)乳腺癌進(jìn)展[16]。NBEAL2是一種含有BEACH結(jié)構(gòu)域的蛋白質(zhì)家族成員,Wu等[17]曾研究報(bào)道NBEAL2與HNSCC的總體生存相關(guān)。TJP3是MAGUK蛋白超家族成員[18],被認(rèn)為可促進(jìn)卵巢癌細(xì)胞的侵襲和遷移[19]。SASH1是一種腫瘤抑制因子,抑制腫瘤細(xì)胞增殖、遷移及上皮間質(zhì)轉(zhuǎn)化[20];另一方面受MiR-9的負(fù)調(diào)控,促進(jìn)HNSCC的腫瘤侵襲[21]。MGST2是谷胱甘肽S-轉(zhuǎn)移酶(GST)基因家族成員[22],在UVB照射下發(fā)生甲基化并促進(jìn)皮膚癌的進(jìn)展[23]。乙醇脫氫酶(ADH)家族成員ADH7與非小細(xì)胞肺癌預(yù)后顯著相關(guān)[24]。因此,本研究預(yù)后模型中的基因與多種癌癥有關(guān),以期為HNSCC的分子基礎(chǔ)提供新的認(rèn)識(shí)。

GSEA富集分析顯示許多生物途徑影響著腫瘤的發(fā)展。抑制蛋白酶體可抑制HNSCC細(xì)胞的生長(zhǎng)[25]。循環(huán)外泌體中的嘌呤代謝水平影響著HNSCC的進(jìn)展[26]。由ASCT2和GLUD組成的谷氨酰胺代謝相關(guān)基因可作為HNSCC靶向治療的生物標(biāo)志物[27]。抑制JAK/STAT通路抑制癌細(xì)胞增殖,可作為HNSCC有效的治療手段[28]。因此,可進(jìn)一步評(píng)估本研究中GSEA富集的通路在HNSCC中的作用。

綜上所述,本研究中的風(fēng)險(xiǎn)模型可作為HNSCC的生物標(biāo)志物,建立的列線圖可指導(dǎo)臨床醫(yī)生的診斷及治療。然而,本研究存在一定的局限性,需要臨床試驗(yàn)評(píng)估風(fēng)險(xiǎn)模型的預(yù)后價(jià)值。