基于WGCNA分析及風險評分構建頭頸部鱗狀細胞癌的預后模型

金祥

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

1 資料與方法

研究方法及技術路線如圖1。

圖1 預后模型構建的流程圖

1.1 數據準備

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

1.2 獲取GEO和ArrayExpress數據集的交集基因

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

1.3 在TCGA隊列中建立預后模型

首先利用單因素Cox分析確定交集基因中與HNSCC生存相關的基因。利用R包glmnet進行LASSO回歸分析建立風險模型,每個患者的風險值運用以下公式計算:

Coefi表示模型中基因的回歸系數,Expi表示基因的表達量。以中位風險值為標準將患者分為高低風險兩組。

1.4 基因預后模型的驗證

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

1.5 統計學方法

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

2 結 果

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

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

圖2 基因差異分析及共表達網絡分析

圖3 交集基因及富集分析

2.2 在TCGA中建立預后模型

利用單因素回歸分析確定60 個交集基因中有11 個基因(EXPH5,ELOVL6,ACOX3,NAGK,LPIN1,TJP3,NBEAL2,SASH1,MGST2,TSPAN6,ADH7)與HNSCC的生存有關(圖4A);其中相關性分析(圖4B)顯示NBEAL2與EXPH5及SASH1顯著正相關(R=0.66)。熱圖可視化11 個基因在TCGA正常及腫瘤樣本中的表達量(圖4C)。通過LASSO回歸建立11基因的預后模型(圖5A～B),利用基因回歸系數及表達量計算風險值: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基于中位風險值,498 例患者被分為高低風險兩組。圖5C顯示風險值和生存狀態的分布及預后基因隨風險值的表達量。繪制熱圖評估高低風險組中預后基因表達量與臨床變量的相關性(圖6),T分期(P<0.001),臨床分期(P<0.01)及N分期(P<0.05)在高低風險兩組中有顯著差異。

圖4 預后相關基因分析 圖5 構建LASSO回歸風險模型

圖6 臨床變量與預后基因相關性的熱圖

2.3 驗證模型的有效性

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

圖7 預后模型的驗證與列線圖的構建

圖8 GSEA富集分析及HPA數據庫驗證

3 討 論

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

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

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

綜上所述,本研究中的風險模型可作為HNSCC的生物標志物,建立的列線圖可指導臨床醫生的診斷及治療。然而,本研究存在一定的局限性,需要臨床試驗評估風險模型的預后價值。