R語言下食管鱗狀細胞癌關鍵驅動基因富集分析的生物信息學研究

王子明，李新陽，姚 歌，王鈺鯤，王新帥，原 翔，張廣平

食管癌是常見的消化系統惡性腫瘤之一，其發病率和病死率在全球惡性腫瘤中分別居第 8 位和第6位，每年死亡人數可達400 000人[1-2]。食管癌在世界范圍內呈現出顯著的地域性和一定的民族差異，在歐美國家以腺癌居多，在亞洲國家以鱗狀細胞癌為主，食管鱗狀細胞癌(esophageal squamous cell carcinoma,ESCC)是食管癌的主要組織學亞型，約占90%[3]。盡管對食管癌的多種治療方法已取得一定的進展，但食管癌仍然是一種惡性疾病，在西方國家總體5 a生存率為12%～20%[4]。近年來，不斷有研究嘗試食管癌生物靶向治療手段，探索相關分子作用的靶點及相關信號，如針對表皮生長因子受體(epidermalgrowth factor receptor,EGFR)、針對血管內皮細胞生長因子(vascularendothelial growth factor,VEG)及VEGF受體(vascularendothelial growth factor receptor,VEGFR)、針對原癌基因人類表皮生長因子受體2(human epidermal growth factor receptor 2,HER2)等。但是，目前根據已發現的靶點而研究出的靶向藥物較少。因此，本研究通過運用R語言及其相關軟件包，在基因組水平上探索食管癌患者發生、發展的潛在分子靶點及相關信號途徑，為進一步發現新的食管癌治療靶標提供重要的理論依據。

1 材料與方法

