基于生存分析的生物信息學(xué)方法篩選乳腺癌樞紐基因和關(guān)鍵通路

陳 思，劉春良，趙 倩，孫海鵬，劉云霞

上海交通大學(xué)基礎(chǔ)醫(yī)學(xué)院病理生理學(xué)系，細胞分化與凋亡教育部重點實驗室，上海200025

乳腺癌是最為常見的惡性腫瘤之一，也是女性中發(fā)病率最高的惡性腫瘤[1]。根據(jù)雌激素受體（estrogen receptor，ER）、孕激素受體（progesterone receptor，PR）和人類表皮生長因子受體2（human epidermal growth factor receptor-2，HER2）的表達情況，可以將乳腺癌分為4種亞型（Luminal A、Luminal B、HER2陽性和三陰乳腺癌）。根據(jù)不同乳腺癌亞型的生物學(xué)特征和臨床病理分期制定相應(yīng)的個體化治療策略，可以使乳腺癌5年生存率達90%以上[2]。然而，乳腺癌的復(fù)發(fā)和轉(zhuǎn)移仍舊是一大難題。除此之外，一些乳腺癌亞型，例如三陰乳腺癌，由于缺乏有效的治療靶點，一直以來是臨床治療的一個瓶頸。因此，明確乳腺癌發(fā)生發(fā)展相關(guān)的關(guān)鍵基因和通路，有助于認識乳腺癌潛在的發(fā)病機制，或?qū)榕R床尋找更多的診斷和治療靶點提供參考。

自基因芯片技術(shù)和高通量測序技術(shù)問世以來，生物信息學(xué)迅速發(fā)展，目前已發(fā)現(xiàn)了許多疾病的生物學(xué)標志物[3-4]。從信息豐富的公共數(shù)據(jù)庫如GEO數(shù)據(jù)庫（Gene Expression Omnibus，基因表達匯編）和TCGA數(shù)據(jù)庫（The Cancer Genome Atlas，癌癥基因組圖譜）中可獲得基因表達數(shù)據(jù)，通過生物信息學(xué)方法，對數(shù)據(jù)庫中的基因表達數(shù)據(jù)進行聚類分析、統(tǒng)計分析、通路分析和可視化作圖等，能夠預(yù)測基因的功能以及基因間的相互作用，了解疾病基因?qū)用娴陌l(fā)病機制，發(fā)現(xiàn)潛在的生物學(xué)標志物，從而為疾病的分子靶向藥物研發(fā)和精準治療提供理論依據(jù)。本研究將乳腺癌基因表達數(shù)據(jù)與臨床生存分析相結(jié)合，以篩選樞紐基因和關(guān)鍵信號通路；基于生存期篩選出來的基因可能更具有臨床意義，或能為乳腺癌的診斷和治療提供新的思路。

1 材料與方法

1.1 數(shù)據(jù)采集

為了獲取乳腺癌的基因表達數(shù)據(jù)集，本研究在GEO數(shù)據(jù)庫下載了3個乳腺癌數(shù)據(jù)集（http://www.ncbi.nlm.nih.gov/geo/），分別為GSE54002、GSE29431和GSE61304，這3個數(shù)據(jù)集都是基于GPL570平臺。GSE54002包含417例乳腺癌樣本和16例正常樣本。GSE29431包含了54例乳腺癌樣本和12例正常樣本：54例乳腺癌樣本中有15例HER2免疫組織化學(xué)評分為3+，且伴有HER2基因擴增；26例評分為2+，其中13例伴有HER2基因擴增，13例不伴有HER2基因擴增；13例評分為0/1+，且不伴有HER2基因擴增。GSE61304包含了58例乳腺癌樣本和4例正常樣本：58例乳腺癌樣本中18例為ER+PR+，19例為ER-PR-，4例為ER+PR-，1例為ER-PR+，其他16例樣本未說明。

1.2 篩選差異表達基因

