肺腺癌EMT相關(guān)基因篩選及其ceRNA網(wǎng)絡(luò)構(gòu)建*

馬楚涵,孫家琰,趙雨馨,趙幼敏,劉奕墨,趙子涵,王子劍,尹智華△

中國醫(yī)科大學(xué):1.第一臨床學(xué)院;2.第四臨床學(xué)院;3.公共衛(wèi)生學(xué)院;4.健康管理學(xué)院,遼寧沈陽 110122

肺癌是我國發(fā)病率與死亡率最高的惡性腫瘤[1]。腺癌是肺癌最常見的病理類型[2]。中晚期腺癌患者5年生存率較低,其原因主要是腫瘤細(xì)胞的淋巴及血道轉(zhuǎn)移導(dǎo)致術(shù)后復(fù)發(fā)。因此,對于肺腺癌侵襲、轉(zhuǎn)移的研究具有重要的臨床意義。上皮間質(zhì)轉(zhuǎn)化(EMT)是指上皮來源的細(xì)胞向間充質(zhì)細(xì)胞轉(zhuǎn)化的過程,上皮細(xì)胞通過EMT失去極性及細(xì)胞連接,獲得富有運(yùn)動能力的間充質(zhì)表型,從而突破基膜進(jìn)一步侵襲轉(zhuǎn)移。因此EMT與腫瘤轉(zhuǎn)移密切相關(guān)[3]。本研究利用基因表達(dá)圖譜(GEO)及腫瘤基因組圖譜(TCGA)數(shù)據(jù)庫,篩選肺腺癌EMT相關(guān)基因,并對其進(jìn)行了生物信息學(xué)分析及關(guān)鍵基因競爭性內(nèi)源RNA(ceRNA)網(wǎng)絡(luò)構(gòu)建,從基因組層面系統(tǒng)描述了EMT過程,為相關(guān)研究提供了新思路。

1 材料與方法