1.1 GEO數據庫及數據來源和處理GEO數據(http://www.ncbi.nlm.nih.gov/geo/)庫全稱GENE EXPRESSION OMNIBUS，是由美國國立生物技術信息中心NCBI創建并維護的基因表達數據庫。它創建于2000年，收錄了世界各國研究機構提交的高通量基因表達數據。本研究運用R語言平臺(3.6.0版本)下的Rstudio(3.6.0版本)和相關的軟件包，從GEO數據庫中獲得數據GSE38129。GSE38129平臺注釋信息：GPL571[HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array。對微表達矩陣芯片數據進行去除缺失值、log2處理、標準化、取平均表達量最高的探針用于探針ID和基因名匹配，將篩選得到的基因表達矩陣數據用作本研究的數據來進行分析。

1.2 差異表達基因矩陣選定運用R語言下limma包(3.40.2版本)對腫瘤和正常樣本組織篩選得到的微表達矩陣芯片數據進行進一步篩選，以“FDR<0.05及logFC≥2或log2FC≤-2”為篩選閾值，得到顯著差異表達基因矩陣，并用R語言下ggplot2包(3.1.1版本)繪制火山圖和熱圖來驗證已得到的差異表達基因矩陣。

1.3 富集分析分析差異表達基因本研究對差異表達基因分別進行包括分子功能(molecular function)、細胞組成(cellular component)及生物過程(biological process)的GO功能富集分析和KEGG通路富集分析，并將所得數據分別通過ggplot2包和pathview包(1.24.0版本)進行GO功能富集分析可視化展示和KEGG通路富集分析可視化圖示。

1.4 顯著表達差異基因蛋白—蛋白相互作用(PPI)網絡構建和分析STRING數據庫(https://string-db.org/)是一個用于構建蛋白之間的相互作用網絡并分析和探索預測的蛋白相互作用網絡在線數據庫，本研究將篩選得到的差異基因上傳到STRING數據庫(10.5版本)和Cytoscape(3.7.2版本)及MCODE插件用于差異基因PPI網絡可視化分析，并以“combined score>0.4、MCODE>2為篩選閾值，得到相應的功能模塊網絡圖，為了進一步分析關鍵基因，選取其中MCODE Score>5，得到2個功能模塊網絡圖，從而得到潛在的關鍵基因。

1.5 生存預后分析在本研究中，所得關鍵基因均通過運用以The Cancer Genome Atlas(TCGA)數據庫中81例ESCC患者數據為參考數據的Kaplan-Meier plotter(http://kmplot.com/analysis/)ESCC數據庫做生存預后分析。

1.6 統計方法P<0.05被認為具有顯著的統計學差異，多檢驗假設采用false discovery rate (FDR)方法對P值進行調整，且FDR<0.05作為篩選差異表達基因的閾值，在進行富集分析時運用超幾何分布方法，同時本研究使用ESCC Kaplan-Meier曲線并用log-rank檢查進行了比較。

2 結果

2.1 ESCC組織樣本和正常細胞組織樣本數據清洗本研究運用R語言相關函數從GEO數據庫獲得GSE38129數據(30個正常樣本組織，30個ESCC樣本組織)，共得到22 277個微表達矩陣芯片數據，通過對數據進行初步清洗，得到12 403個基因表達數據矩陣。

2.2 表達基因矩陣選定對所得12 403個基因表達矩陣通過運用R語言下limma包進行差異表達分析，得到12 403個差異基因表達矩陣(上調基因或者下調基因)。以“FDR<0.05及logFC≥2或log2FC≤-2”為閾值，共得到132個顯著差異表達基因，其中51個上調基因、81個下調基因(圖1)，同時將差異表達基因進行可視化做圖(圖2)。

紅色部分代表51個上調基因，藍色部分代表81個下調基因。圖1 12 403個差異表達的火山圖

圖2 差異表達基因的熱圖

2.3 富集分析差異表達基因通過clusterProfiler包對篩選得到的差異表達基因進行富集分析。首先分別對差異表達基因進行分子功能、細胞組成、生物過程富集分析，發現這些基因表達產物分別表現肽鏈內切酶絲氨酸型肽酶活性、細胞外基質結構成分、絲氨酸型肽酶活性、絲氨酸水解酶活性等分子功能上(圖3、表1)，細胞外基質、纖維膠原三聚體、編碼角質化包膜、紡錘體等細胞組成中(圖4、表2)，顯著富集在表皮生長、細胞外基質組織、膠原蛋白分解代謝過程、女性減數分裂核分裂等生物過程內(圖5、表3)。其次，對差異基因進行KEGG通路富集分析，發現這些基因顯著富集在蛋白質消化吸收通路、Amoebiasis通路、類風濕性關節炎通路上(表4)。

圖3 正常組織和ESCC組織的差異表達基因分子功能富集分析可視化熱圖

表1 正常組織和ESCC組織的差異表達基因(FDR<0.05及logFC≥2或log2FC≤-2)分子功能富集分析結果

注：FDR：false discovery rate。

圖4 正常組織和ESCC組織的差異表達基因細胞組成富集分析可視化網絡圖

GO Path編號細胞組成PFDR基因數GO:0031012細胞外基質1.50E-133.00E-1123GO:0005583纖維膠原三聚體6.75E-094.50E-075GO:0098643帶狀膠原原纖維6.75E-094.50E-075GO:0062023含膠原細胞外基質1.84E-089.20E-0716GO:0044420細胞外基質成分3.35E-081.34E-067GO:0098644膠原三聚體復合物1.20E-074.02E-065GO:0016324頂端質膜1.26E-063.60E-0512GO:0045177細胞頂端部分1.63E-064.08E-0513GO:0001533編碼角質化包膜1.32E-052.93E-045GO:0030496中間體2.64E-055.27E-048GO:0005581膠原蛋白三聚體3.02E-055.49E-046GO:0005819紡錘體4.67E-047.78E-039

注：FDR：false discovery rate。

表3 正常組織和ESCC組織的差異表達基因(FDR<0.05及logFC≥2或log2FC≤-2)生物過程富集分析結果

注：FDR：false discovery rate。

圖5 正常組織和ESCC組織的差異表達基因生物過程富集分析可視化柱狀圖

編 號通路名稱PFDR基因數hsa04974蛋白質消化吸收7.40E-068.88E-047hsa05323類風濕性關節炎7.12E-044.27E-025hsa05146Amoebiasis1.08E-034.33E-025

注：FDR：false discovery rate。

2.4 差異表達基因蛋白—蛋白相互作用網絡構建和分析本研究將得到的132個差異表達基因上傳到STRING數據庫進行差異基因蛋白—蛋白相互作用網絡構建(圖7A)。為了進一步識別關鍵基因，本研究將STRING數據庫得到的網絡圖數據上傳至Cytoscape軟件。將combined score>0.4、MCODE>2設定為篩選閾值，得到4個功能模塊網絡圖，選取其中MCODE Score>5模塊，得到2個功能模塊網絡圖(圖7B、C)，根據連接度和MCODE評分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9個關鍵基因。

A：131個差異基因蛋白相互作用網絡圖；B：由15個節點、105個邊組成的蛋白相互作用網絡圖；C：由10個節點、45個邊組成的蛋白相互作用網絡圖。圖7 差異表達基因蛋白—蛋白相互作用圖

2.5 生存預后分析所篩選出來的差異表達基因均以TCGA數據庫中ESCC生存分析數據為參考數據，通過Kaplan-Meier plotter進行生存預后分析，根據連接度和MCODE評分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9個關鍵基因與生存預后關系(圖8)。

3 討論

ESCC惡性程度較高，預后較差，目前ESCC發生、發展的確切分子作用機制及信號通路尚不清楚[4]。本研究利用ESCC和正常組織微表達矩陣芯片的表達數據，進行ESCC中的特異性表達基因的生物學分析。

在本研究中，通過以“FDR<0.05及logFC≥2或log2FC≤-2”為選定閾值，篩選出來51個上調基因和81個下調基因，并把這些基因進行生物信息學分析。首先，對這些基因進行富集分析、分子功能富集分析表明在細胞外基質結構成分、絲氨酸型肽酶活性等處顯著富集(表1)，VCAN、MMP1、CXCL12、ATAD2等差異基因參與了這些分子功能。ATAD2基因可能是核受體ESR 1的轉錄輔助激活因子，可誘導雌激素靶基因的表達，如CCND1、MYC和E2F1，可參與雌激素誘導的乳腺癌細胞增殖和細胞周期進程，可能在某些ESR1靶基因啟動子中參與CREBBP的合成或占據[5]，可能導致ESCC的發生、發展。細胞組成富集分析顯示在細胞外基質、紡錘體、染色體分離等處顯著富集(表2)，VCAN、TPO2A、MMP1、SSP1等差異基因參與了這些細胞組成，SPP1基因編碼的蛋白是一種細胞因子，可上調干擾素-γ和白細胞介素-12的表達。該基因對細胞-基質相互作用很重要，且作為細胞因子參與促進干擾素-γ和白細胞介素-12的產生和減少、白細胞介素-10的產生，在導致Ⅰ型免疫的途徑中扮演重要角色[6-7]。生物過程富集分析表明在細胞外基質組織、膠原蛋白分解代謝過程等處顯著富集(表3)，DLGAP5、VCAN、TPO2A、MMP1、MMP13、SSP1等差異基因參與了這些生物過程，DLGAP 5是一個蛋白質編碼基因。該基因在泛素-蛋白酶體通路調控有絲分裂磷蛋白，是細胞黏附、連接、完整性和分化的關鍵調節因子，可能是在細胞癌變過程中起作用的潛在細胞周期調節因子[8-9]，從而可能導致ESCC細胞過度增殖分化。MMP1基因是11號染色體上MMP基因簇的一部分，其編碼基質金屬蛋白酶(MMPs)M10家族的一個成員。在正常的生理過程中，如胚胎發育、繁殖和組織重塑，以及關節炎和轉移等疾病過程中，該家族的蛋白質參與細胞外基質的破壞。先前研究表明MMP1可能與膀胱癌、口腔鱗狀細胞癌等惡性腫瘤相關[10]，也可能參與ESCC發生過程。

同時，把這些基因進行KEGG通路分析，得到蛋白質消化吸收通路、Amoebiasis通路、類風濕性關節炎通路3條通路，其中COL11A1、COL1A1、COL1A2、COL5A2等COL基因家族參與了蛋白質消化吸收通路(表4)，COL1A1是一種蛋白質編碼基因，該基因編碼Ⅰ型膠原蛋白的親α1鏈，其三螺旋結構由兩條α1鏈和一條α2鏈組成。該基因和血小板衍生生長因子β基因所在的第17和22號染色體之間的相互易位與一種稱為突起皮膚纖維肉瘤的特殊類型的皮膚腫瘤有關，這種腫瘤是由于生長因子的不受管制的表達而引起的，其可能是ESCC發生的潛在基因[11-12]。

A：DLGAP5，P=0.000 19；B：MMP1，P=0.12；C：TOP2A，P=0.018；D：MMP13，P=0.12；E：CXCL12，P=0.22；F：SPP1，P=0.089；G：ATAD2，P=0.001 9；H：PBK，P=0.000 62；I：VCAN，P=0.003 7。紅色代表基因高表達，黑色代表基因低表達。圖8 所有差異表達基因在以TCGA數據庫中81例ESCC患者做Kaplan-Meier總生存分析

其次，把這些差異表達基因上傳至STRING數據庫和 Cytoscape軟件(圖7A、B、C)，根據連接度和MOCDE評分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9個關鍵基因。最后，9個關鍵基因在Kaplan-Meier plotter數據庫做生存分析(圖8)。PBK、VCAN基因與食管癌患者生存期明顯相關(P<0.01)(圖8H、I)，到目前為止，沒有文章或者報道出PBK、VCAN和ESCC發生、發展相關，PBK基因編碼絲氨酸/蘇氨酸蛋白激酶，與雙特異性絲裂原活化蛋白激酶(MAPKK)家族相關。這一基因的過度表達與腫瘤的發生有關，有關研究已證實PBK參與了血液學腫瘤發生[13-14]，可能參與了ESCC發生。VCAN基因編碼的蛋白是一種大型硫酸軟骨素蛋白多糖，是細胞外基質的主要組成部分。該基因編碼的蛋白參與細胞黏附、增殖、遷移和血管生成，在組織形態發生和維持中起核心作用。此基因可能在細胞間信號傳遞和連接細胞與細胞外基質中發揮作用，可能參與細胞運動、生長和分化的調控[15-16]，可能對ESCC的發生起到了調控作用。

除此之外，DLGAP5、TPO2A、ATAD2也與ESCC患者生存期顯著相關(P<0.01)(圖8A、C、G)，既往研究表明，DLGAP5與前列腺癌、結直腸癌生存預后相關[17-18]，可能參與了惡性膠質母細胞瘤、肝細胞癌等惡性腫瘤發生、發展過程[19-20]。ATAD2可以調控眾多經典的腫瘤相關信號通路，其中包括誘導MYC、E2F和甲狀腺激素受體和視黃醛激活因子，ATAD2表達與腫瘤的組織學分級以及腫瘤轉移、疾病復發和患者生存期等相關，可作為腫瘤預后評價分子，具有重要的臨床應用價值，且在前列腺癌、乳腺癌、肺腺癌等多種惡性腫瘤的發生和發展過程中發揮重要作用[21]。TOP2A基因編碼一種DNA拓撲異構酶，這種核酶參與了染色體的濃縮、染色單體的分離以及DNA轉錄和復制過程中扭轉應激的緩解。它催化了兩條雙鏈DNA的瞬間斷裂和重新連接，使兩股DNA相互通過，從而改變了DNA的拓撲結構[22-23]。TOP2A參與了結腸癌、惡性膠質母細胞瘤、膽管癌、乳腺癌等多種惡性腫瘤發生過程[24-25]。MMP1、MMP13、CXCL12、SPP1基因高表達與ESCC患者生存期呈負相關(P>0.05)(圖8B、D、E、F)，無統計學意義，可能與本研究所選ESCC患者樣本數量較少有關。

總之，本研究通過ESCC與正常組織的基因表達譜綜合生物信息學分析，表明PBK、VCAN、DLGAP5、ADAT2、TOP2A可能是ESCC發生、發展和預后的關鍵基因。雖然數據挖掘和整合是預測惡性腫瘤生物標志物的一種很有前途的工具，但是本研究只是一個初步的結論，且本研究中的樣本數量是有限的。由于腫瘤生物標志物只有與臨床數據相結合才有意義，因此需要進行進一步的實驗來證實本研究的結論。