在Rstudio軟件中（版本3.4.0）加載并下載Bioconductor網(wǎng)站上軟件包來分析上述3個乳腺癌數(shù)據(jù)集。首先使用affy包導(dǎo)入CEL文件，使用simpleaffy包評估微陣列數(shù)據(jù)質(zhì)量[5]，gcrma包中的RMA算法預(yù)處理原始數(shù)據(jù)[6]，genefilter包過濾非特異性結(jié)合的探針和數(shù)據(jù)質(zhì)量低的探針，limma包進行差異基因表達的統(tǒng)計學(xué)驗證[7]。基因表達變化倍數(shù)的對數(shù)值的絕對值（|log2FC|） ＞1且P＜0.05認定為差異表達基因。最后，為了提高差異表達基因的穩(wěn)健性，使用Funrich軟件（版本3.1.3）獲得3個數(shù)據(jù)集中均上調(diào)或下調(diào)的基因，用于下一步的分析。

1.3 篩選與總體生存期相關(guān)的差異表達基因

為了分析差異表達基因?qū)θ橄侔┗颊呖傮w生存期的影響，在Kaplan-Meier plotter數(shù)據(jù)庫（www.kmplot.com）中根據(jù)基因的中位表達值將患者樣本分為高表達組和低表達組。使用默認參數(shù)，計算每個基因高表達組和低表達組的中位生存期；若log-rank P＜0.05，則該基因被視為與總體生存期相關(guān)的差異表達基因。

1.4 分析與總體生存期相關(guān)的差異表達基因的功能

基因本體數(shù)據(jù)庫（Gene Ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genome，KEGG）通路富集分析被廣泛用于識別基因的功能和通路。本研究使用Rstudio軟件中clustProfiler包對與乳腺癌患者總體生存期相關(guān)的差異表達基因進行GO分析和KEGG分析，P＜0.05認為具有統(tǒng)計學(xué)意義[8]。

1.5 蛋白 - 蛋白相互作用網(wǎng)絡(luò)構(gòu)建及篩選樞紐基因

蛋白與蛋白相互作用在調(diào)節(jié)生物學(xué)過程中起著至關(guān)重要的作用。這種關(guān)系可以通過蛋白 - 蛋白相互作用（proteinprotein interaction，PPI）網(wǎng)絡(luò)表示，每個節(jié)點代表一個蛋白，邊代表蛋白質(zhì)之間的相互作用。緊密相連的區(qū)域可以作為富集功能群。STRING數(shù)據(jù)庫（https://string-db.org/）包含了豐富的蛋白質(zhì)之間相互作用的信息[9]。為了評估與總體生存期相關(guān)的差異表達基因之間的相互關(guān)系，將差異表達基因列表導(dǎo)入STRING數(shù)據(jù)庫，并設(shè)定信度為0.4；接著將PPI表格數(shù)據(jù)導(dǎo)入Cytoscape軟件中構(gòu)建PPI網(wǎng)絡(luò)，使用軟件中的插入式分子復(fù)合物檢測（MCODE評分＞4分，節(jié)點數(shù)＞5個）篩選出PPI網(wǎng)絡(luò)中的樞紐模塊，最后通過CytoHubba插件計算網(wǎng)絡(luò)中每一個基因的最大團中心性（maximal clique centrality，MCC）分數(shù)，將得分前10的基因作為樞紐基因[10-12]。

1.6 樞紐基因驗證

1.6.1 數(shù)據(jù)庫驗證樞紐基因表達 利用Oncomine數(shù)據(jù)庫（https://www.oncomine.org/resource/main.html/）和人類蛋白質(zhì)圖譜（Human Protein Atlas，HPA）數(shù)據(jù)庫（https://www.proteinatlas.org/）對樞紐基因在乳腺癌腫瘤組織和正常組織間的mRNA和蛋白水平的表達進行驗證[13]。

