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

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

1 材料與方法

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

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

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

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

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

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

2 結(jié)果

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

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

紅色部分代表51個上調(diào)基因，藍(lán)色部分代表81個下調(diào)基因。圖1 12 403個差異表達(dá)的火山圖

圖2 差異表達(dá)基因的熱圖

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

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

表1 正常組織和ESCC組織的差異表達(dá)基因(FDR<0.05及l(fā)ogFC≥2或log2FC≤-2)分子功能富集分析結(jié)果

注：FDR：false discovery rate。

圖4 正常組織和ESCC組織的差異表達(dá)基因細(xì)胞組成富集分析可視化網(wǎng)絡(luò)圖

GO Path編號細(xì)胞組成PFDR基因數(shù)GO:0031012細(xì)胞外基質(zhì)1.50E-133.00E-1123GO:0005583纖維膠原三聚體6.75E-094.50E-075GO:0098643帶狀膠原原纖維6.75E-094.50E-075GO:0062023含膠原細(xì)胞外基質(zhì)1.84E-089.20E-0716GO:0044420細(xì)胞外基質(zhì)成分3.35E-081.34E-067GO:0098644膠原三聚體復(fù)合物1.20E-074.02E-065GO:0016324頂端質(zhì)膜1.26E-063.60E-0512GO:0045177細(xì)胞頂端部分1.63E-064.08E-0513GO:0001533編碼角質(zhì)化包膜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組織的差異表達(dá)基因(FDR<0.05及l(fā)ogFC≥2或log2FC≤-2)生物過程富集分析結(jié)果

注：FDR：false discovery rate。

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

編 號通路名稱PFDR基因數(shù)hsa04974蛋白質(zhì)消化吸收7.40E-068.88E-047hsa05323類風(fēng)濕性關(guān)節(jié)炎7.12E-044.27E-025hsa05146Amoebiasis1.08E-034.33E-025

注：FDR：false discovery rate。

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

A：131個差異基因蛋白相互作用網(wǎng)絡(luò)圖；B：由15個節(jié)點、105個邊組成的蛋白相互作用網(wǎng)絡(luò)圖；C：由10個節(jié)點、45個邊組成的蛋白相互作用網(wǎng)絡(luò)圖。圖7 差異表達(dá)基因蛋白—蛋白相互作用圖

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

3 討論

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

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

同時，把這些基因進(jìn)行KEGG通路分析，得到蛋白質(zhì)消化吸收通路、Amoebiasis通路、類風(fēng)濕性關(guān)節(jié)炎通路3條通路，其中COL11A1、COL1A1、COL1A2、COL5A2等COL基因家族參與了蛋白質(zhì)消化吸收通路(表4)，COL1A1是一種蛋白質(zhì)編碼基因，該基因編碼Ⅰ型膠原蛋白的親α1鏈，其三螺旋結(jié)構(gòu)由兩條α1鏈和一條α2鏈組成。該基因和血小板衍生生長因子β基因所在的第17和22號染色體之間的相互易位與一種稱為突起皮膚纖維肉瘤的特殊類型的皮膚腫瘤有關(guān)，這種腫瘤是由于生長因子的不受管制的表達(dá)而引起的，其可能是ESCC發(fā)生的潛在基因[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。紅色代表基因高表達(dá)，黑色代表基因低表達(dá)。圖8 所有差異表達(dá)基因在以TCGA數(shù)據(jù)庫中81例ESCC患者做Kaplan-Meier總生存分析

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

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

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