1.1標(biāo)本來源 在GEO數(shù)據(jù)庫(http://www.ncbi.nlm.nih.gov/geo/)中下載GSE10072數(shù)據(jù)集進(jìn)行分析,該數(shù)據(jù)集包含58例肺腺癌標(biāo)本及49例正常組織標(biāo)本,采用Affymetrix Human Genome U133A Array芯片檢測全基因組基因。在TCGA數(shù)據(jù)庫中下載肺腺癌項(xiàng)目數(shù)據(jù),包括516例肺腺癌標(biāo)本、56例正常組織標(biāo)本的基因表達(dá)數(shù)據(jù),使用R語言進(jìn)行后續(xù)分析。本研究所選標(biāo)本均來源于44歲以上患者,從經(jīng)組織病理學(xué)證實(shí)的肺腺癌Ⅰ～Ⅳ期患者術(shù)中切除所得,未區(qū)分性別與是否吸煙。

1.2方法

1.2.1差異表達(dá)分析 本研究使用GEO2R篩選GEO數(shù)據(jù)庫中的差異表達(dá)基因;構(gòu)建TCGA數(shù)據(jù)庫基因表達(dá)矩陣,采用R語言中的limma包進(jìn)行差異分析(P<0.05且|logFC|>1),取兩次分析結(jié)果的交集作為最終的差異表達(dá)基因。

1.2.2篩選肺腺癌EMT相關(guān)基因 將差異表達(dá)基因的Gene Symbol導(dǎo)入GenClip3(http://cismu.net/genclip3/analysis.php)。軟件將從Pubmed中檢索相關(guān)文獻(xiàn),并基于關(guān)鍵詞對于導(dǎo)入的基因進(jìn)行聚類分析[4]。

1.2.3基因本體論(GO)及京都基因和基因組百科全書(KEGG)富集分析 采用Metascape(http://metascape.org/)對EMT相關(guān)基因進(jìn)行功能注釋及GO與KEGG富集分析[5]。

1.2.4構(gòu)建蛋白互作網(wǎng)絡(luò)(PPI)及篩選關(guān)鍵基因 采用STRING在線工具(https://cn.string-db.org/)分析EMT相關(guān)基因編碼的蛋白質(zhì)之間的相互作用[6],得到PPI,將其中結(jié)合分?jǐn)?shù)值>0.4的蛋白相互作用對導(dǎo)入Cytoscape軟件中,去除離散蛋白,通過cytoNCA分析蛋白之間的相關(guān)程度。根據(jù)degree排序,選擇top 10的基因作為關(guān)鍵基因,并針對其繪制PPI網(wǎng)絡(luò)圖。

1.2.5預(yù)測EMT關(guān)鍵基因潛在治療中藥 將關(guān)鍵基因與肺腺癌、EMT共同導(dǎo)入Coremine medical數(shù)據(jù)庫(http://www.coremine.com/),以P<0.05為標(biāo)準(zhǔn)篩選治療性中藥。將針對10個(gè)關(guān)鍵基因的中藥取交集,并繪制中藥與關(guān)鍵基因相關(guān)性的氣泡圖。

1.2.6Kaplan-Meier分析 通過Kaplan-Meier(www.kmplot.com)對GEO、歐洲基因組表型檔案(EGA)及TCGA中的患者數(shù)據(jù)進(jìn)行分析,以50個(gè)月的生存率評估患者預(yù)后[7]。最后,選取Log-rankP<0.05的基因作進(jìn)一步研究。

1.2.7根據(jù)ceRNA假說預(yù)測關(guān)鍵基因上游微小RNA(miRNA)及miRNA上游的長鏈非編碼RNA(lncRNA) 根據(jù)ceRNA假說,lncRNA能夠競爭性抑制miRNA對靶基因的調(diào)控作用,因此,靶基因的表達(dá)水平與miRNA呈負(fù)相關(guān)而與lncRNA呈正相關(guān)。使用miRTarbase(https:// mirtarbase.cuhk.edu.cn/)預(yù)測關(guān)鍵基因上游miRNA,為提高研究結(jié)果可信度,只選擇研究報(bào)道經(jīng)蛋白質(zhì)印跡法、實(shí)時(shí)熒光定量PCR證實(shí)的miRNA納入研究[8]。使用ENCORI(http://starbase.sysu.edu.cn/)計(jì)算logFC,評估m(xù)iRNA及關(guān)鍵基因在正常組織與腫瘤組織中的差異表達(dá)情況,并繪制生存曲線,以Log-rankP<0.05判定基因表達(dá)與患者預(yù)后的關(guān)系。選擇與患者預(yù)后相關(guān)并與靶基因在腫瘤組織與正常組織間差異表達(dá)趨勢相反的miRNA進(jìn)行下一步研究。使用miRNet database(https://www.mirnet.ca/)預(yù)測miRNA上游lncRNA[9],通過ENCORI選取與患者預(yù)后相關(guān)并與靶基因在腫瘤組織與正常組織間差異表達(dá)趨勢相同的lncRNA作進(jìn)一步研究。

1.2.8構(gòu)建ceRNA網(wǎng)絡(luò) 根據(jù)上一步預(yù)測的miRNA及l(fā)ncRNA,使用ENCORI繪制散點(diǎn)圖,通過Pearson相關(guān)性檢驗(yàn)判斷基因表達(dá)的相關(guān)性,并繪制關(guān)鍵基因的ceRNA調(diào)控網(wǎng)絡(luò)圖[10]。

1.3統(tǒng)計(jì)學(xué)處理 本研究采用R軟件(4.1.0版)和Microsoft Excel 2019對TCGA數(shù)據(jù)進(jìn)行整理。應(yīng)用R軟件(4.1.2版)對RNA表達(dá)數(shù)據(jù)進(jìn)行差異分析(P<0.05且|logFC|>1),并將其與GEO2R所得差異表達(dá)基因取交集獲得最終結(jié)果。然后,使用GenClip3軟件進(jìn)行聚類分析,并采用Metascape對目標(biāo)基因進(jìn)行GO和KEGG富集分析。同時(shí),采用STRING在線工具預(yù)測蛋白質(zhì)之間的相互作用,通過Cytoscape軟件選擇連接度最高的10個(gè)基因作為關(guān)鍵基因。在此基礎(chǔ)上使用Coremine medical數(shù)據(jù)庫篩選治療性中藥(P<0.05),并且結(jié)合生存數(shù)據(jù)進(jìn)行Kaplan-Meier(Log-rankP<0.05)分析。最后,使用在線預(yù)測網(wǎng)站構(gòu)建ceRNA,通過Pearson相關(guān)性檢驗(yàn)判斷基因表達(dá)的相關(guān)性。

2 結(jié) 果

2.1差異表達(dá)分析 使用GEO2R在線分析工具分析,得到差異表達(dá)基因662個(gè);在TCGA數(shù)據(jù)庫中使用limma包進(jìn)行差異表達(dá)分析,得到差異表達(dá)基因8 295個(gè)。取交集得到差異表達(dá)基因565個(gè),見圖1。

圖1 GEO和TCGA數(shù)據(jù)庫差異表達(dá)基因的韋恩圖

2.2EMT相關(guān)基因篩選 將565個(gè)差異表達(dá)基因?qū)隚enClip3中,結(jié)果顯示所有差異表達(dá)基因均可在文獻(xiàn)中檢索到,使用GenClip3對差異表達(dá)基因進(jìn)行聚類分析,得到156個(gè)EMT相關(guān)基因。

2.3GO、KEGG富集分析 GO分析發(fā)現(xiàn),這些基因主要與血管發(fā)育、激素應(yīng)答、生長因子刺激的細(xì)胞應(yīng)答、細(xì)胞遷移的正向調(diào)節(jié)、細(xì)胞外基質(zhì)等生物過程、細(xì)胞組件、分子功能相關(guān),見圖2。KEGG通路分析發(fā)現(xiàn),這些基因主要與癌通路、糖尿病并發(fā)癥中的高級糖基化終末產(chǎn)物-受體(AGE-RAGE)通路、缺氧誘導(dǎo)因子-1(HIF-1)信號通路、腫瘤中的miRNA通路、細(xì)胞黏附分子通路等相關(guān),見圖3。

圖2 EMT相關(guān)基因的GO富集圖

圖3 EMT相關(guān)基因的KEGG富集圖

2.4PPI構(gòu)建 根據(jù)degree排序,得到的10個(gè)關(guān)鍵基因,分別為鈣黏蛋白1(CDH1)、白細(xì)胞介素-6(IL-6)、基質(zhì)金屬蛋白酶9(MMP9)、血小板內(nèi)皮細(xì)胞黏附分子(PECAM1)、細(xì)胞周期蛋白依賴性激酶抑制劑2A(CDKN2A)、α1-Ⅰ型膠原基因(COL1A1)、分泌型磷酸蛋白1(SPP1)、TIMP基質(zhì)金屬蛋白酶抑制因子1(TIMP1)、小窩蛋白1(CAV1)、zeste 基因增強(qiáng)子同源物(EZH2),見表1。從PPI網(wǎng)絡(luò)圖中可以看出,關(guān)鍵基因編碼的蛋白之間相互作用明顯,見圖4。

表1 關(guān)鍵基因及degree

圖4 關(guān)鍵基因的PPI

2.5關(guān)鍵基因潛在治療性中藥 10個(gè)關(guān)鍵基因預(yù)測到紫草、溫莪術(shù)、片姜黃等多味中藥,見表2。將10個(gè)關(guān)鍵基因的潛在治療性中藥取交集,得到丹參、人參蘆、人參、人參葉、人參花共五味中藥,以P值為參考,使用R軟件的ggpubr包繪制關(guān)鍵基因與中藥相關(guān)性的氣泡圖,見圖5。

表2 針對關(guān)鍵基因的部分治療性中藥

圖5 中藥與關(guān)鍵基因的相關(guān)性氣泡圖

2.6關(guān)鍵基因潛在治療性中藥的Kaplan-Meier分析 結(jié)果顯示,10個(gè)關(guān)鍵基因均與患者預(yù)后相關(guān)(Log-rankP<0.05),其中高表達(dá)CDH1、PECAM1、CAV1的生存曲線與對照的風(fēng)險(xiǎn)比(HR)<1,表明這些基因的表達(dá)可以降低患者死亡風(fēng)險(xiǎn);高表達(dá)IL-6、MMP9、CDKN2A、COL1A1、SPP1、TIMP1的生存曲線與對照的HR>1,表明這些基因的表達(dá)增加患者的死亡風(fēng)險(xiǎn)。見圖6。

圖6 關(guān)鍵基因的Kaplan-Meier分析結(jié)果森林圖

2.7根據(jù)ceRNA假說預(yù)測關(guān)鍵基因上游miRNA及miRNA上游lncRNA 差異分析結(jié)果顯示,IL-6(logFC=-2.64)在腫瘤組織中呈低表達(dá),而MMP9(logFC=2.14)、COL1A1(logFC=3.05)、SPP1(logFC=4.95)、EZH2(logFC=2.70)在腫瘤組織中呈高表達(dá),同理可知IL-6的上游miRNA hsa-miR-9-5p,MMP9的上游miRNA hsa-miR-211-5p,COL1A1的上游miRNA hsa-miR-133a-3p,SSP1的上游miRNA hsa-miR-299-5p,EZH2的上游miRNA hsa-miR-101-3p、hsa-miR-30d-5p與其靶基因的差異表達(dá)趨勢相反,符合ceRNA假說;并且上述miRNA的Log-rankP<0.05,與患者預(yù)后相關(guān),因此將其納入下一步研究。由logFC可以看出,hsa-miR-9-5p的上游lncRNA CCDC13-AS1、LINC00996、SH3BP5-AS1、LINC00921,hsa-miR-211-5p的上游lncRNA EBLN3P、DHRS4-AS1、FAM30A、ZNF674-AS1,hsa-miR-133a-3p的上游lncRNA LINC00847、LINC02535、MIR4697HG,hsa-miR-299-5p的上游lncRNA GUSBP11,hsa-miR-101-3p的上游lncRNA LINC01579、EBLN3P、GSEC,hsa-miR-30d-5p的上游lncRNA HCG18、LINC02535與靶基因的差異表達(dá)趨勢相同,而與上游miRNA差異表達(dá)趨勢相反,符合ceRNA假說;并且上述lncRNA的Log-rankP<0.05,與患者預(yù)后相關(guān),因此納入下一步研究。見圖7～9。

圖7 mRNA-miRNA-lncRNA相互作用的桑基圖

圖8 基因差異表達(dá)分析的偏差圖

圖9 基因HR的棒棒糖圖

2.8構(gòu)建ceRNA網(wǎng)絡(luò) 在預(yù)測的miRNA及l(fā)ncRNA中,同時(shí)滿足mRNA的表達(dá)水平與miRNA呈負(fù)相關(guān)而與lncRNA呈正相關(guān),lncRNA的表達(dá)水平與miRNA呈負(fù)相關(guān)的基因被納入ceRNA網(wǎng)絡(luò)。結(jié)果表明,EZH2/hsa-miR-101-3p/GSEC符合ceRNA假說,見圖10～16。

圖10 EZH2在腫瘤組織和正常組織中表達(dá)水平差異的箱線圖

圖11 EZH2的Kaplan-Meier分析曲線

圖12 hsa-miR-101-3p在腫瘤組織和正常組織中表達(dá)水平差異的箱線圖

圖13 hsa-miR-101-3p的生存曲線

圖14 GSEC在腫瘤組織和正常組織中表達(dá)水平差異的箱線圖

圖15 GSEC的生存曲線

圖16 ceRNA網(wǎng)絡(luò)及各組分之間表達(dá)水平的相關(guān)性

3 討 論

本研究獲得了565個(gè)肺腺癌差異表達(dá)基因及156個(gè)EMT相關(guān)基因,并對其進(jìn)行了功能富集分析。隨后繪制了EMT相關(guān)基因的PPI網(wǎng)絡(luò)圖,并根據(jù)degree獲得10個(gè)EMT關(guān)鍵基因。其中,CDH1編碼E-鈣黏蛋白(E-cadherin),E-cadherin下調(diào)被廣泛報(bào)道為EMT的標(biāo)志,其表達(dá)下降與腫瘤患者的不良預(yù)后密切相關(guān),在正常細(xì)胞中,E-cadherin可以阻止β-catenin,與參與wnt信號通路的淋巴樣增強(qiáng)因子/T細(xì)胞因子結(jié)合發(fā)揮腫瘤抑制作用[11]。IL-6可以通過Janus激酶/信號轉(zhuǎn)導(dǎo)子和轉(zhuǎn)錄激活子信號通路誘導(dǎo)EMT[12],但也有研究顯示IL-6可以動員免疫反應(yīng)抑制腫瘤生長。MMP9通過參與細(xì)胞外基質(zhì)分解促進(jìn)腫瘤侵襲轉(zhuǎn)移,有研究顯示STAT3過表達(dá)時(shí)促進(jìn)了EMT過程,同時(shí)MMP9的表達(dá)水平也有所升高,提示MMP9在EMT中發(fā)揮一定的作用[13]。PECAM1又名CD31,其編碼的蛋白質(zhì)存在于單核細(xì)胞、血小板、中性粒細(xì)胞和某些類型的T細(xì)胞表面,并構(gòu)成細(xì)胞連接,CD31可以促進(jìn)EMT從而促進(jìn)結(jié)腸癌腹膜轉(zhuǎn)移[14],也有研究將CD31作為血管生成的標(biāo)志物[15],說明CD31也可能參與腫瘤血道轉(zhuǎn)移。CDKN2A是一種抑癌基因,但在結(jié)腸癌中CDKN2A的表達(dá)增加促進(jìn)EMT及腫瘤轉(zhuǎn)移[16],并且可以與趨化素樣因子超家族8、整合素連接激酶等通過轉(zhuǎn)化生長因子β通路、Notch通路等形成免疫抑制微環(huán)境,從而影響患者的預(yù)后[17]。COL1A1是Ⅰ型膠原蛋白的主要成分,通過結(jié)合integrin的β1亞基激活局部黏著斑激酶(FAK)。磷酸化的FAK促進(jìn)E-cadherin/catenin復(fù)合物的分解,從而激活EMT[18]。SPP1上調(diào)在前列腺癌、結(jié)直腸癌中都有激活EMT,促進(jìn)腫瘤轉(zhuǎn)移的作用,可能與SPP1激活細(xì)胞外信號相關(guān)激酶1/2信號通路有關(guān)[19-20]。TIMP1為MMP1的抑制劑,但其激活EMT的作用是通過與CD63結(jié)合誘導(dǎo)TWIST1的表達(dá),與MMP抑制功能無關(guān)。CAV1可以通過Wnt/β-catenin通路誘導(dǎo)EMT過程,促進(jìn)腫瘤轉(zhuǎn)移。EZH2促進(jìn)絕大多數(shù)腫瘤的侵襲轉(zhuǎn)移,但有時(shí)也發(fā)揮抗腫瘤作用,這種矛盾現(xiàn)象與EMT密切相關(guān)[21]。對10個(gè)關(guān)鍵基因的Kaplan-Meier分析顯示,IL-6、MMP9、CDKN2A、COL1A1、SPP1、TIMP1的表達(dá)與患者預(yù)后不良有關(guān),而上述研究顯示這些基因?qū)τ贓MT均起促進(jìn)作用,可能是其影響患者生存期的原因之一。并使用Coremine medical數(shù)據(jù)庫預(yù)測了10個(gè)關(guān)鍵基因的治療性中藥,人參、人參葉、人參花、人參蘆、丹參這五味中藥對于所有關(guān)鍵基因均有一定效果,本研究還根據(jù)P值繪制了這五味中藥與關(guān)鍵基因的相關(guān)性氣泡圖。從氣泡圖中分析,人參、人參葉、人參花、人參蘆與EMT的關(guān)系較丹參更為密切,而人參、人參葉、人參花、人參蘆來源于同一種植物,因此人參中某種成分可能有效抑制EMT從而延緩腫瘤進(jìn)展。有研究表明,人參皂苷可以抑制EMT從而抑制腫瘤生長與轉(zhuǎn)移[22],可能為其潛在機(jī)制之一。一項(xiàng)針對隨機(jī)對照臨床試驗(yàn)的Meta分析表明,人參的主要成分人參皂苷和多糖作為化療輔助用藥能夠提高非小細(xì)胞肺癌患者的客觀緩解率和疾病控制率,降低患者不良反應(yīng)發(fā)生率[23]。根據(jù)ceRNA假說,本研究構(gòu)建了EZH2/hsa-miR-101-3p/GSEC調(diào)控網(wǎng)絡(luò)。研究表明hsa-miR-101-3p可以通過調(diào)控EZH2的表達(dá)抑制腫瘤進(jìn)展,在肺鱗癌中hsa-miR-101-3p可以通過下調(diào)EZH2抑制腫瘤活性,促進(jìn)腫瘤細(xì)胞凋亡[24];在腎癌中,hsa-miR-101-3p可以通過調(diào)控EZH2抑制腫瘤的侵襲轉(zhuǎn)移[25];hsa-miR-101-3p還可以與Syn-cal14.1a協(xié)同作用抑制EZH誘導(dǎo)的乳腺癌進(jìn)展[26]。GSEC/hsa-miR-101-3p與肺腺癌細(xì)胞的鐵死亡相關(guān),上調(diào)GSEC或下調(diào)hsa-miR-101-3p的表達(dá)可以抑制肺腺癌細(xì)胞的增殖與遷移能力[27]。然而EZH2與GSEC之間的相互作用及EZH2/hsa-miR-101-3p/GSEC在調(diào)控EMT中的作用仍需進(jìn)一步研究。

綜上所述,本研究利用GEO、TCGA數(shù)據(jù)庫及GenClip3文獻(xiàn)挖掘平臺等分析工具,篩選肺腺癌EMT過程的相關(guān)基因并對其進(jìn)行了功能富集分析,以及PPI的構(gòu)建,從基因組層面對EMT過程進(jìn)行了系統(tǒng)描述,加深了對于EMT過程的認(rèn)識。同時(shí),篩選出了EMT過程中的關(guān)鍵基因并預(yù)測了人參等針對關(guān)鍵基因的治療性中藥,從患者預(yù)后角度構(gòu)建了EZH2/hsa-miR-101-3p/GSEC調(diào)控網(wǎng)絡(luò),為肺腺癌侵襲、轉(zhuǎn)移的治療提供了新思路,為肺腺癌預(yù)后標(biāo)志物的研究提供了參考。