1.6.2 RNA抽提及實時熒光定量PCR 人乳腺癌細胞MDA-MB-231和人正常乳腺上皮細胞MCF-10A（購自中國科學(xué)院干細胞庫）的總RNA抽提按照Invitrogen公司TRIzol試劑盒提供的方法。反轉(zhuǎn)錄以1 000 ng RNA為模板，按TaKaRa公司M-MLV Reverse Transcription Kit試劑盒說明配制反轉(zhuǎn)錄反應(yīng)液，進行反轉(zhuǎn)錄反應(yīng)。實時熒光定量PCR（quantitative real-time PCR，qPCR）使用ABI公司7500 Real-Time PCR System試劑盒。每份采用10 μL體系，cDNA模板1 μL，每份樣品做3個復(fù)孔，求平均值；以18S rRNA為內(nèi)參，引物序列見表1。

表1 qPCR引物序列Tab 1 Primer sequences for qPCR

1.7 統(tǒng)計學(xué)分析

應(yīng)用OriginPro 2017C和Adobe Illustraor CS6軟件進行統(tǒng)計學(xué)分析和作圖，qPCR數(shù)值用x—±s表示，組間比較采用獨立樣本t檢驗。P＜0.05認為差異具有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 篩選差異表達基因

在|log2FC|＞1且P＜0.05的篩選條件下，從GSE54002中得到差異表達基因3 389個，其中上調(diào)基因1 561個，下調(diào)基因1 828個；GSE29431中得到差異表達基因3 660個，其中上調(diào)基因1 097個，下調(diào)基因2 563個；GSE61304中得到差異表達基因1 828個，其中上調(diào)基因821個，下調(diào)基因1 007個。然后用Funrich軟件的Vene圖篩選得到了211個共同上調(diào)和374個共同下調(diào)的差異表達基因，如圖1所示。

圖1 GSE29431、GSE54002和GSE61304中篩選得到的差異表達基因的Venn圖Fig 1 Venn diagrams of the differentially expressed genes screened in GSE29431, GSE54002 and GSE61304

2.2 乳腺癌中與總體生存期相關(guān)的差異表達基因

為了篩選與總體生存期相關(guān)的基因，使用Kaplan-Meier plotter數(shù)據(jù)庫對這585個差異表達基因進行總體生存分析，計算每個基因的高表達組和低表達組的中位生存期和log-rankP值。結(jié)果顯示，有262個基因的高表達組和低表達組的總體生存期之間的差異有統(tǒng)計學(xué)意義。表2列出了與生存相關(guān)最顯著的前10個差異表達基因。這262個與總體生存期相關(guān)的差異表達基因?qū)⒂糜谙乱徊降幕蚬δ芊治觥?/p>

表2 與總體生存期相關(guān)最顯著的前10個差異表達基因Tab 2 List of the top 10 overall survival-related differentially expressed genes

2.3 乳腺癌中與總體生存期相關(guān)的差異表達基因的功能分析

為進一步了解這些基因的功能，利用Rstudio軟件中的clusterProfiler包對得到的生存期相關(guān)的差異表達基因進行功能分析，P＜0.05認為具有統(tǒng)計學(xué)意義。GO功能分析結(jié)果顯示，這些基因的功能主要與細胞分裂（cell division）、核分裂（nuclear division）、細胞器分裂（organelle fission）、細胞周期的調(diào)控（regulation of mitotic cell cycle）、負向調(diào)控細胞有絲分裂（mitotic nuclear division）以及染色體分離（chromosome segregation）等生物學(xué)過程有關(guān)（圖2A）。KEGG通路分析顯示這些基因參與細胞周期（cell cycle）、人T細胞白血病病毒1感染（human T-cell leukemia virus 1 infection）、FoxO信號通路（FoxO signaling pathway）、卵母細胞減數(shù)分裂（oocyte meiosis）、細胞衰老（cellular senescence）以及孕酮介導(dǎo)的卵母細胞成熟（progesterone-mediated oocyte maturation）等信號通路（圖2B）。

圖2 與總體生存期相關(guān)的差異表達基因GO功能分析和KEGG通路分析Fig 2 GO function analysis and KEGG pathway analysis of the overall survival-related differentially expressed genes

