肺腺癌EMT相關基因篩選及其ceRNA網(wǎng)絡構建*

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

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

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

1 材料與方法

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

1.2方法

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

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

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

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

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

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

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

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

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

2 結 果

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

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

2.2EMT相關基因篩選 將565個差異表達基因導入GenClip3中,結果顯示所有差異表達基因均可在文獻中檢索到,使用GenClip3對差異表達基因進行聚類分析,得到156個EMT相關基因。

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

圖2 EMT相關基因的GO富集圖

圖3 EMT相關基因的KEGG富集圖

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

表1 關鍵基因及degree

圖4 關鍵基因的PPI

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

表2 針對關鍵基因的部分治療性中藥

圖5 中藥與關鍵基因的相關性氣泡圖

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

圖6 關鍵基因的Kaplan-Meier分析結果森林圖

2.7根據(jù)ceRNA假說預測關鍵基因上游miRNA及miRNA上游lncRNA 差異分析結果顯示,IL-6(logFC=-2.64)在腫瘤組織中呈低表達,而MMP9(logFC=2.14)、COL1A1(logFC=3.05)、SPP1(logFC=4.95)、EZH2(logFC=2.70)在腫瘤組織中呈高表達,同理可知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與其靶基因的差異表達趨勢相反,符合ceRNA假說;并且上述miRNA的Log-rankP<0.05,與患者預后相關,因此將其納入下一步研究。由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與靶基因的差異表達趨勢相同,而與上游miRNA差異表達趨勢相反,符合ceRNA假說;并且上述lncRNA的Log-rankP<0.05,與患者預后相關,因此納入下一步研究。見圖7～9。

圖7 mRNA-miRNA-lncRNA相互作用的?；鶊D

圖8 基因差異表達分析的偏差圖

圖9 基因HR的棒棒糖圖

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

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

圖11 EZH2的Kaplan-Meier分析曲線

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

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

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

圖15 GSEC的生存曲線

圖16 ceRNA網(wǎng)絡及各組分之間表達水平的相關性

3 討 論

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

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