2.4 PPI網(wǎng)絡(luò)構(gòu)建及模塊分析

將與總體生存期相關(guān)的差異表達基因列表上傳至STRING，并設(shè)定信度0.4作為判斷相互作用是否有意義的標準，構(gòu)建了PPI網(wǎng)絡(luò)（圖3）。Cytoscape軟件根據(jù)MCODE評分排序篩選出了3個樞紐模塊（圖4），并對樞紐模塊中的基因進行GO功能分析，結(jié)果顯示模塊1和模塊2的基因主要富集在細胞分裂（cell division）、細胞核分裂（nuclear division）和細胞器分裂（organelle fission）等生物學(xué)過程，模塊3的基因主要集中在血小板脫顆粒（platelet degranulation）、負向調(diào)控線粒體細胞色素c的釋 放（negative regulation of release of cytochrome c from mitochondria）等生物學(xué)過程（表3）。

圖3 與總體生存期相關(guān)的差異表達基因的PPI網(wǎng)絡(luò)Fig 3 PPI network of overall survival-related differentially expressed genes

圖4 MCODE評分排序篩選出的3個樞紐模塊Fig 4 Three modules selected by MCODE scoring sorting

表3 3個模塊基因的GO功能分析Tab 3 GO functional analysis of the differentially expressed genes in three modules

2.5 乳腺癌中樞紐基因的鑒定

樞紐基因是一類在生物學(xué)過程中發(fā)揮至關(guān)重要作用的基因，在相關(guān)通路中其他非樞紐基因的調(diào)控往往要受到這類基因的影響，因此樞紐基因有可能成為乳腺癌的生物學(xué)標志物和治療靶標。使用Cytoscape中的插件Cytohubba，通過MCC法得到了得分前10的樞紐基因，分別是NDC80、BUB1、CDCA8、BUB1B、BIRC5、CCNB1、KIF2C、CENPF、MAD2L1和CDC20（圖5）。

圖5 通過MCC法得到的10個樞紐基因及其相互作用Fig 5 Ten hub genes and their interactions by MCC

2.6 樞紐基因的驗證

通過Oncomine數(shù)據(jù)庫和HPA數(shù)據(jù)庫對篩選得到的10個樞紐基因進行表達驗證。Oncomine數(shù)據(jù)庫結(jié)果顯示，BIRC5、CDC20、NDC80、CENPF、MAD2L1、CDCA8、KIF2C、BUB1、CCNB1和BUB1B的mRNA水平在乳腺癌組織中明顯上調(diào)（圖6）。HPA數(shù)據(jù)庫免疫組織化學(xué)檢測結(jié)果顯示，腫瘤 組 織 中CDCA8、BIRC5、CDC20、CENPF、MAD2L1和CCNB1的蛋白表達水平高于正常乳腺組織，另外4個基因NDC80、KIF2C、BUB1、BUB1B未被HPA數(shù)據(jù)庫收錄。以上結(jié)果提示，篩選得到的樞紐基因具有較好的穩(wěn)健性。

圖6 通過Oncomine數(shù)據(jù)庫驗證樞紐基因轉(zhuǎn)錄水平的表達Fig 6 Validation of the expression of hub genes at mRNA level using Oncomine database

此外，通過提取人乳腺癌細胞MDA-MB-231和人正常乳腺上皮細胞MCF-10A的RNA進行反轉(zhuǎn)錄和qPCR，驗證了這些樞紐基因的表達水平。如圖7結(jié)果所示，BIRC5、NDC80、MAD2L1、CENPF、CCNB1、BUB1、BUB1B、KIF2C、CDC20和CDCA8在乳腺癌細胞中的表達均高于人正常乳腺上皮細胞；同時，雖然2個樞紐基因MAD2L1與BUB1的表達差異無統(tǒng)計學(xué)意義，但乳腺癌細胞中水平仍略高于正常乳腺上皮細胞。以上結(jié)果在細胞水平驗證了各樞紐基因表達。

圖7 qPCR檢測乳腺癌細胞樞紐基因的表達Fig 7 qPCR analysis of hub genes in breast cancer cells

3 討論

本研究利用基于生存分析的生物信息學(xué)方法，將基因表達數(shù)據(jù)與臨床生存分析相結(jié)合，篩選乳腺癌中的樞紐基因和關(guān)鍵通路。首先，我們從GEO數(shù)據(jù)庫下載了3個乳腺癌組織和癌旁組織的基因表達數(shù)據(jù)集，篩選出了585個差異表達基因，其中上調(diào)的基因有211個，下調(diào)的基因有374個。然后用Kaplan-Meier plotter數(shù)據(jù)庫篩選出了262個與乳腺癌患者總體生存期相關(guān)的差異表達基因。接著對這262個基因進行了GO功能分析，結(jié)果顯示這些基因主要與細胞分裂、細胞周期的調(diào)控以及染色體分離等生物學(xué)過程相關(guān)；KEGG通路分析結(jié)果表明，這些與總體生存期相關(guān)的差異表達基因主要富集在細胞周期、FoxO信號通路和卵母細胞減數(shù)分裂等通路上。此外我們發(fā)現(xiàn)人類T細胞白血病病毒1感染信號通路在乳腺癌中失調(diào)，但目前并無關(guān)于人類T細胞白血病病毒1感染信號通路與乳腺癌的報道；由于其在乳腺癌中的分子機制仍不明確，因此需要進一步研究。以上篩選出來的通路可以為乳腺癌發(fā)生發(fā)展的分子機制研究提供依據(jù)。

此外，我們通過構(gòu)建PPI網(wǎng)絡(luò)篩選出了10個乳腺癌的樞紐基因，分別為BIRC5、CDC20、NDC80、CENPF、MAD2L1、CDCA8、KIF2C、BUB1、CCNB1和BUB1B，經(jīng)Oncomine數(shù)據(jù)庫和HPA數(shù)據(jù)庫以及qPCR驗證，它們在乳腺癌中均高表達，并與乳腺癌患者較差的生存期相關(guān)。通過查閱文獻發(fā)現(xiàn)這些基因主要參與了細胞有絲分裂染色體的分離和細胞周期的調(diào)控。我們還嘗試在Oncomine數(shù)據(jù)庫中研究這些樞紐基因的表達是否與乳腺癌分子分型、分期及惡性程度有關(guān)；但由于數(shù)據(jù)庫中數(shù)據(jù)的一些局限，無法直接分析樞紐基因在不同分子分型乳腺癌的表達情況，因此無法確認樞紐基因的表達與分子分型的相關(guān)性。但我們發(fā)現(xiàn)樞紐基因在ER陰性乳腺癌和PR陰性乳腺癌的表達量分別高于ER陽性乳腺癌和PR陽性乳腺癌的表達量；在HER2陽性或陰性的乳腺癌中，這些基因的表達卻沒有明顯差別（數(shù)據(jù)未展示），提示樞紐基因的表達可能與分子分型有關(guān)。通過分析樞紐基因在不同分期的乳腺癌中的表達情況，我們發(fā)現(xiàn)這些基因在ⅠB期的表達量明顯低于其他分期的表達量，而BUB1B在Ⅲ期的表達量顯著高于其他分期（數(shù)據(jù)未展示）。MAD2L1、CDC20、CENPF、KIF2C、CCNB1、NDC80、BUB1、CDCA8和KIF2C在各分期的表達量沒有明顯差異。上述發(fā)現(xiàn)說明樞紐基因可能可以作為特定的分子分型和分期治療的判斷依據(jù)。

有 研 究 報 道，NDC80、MAD2L1、CDC20、BUB1、BIRC5和CCNB1在乳腺癌的發(fā)生發(fā)展中發(fā)揮了重要功能。NDC80基因編碼的NDC80蛋白參與構(gòu)成微管與動粒連接復(fù)合體，為染色體正常分離所必需，在細胞有絲分裂過程中起著至關(guān)重要的作用；NDC80的異常表達會造成染色體的異常分離，從而使染色體不穩(wěn)定，最終導(dǎo)致腫瘤的發(fā)生[14]。多項研究[15-16]表明NDC80在多種腫瘤中高表達并參與腫瘤的發(fā)生發(fā)展。NDC80抑制劑TAI-95可以在體外和體內(nèi)抑制乳腺癌腫瘤生長，NDC80有望成為乳腺癌的治療靶點[17]。MAD2L1和CDC20蛋白是紡錘體檢查點（spindle assembly checkpoint，SAC）復(fù)合體的組成成分，SAC負責(zé)確保姐妹染色單體的動粒和紡錘體結(jié)合并準確分離到2個子細胞中[18]。研究[19]顯示通過shRNA敲除MAD2L1可以抑制乳腺癌細胞的生長和侵襲能力。CDC20的高表達與乳腺癌患者較差的預(yù)后相關(guān)，且由它介導(dǎo)的BTG3相關(guān)核蛋白 （BTG3 associated nuclear protein，SMAR1）降解可以促進乳腺癌細胞遷移和侵襲[20-21]。BUB1是一類絲氨酸/蘇氨酸蛋白激酶，在有絲分裂和DNA損傷應(yīng)答中發(fā)揮著重要的功能[22]。Han等[23]發(fā)現(xiàn)，敲低BUB1會降低腫瘤干細胞的潛能，BUB1可能成為針對乳腺癌干細胞的治療靶點。BIRC5編碼的蛋白是凋亡抑制（inhibitor of apoptosis，IAP）家族的成員之一，BIRC5具有明顯的抑制細胞凋亡的作用。Li等[24]認為BIRC5與乳腺癌患者進展及不良預(yù)后相關(guān)；研究[25-27]發(fā)現(xiàn)，BIRC5是miR-485-5p、ZIC家族成員1（Zic family member 1，ZIC1）和半乳凝素1（galectin-1）的靶點，可能在乳腺癌中發(fā)揮重要作用。CCNB1蛋白是調(diào)控細胞周期G2/M過渡階段所必需的；研究[28]表明，CCNB1蛋白與侵襲程度（腫瘤分級、腫瘤體積和淋巴結(jié)狀態(tài)）相關(guān)，是乳腺癌的預(yù)后因素之一。

目前KIF2C、CENPF、CDCA8和BUB1B這4個樞紐基因與乳腺癌相關(guān)的研究還比較少，具體分子機制有待進一步研究。KIF2C和CENPF蛋白參與有絲分裂的染色體分離。研究[29]報道，KIF2C在多種上皮型腫瘤中高表達，并與腫瘤的分級、分期和預(yù)后相關(guān)。CENPF基因編碼一種與著絲粒 - 動粒復(fù)合體相關(guān)的蛋白，其表達水平在有絲分裂前期增加，并在后期開始時發(fā)生蛋白水解[30-31]；CENPF的上調(diào)與乳腺癌的不良預(yù)后相關(guān)，可能是有臨床意義的治療靶點[32]。CDCA8是細胞周期分裂相關(guān)蛋白家族的一員，此蛋白受細胞周期調(diào)控，對染色質(zhì)誘導(dǎo)的微管穩(wěn)定和紡錘體形成是必需的；CDCA8的高表達與乳腺癌患者的低生存率和較差預(yù)后相關(guān)[33]。BUB1B基因編碼的蛋白是SAC復(fù)合體的重要成員，BUB1B表達增加是惡性程度高的乳腺癌的特征之一[34]。

綜上所述，本研究基于與總生存期相關(guān)的差異表達基因，篩選了乳腺癌中的樞紐基因和關(guān)鍵通路；基于生存期的樞紐基因可能更具有臨床意義，但這些樞紐基因的功能仍需進一步研究。