基于加權(quán)基因共表達網(wǎng)絡(luò)分析識別心肌肥厚中的樞紐基因

付一樂,陶 亮,彭富華,鄭寧澤,林 青,蔡少儀,王 琴(中山大學(xué)中山醫(yī)學(xué)院藥理學(xué)教研室,廣州 510080;共同第一作者;通訊作者,E-mail:wangqin6@mail.sysu.edu.cn)

心肌肥厚(cardiac hypertrophy,CH)是心肌細胞對多種刺激因素所產(chǎn)生的適應(yīng)性反應(yīng)[1]。這是一個復(fù)雜的病理過程,其病理變化包括心肌細胞肥大、心肌間質(zhì)細胞增殖以及心臟細胞外基質(zhì)改建等多方面的改變,即心肌重構(gòu)[2]。心肌肥厚是許多心臟疾病發(fā)展的一個重要階段,但是心肌肥厚發(fā)生的分子機制非常復(fù)雜,人們對心肌肥厚發(fā)生的確切分子和細胞信號轉(zhuǎn)導(dǎo)機制仍知之甚少[3]。

近年來,隨著分子生物學(xué)和細胞生物學(xué)等技術(shù)的應(yīng)用,有關(guān)心肌肥厚發(fā)生發(fā)展中機制方面的研究取得了一定進展,主要集中在JAK-STAT信號通路,Ca2+介導(dǎo)的信號轉(zhuǎn)導(dǎo)通路、有絲分裂原活化蛋白激酶信號轉(zhuǎn)導(dǎo)通路等[4,5]。但是,僅僅從局部關(guān)注單個或某幾個基因已經(jīng)不能滿足這種具有高度復(fù)雜性的調(diào)控研究。立足于調(diào)控網(wǎng)絡(luò)整體,心肌肥厚中某些基因異常表達且與其他多個基因關(guān)系密切,它們的表達可能在心肌肥厚的發(fā)生發(fā)展過程中扮演重要的角色,這類處于調(diào)控網(wǎng)絡(luò)樞紐位置的基因稱作樞紐基因(hub genes)[6,7]。近年來,隨著基因芯片、RNA-Seq等高通量測序技術(shù)的興起,產(chǎn)生了大規(guī)模組學(xué)數(shù)據(jù),進一步數(shù)據(jù)分析與挖掘?qū)⒂欣谙到y(tǒng)性地研究多個基因表達關(guān)聯(lián)性。因此,根據(jù)基因芯片數(shù)據(jù)來研究心肌肥厚發(fā)生與發(fā)展過程中的樞紐基因具有明顯意義[8]。

在生物信息學(xué)領(lǐng)域中,網(wǎng)絡(luò)分析的應(yīng)用日漸成熟。其中,作為系統(tǒng)生物學(xué)分析方法中極具潛力的一種新方法,加權(quán)基因共表達網(wǎng)絡(luò)分析(weighted gene co-expression network analysis,WGCNA)以其在描述解析分子作用機制和網(wǎng)絡(luò)關(guān)系等方面具有的獨特優(yōu)勢迅速脫穎而出[9]。自開發(fā)以來,WGCNA已成功應(yīng)用于多種生物環(huán)境,并取得了一系列的研究成果。

本研究采用WGCNA進行研究,從大量基因中篩選與心肌肥厚密切相關(guān)的樞紐基因。結(jié)合GO富集分析,進行功能注釋,然后對基因模塊進行KEGG信號通路分析,研究樞紐基因與心肌肥厚發(fā)生與發(fā)展間的關(guān)聯(lián)性,這些工作能夠為某些信號通路的調(diào)控機制提供參考,有利于構(gòu)建評估心肌肥厚早期診斷預(yù)測平臺,降低人群心肌肥厚發(fā)病率及死亡率。有理由相信,隨著基因篩選方法的不斷發(fā)展,WGCNA有望揭示心肌肥厚發(fā)生與發(fā)展的分子機制,并在識別潛在的治療靶點等領(lǐng)域發(fā)揮越來越大的作用。

1 材料與方法

1.1 GEO全基因組表達數(shù)據(jù)

從NCBI的GEO下載心肌肥厚的全基因組表達數(shù)據(jù)GSE36961。該數(shù)據(jù)集發(fā)表于2012年9月,包括106個肥厚型心肌病患者的心臟組織和39個正常對照組的心臟組織,提取mRNA后,用Illumina HumanHT-12 V3.0芯片檢測全基因組基因表達。

1.2 共表達網(wǎng)絡(luò)的構(gòu)建和模塊的識別

利用R-3.6.2軟件運行“WGCNA”包[9]。為降低運算量,計算每個基因的倍數(shù)差異值(取對數(shù)變換),并推斷差異發(fā)生的概率P,當(dāng)倍數(shù)差異值大于預(yù)設(shè)倍數(shù),且篩選P<0.05時即判定為差異基因(即該基因表達無差異的概率小于0.05)。將篩選出的差異基因用于共表達網(wǎng)絡(luò)的構(gòu)建。采用樣本聚類樹方法,根據(jù)聚類圖剔除離群樣本來保證構(gòu)建穩(wěn)定的共表達網(wǎng)絡(luò)。構(gòu)建無尺度網(wǎng)絡(luò)使基因共表達網(wǎng)絡(luò)符合無尺度現(xiàn)象,以無尺度網(wǎng)絡(luò)指數(shù)(R2)=0.9作為滿足無尺度條件的標(biāo)準(zhǔn),同時根據(jù)平均連接度確定軟閾值(β)。利用拓?fù)渲丿B矩陣(TOM)、相異常度矩陣計算基因與基因間的關(guān)聯(lián)程度;對基因構(gòu)建層次聚類樹圖形,采用動態(tài)剪枝法計算基因模塊顏色。計算基因模塊的特征值(ME),引入臨床信息,對ME進行分層聚類并繪制樹狀圖,設(shè)置高度值0.25為分割線,合并相似程度較高的基因模塊,再用剪切后的模塊繪制新的聚類樹和模塊圖。

1.3 基因富集分析

基因本體(GO)及京都基因與基因組百科全書(KEGG)通路富集廣泛應(yīng)用于高通量基因數(shù)據(jù)功能分析。DAVID數(shù)據(jù)庫為研究人員提供了一套全面的功能注釋工具,以了解大量基因背后的生物學(xué)意義[10]。本研究采用DAVID數(shù)據(jù)庫對關(guān)鍵網(wǎng)絡(luò)模塊進行GO富集注釋分析,其中包含生物過程(biological process,BP)、細胞成分(cellular component,CC)、分子功能(molecular function,MF)以及KEGG通路富集分析。

1.4 蛋白-蛋白互作網(wǎng)絡(luò)構(gòu)建及樞紐基因識別

STRING數(shù)據(jù)庫是用于預(yù)測蛋白質(zhì)-蛋白質(zhì)相互作用(protein-protein interaction,PPI)的生物學(xué)數(shù)據(jù)庫[11],可用于檢索經(jīng)實驗驗證或基于文獻預(yù)測的蛋白/基因相互作用。為了進一步探討關(guān)鍵網(wǎng)絡(luò)模塊的生物學(xué)功能,本研究應(yīng)用STRING數(shù)據(jù)庫構(gòu)建靶點PPI網(wǎng)絡(luò)。在數(shù)據(jù)庫中將蛋白種類設(shè)置為“Homo sapiens”(人類)進行操作,最低相互作用閥值設(shè)為高置信度“high confidence≥0.4”,其余參數(shù)均保持默認(rèn)。進一步將獲得的PPI數(shù)據(jù)導(dǎo)入Cytoscape(Version 3.7.2)進行可視化分析,并應(yīng)用Network Analyzer分析模塊進行網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)計算,依據(jù)節(jié)點度(degree)值識別網(wǎng)絡(luò)樞紐基因,其值越大,表明與該節(jié)點相連邊越多,在維持網(wǎng)絡(luò)整體結(jié)構(gòu)上可能占據(jù)重要地位。

1.5 心肌肥厚模型大鼠的建立

選用由北京維通利華實驗動物技術(shù)有限公司提供的健康狀況良好,清潔級、16周齡、雄性WKY大鼠3只為對照組(WKY組),自發(fā)性高血壓大鼠3只為心肌肥厚模型組(SHR組),合格證號:SYXK(粵)2015-0107,飼養(yǎng)于中山大學(xué)SPF級實驗動物中心,飼養(yǎng)溫度20-25 ℃,濕度40%-60%,光照節(jié)律:06∶00-18∶00,保持通風(fēng)良好。

1.6 左心室質(zhì)量指數(shù)和全心質(zhì)量指數(shù)的測定及心臟組織的HE染色

在取材之前測量大鼠體質(zhì)量(BW),動物麻醉后,迅速取出大鼠心臟,濾紙吸干,沿房室分隔處分離心房及右心室,保留室間隔及左心室,稱左心室質(zhì)量(LVW)和全心質(zhì)量(HW),計算左心室質(zhì)量指數(shù)(LVWI)=左心室質(zhì)量(LVW)/體質(zhì)量(BW)和全心質(zhì)量指數(shù)(HWI)=全心質(zhì)量(HW)/體質(zhì)量(BW)[12]。取2/3心臟組織10%福爾馬林固定,石蠟包埋,沿左室長軸常規(guī)切片,行蘇木素-伊紅染色(HE染色),在光鏡下觀察心臟的病理改變[13]。

1.7 RT-qPCR檢測樞紐基因表達

按照TRIzol試劑盒說明書從SHR模型組大鼠和WKY對照組大鼠左心室區(qū)域提取總RNA,用RT-qPCR檢測APOB、CCND1、SOX9、MYC、APP mRNA水平表達。采用SYBR Premix Ex TaqⅡ試劑盒,應(yīng)用實時定量PCR擴增儀(ABI StepOne Plus)檢測每個擴增循環(huán)后SYBR Green熒光信號水平。PCR擴增條件:94 ℃ 30 s,58 ℃ 15 s,72 ℃ 10 s,共計40個循環(huán)。收集熒光信號,建立熔解曲線,運用熔解曲線檢測擴增產(chǎn)物純度。采用相對定量2-ΔΔCt法對獲得的Ct值進行分析。用Student’st檢驗分析樞紐基因在心肌肥厚與正常對照之間的差異[13]。由上海生工生物工程公司合成引物,序列見表1。

表1 引物序列Table 1 Primer sequences for PCR

1.8 統(tǒng)計學(xué)方法

2 結(jié)果

2.1 GEO全基因組表達數(shù)據(jù)基本信息

GSE36961全基因組表達數(shù)據(jù)包括106個模型組與39個對照組,這些數(shù)據(jù)均已經(jīng)過背景校正和標(biāo)準(zhǔn)化處理。為了減輕噪音和減輕計算機的運行負(fù)擔(dān),通過LIMMA包[14]濾過差異并不顯著的基因,最后保留6 321個差異基因(|logFC|>1,FDR<0.05)作為本文的重點研究對象。

2.2 樣本數(shù)據(jù)初篩

采用層次聚類法得到樣本聚類圖,發(fā)現(xiàn)145個樣本表達譜數(shù)據(jù)中存在4個離群的樣本(見圖1),離群樣本經(jīng)過剪枝去除(剪枝高度為50),最后剩余106個心肌肥厚樣本和35個正常組樣本進入后續(xù)的網(wǎng)絡(luò)分析。

采用層次聚類法得到樣本聚類圖,發(fā)現(xiàn)在聚類高度50(紅線標(biāo)記)時聚類樹形圖發(fā)生明顯偏移,剪枝去除4個離群樣本,最后剩余141個樣本

2.3 閾值選擇與網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)分析

選擇合適的領(lǐng)接矩陣權(quán)重參數(shù)βg以盡量滿足無尺度分布。參數(shù)βg從1-20取值,對某節(jié)點連接度的對數(shù)logi與該節(jié)點出現(xiàn)的概率的對數(shù)logp(i)建立線性模型,參數(shù)βg即系數(shù)R的平方,按照候選的參數(shù)β,返回被檢測的網(wǎng)絡(luò)參數(shù)(見圖2)。我們選取了首次接近0.90時的βg值來構(gòu)建基因網(wǎng)絡(luò),最終我們選取β=7,這樣既保證了網(wǎng)絡(luò)接近于無尺度網(wǎng)絡(luò)同時也是使曲線趨于平滑的最小閾值,并且它使得網(wǎng)絡(luò)的平均連接程度不會太小,這有利于網(wǎng)絡(luò)包含足夠的信息。

圖2 軟閾值參數(shù)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的分析Figure 2 Analysis of network topology for soft thresholding parameters

2.4 TOM聚類鑒定基因模塊

選取β=7來計算網(wǎng)絡(luò)拓?fù)渲丿BTOM,利用層次聚類法得到基因聚類樹。根據(jù)動態(tài)分層剪切樹算法來挖掘基因模塊。設(shè)置每個基因模塊中基因數(shù)目最小為30,設(shè)定基因的聚類高度上限為0.995。最終得到11個網(wǎng)絡(luò)模塊,模塊中的基因個數(shù)從34到1 160不等(見圖3)。另外,Grey模塊表示未分配到任何一個模塊的基因集,其中含有122個基因。

每個分支代表一個基因,下面的每種顏色代表一個網(wǎng)絡(luò)模塊

2.5 模塊關(guān)聯(lián)度分析與鑒定關(guān)鍵模塊

根據(jù)拓?fù)渲丿B程度,選擇全部的基因制作熱圖(見圖4)。通過WGCNA,我們識別了心肌肥厚中相關(guān)度高的基因模塊,即共表達模塊代表的基因在功能上具有緊密的關(guān)系。而且,通過層次聚類法對11個網(wǎng)絡(luò)模塊進行網(wǎng)絡(luò)聚類,結(jié)果表明11個網(wǎng)絡(luò)模塊被分成了兩組(見圖5A)。進一步分析發(fā)現(xiàn),在11個網(wǎng)絡(luò)模塊中,Blue模塊和Turquoise模塊與心肌肥厚組具有高相關(guān)性(見圖5B),因此,我們選擇Blue模塊和Turquoise模塊作為心肌肥厚疾病的關(guān)鍵模塊作進一步的功能分析。

不同的顏色代表不同的模塊,圖中間的黃色密度表示不同模塊之間的關(guān)聯(lián)性

圖B中每一列代表一個組別,每一行對應(yīng)一個模塊特征基因,顏色表示相關(guān)性大小

2.6 關(guān)鍵模塊的功能分析與信號通路分析

經(jīng)過GO富集分析,我們得到了Blue模塊和Turquoise模塊的詳細功能(見圖6A和圖6B),發(fā)現(xiàn)Blue模塊中的基因主要參與肌細胞發(fā)育、表皮細胞增殖、氨基多糖結(jié)合、膠原結(jié)合、TGF-βg受體結(jié)合等分子功能和生物學(xué)過程及信號通路,而Turquoise模塊中的基因主要參與免疫系統(tǒng)調(diào)節(jié)、刺激反應(yīng)、炎癥反應(yīng)、趨化因子和細胞因子介導(dǎo)的炎癥反應(yīng)信號通路等分子功能和生物學(xué)過程,這可能與心肌肥厚的發(fā)生有一定的聯(lián)系。

另外,我們還對Blue模塊和Turquoise模塊進行了KEGG信號通路富集分析(見圖6C和圖6D),發(fā)現(xiàn)Blue模塊富集于10條KEGG信號通路,包括TGF-β信號通路、Wnt信號通路、凋亡信號通路等,而Turquoise模塊富集于64條KEGG信號通路,包括AMPK信號通路、PPAR信號通路、腎素-血管緊張素系統(tǒng)、IL-17信號通路、凋亡通路、縫隙連接、TNF信號通路和吞噬體等。

BP:生物過程;CC:細胞組成;MF:分子功能

2.7 樞紐基因識別與共表達網(wǎng)絡(luò)可視化

通過比較模塊內(nèi)連通性和基因重要性可以幫助找到一些與心肌肥厚相關(guān)的樞紐基因。在模塊中與其他基因關(guān)聯(lián)性越多的基因,往往在模塊中的作用越重要,Blue模塊和Turquoise模塊中在右上角的基因既與其他基因相關(guān)性高又對心肌肥厚重要(見圖7A和7B)。因此我們用Cytoscape軟件分別對Blue模塊和Turquoise模塊構(gòu)建蛋白-蛋白相互作用網(wǎng)絡(luò),特別關(guān)注模塊內(nèi)連通性最大的30個基因。結(jié)果發(fā)現(xiàn),Blue模塊和Turquoise模塊均存在多個樞紐基因,如APOB、SOX9、CCND1等(見圖7C和7D)。

對Blue模塊和Turquoise模塊內(nèi)連通性最大的30個基因構(gòu)建PPI網(wǎng)絡(luò),應(yīng)用Network Analyzer計算節(jié)點度(degree)

2.8 樞紐基因的驗證

將3只WKY大鼠作為正常組,3只SHR大鼠作為模型組,通過比較左心室質(zhì)量指數(shù)(LVWI)和全心質(zhì)量指數(shù)(HWI)來觀察SHR模型組的心肌肥厚情況。結(jié)果發(fā)現(xiàn),模型組的LVWI和HWI分別為2.23和3.56,與正常組比較差異均具有統(tǒng)計學(xué)意義(P<0.05,見圖8A),提示SHR大鼠的左心室已出現(xiàn)心肌肥厚。另外,HE染色光鏡下可見:正常組心肌細胞結(jié)構(gòu)清晰、排列整齊,細胞間隙均勻,未見心肌排列纖維增生;與正常組比較,模型組心肌細胞增大、排列紊亂,細胞間隙增寬,可見大量的膠原纖維增生,炎性細胞浸潤在組織周圍,提示SHR大鼠心肌細胞已出現(xiàn)病理損傷。RT-qPCR實驗結(jié)果顯示(見圖8)。與WKY大鼠相比,SHR大鼠心肌組織中APOB、SOX9、CCND1、MYC和APP mRNA水平分別升高1.57倍、1.81倍、2.64倍、1.99倍和2.61倍,差異均具有統(tǒng)計學(xué)意義(P<0.05)。

與WKY組相比,*P<0.05,**P<0.01

3 討論

目前學(xué)術(shù)界普遍認(rèn)為,心肌肥厚發(fā)生與發(fā)展的情況異常復(fù)雜[15],不僅與壓力負(fù)荷、容量負(fù)荷、機械負(fù)荷相關(guān),還與神經(jīng)體液因子密切相關(guān)[16]。心臟為適應(yīng)各種刺激而產(chǎn)生的心肌細胞體積增大和質(zhì)量增加雖然有一定的代償意義,但過度刺激將造成心肌缺血、心臟收縮力下降和輸出功能障礙等,最終可致擴張性心肌病、心力衰竭和猝死[16]。在此過程中,心肌細胞內(nèi)發(fā)生一系列信號轉(zhuǎn)導(dǎo)方面的變化。近年來,人們發(fā)現(xiàn)Ca2+介導(dǎo)的信號轉(zhuǎn)導(dǎo)通路、有絲分裂活化蛋白激酶信號轉(zhuǎn)導(dǎo)通路、Janus激酶/信號轉(zhuǎn)導(dǎo)分子和轉(zhuǎn)錄激活子信號轉(zhuǎn)導(dǎo)通路在心肌肥厚的發(fā)病過程中發(fā)揮重要作用,但大部分的調(diào)控關(guān)系仍處于長期探索中[4]。

本研究首先從NCBI的GEO數(shù)據(jù)庫下載心肌肥厚的表達數(shù)據(jù)集GSE36961(含145個樣本),利用WGCNA模擬基因簇篩選出了與心肌肥厚相關(guān)的基因模塊并對模塊內(nèi)基因進行GO功能注釋和KEGG通路富集分析,這些基因主要富集在免疫系統(tǒng)調(diào)節(jié)、刺激反應(yīng)、炎癥反應(yīng)、趨化因子和細胞因子介導(dǎo)的炎癥反應(yīng)信號通路等分子功能和生物學(xué)過程。需要指出的是,已經(jīng)有研究表明NF-κB信號通路可能在心肌肥厚的發(fā)生發(fā)展中起重要作用[17]。NF-κB信號通路發(fā)揮作用主要是作為介導(dǎo)炎性反應(yīng)的細胞內(nèi)信號轉(zhuǎn)導(dǎo)通路。炎性細胞因子,如:腫瘤壞死因子α(TNF-α)和IL-1等可能是介導(dǎo)NF-κB與心肌肥厚有關(guān)的介質(zhì)。通過采用NF-κB活性抑制劑阻斷NF-κB的激活,可以顯著抑制心肌細胞肥大的進展[17]。另外需要值得注意的是,部分模塊同時富集到凋亡通路、縫隙連接和吞噬體等信號通路,提示這些通路模塊中的基因與心肌肥厚的發(fā)生發(fā)展也可能有著內(nèi)在關(guān)聯(lián),這為心肌肥厚的發(fā)生發(fā)展提供了新的研究思路。

我們還進一步通過Cytoscape軟件分別對Blue模塊和Turquoise模塊構(gòu)建蛋白-蛋白相互作用網(wǎng)絡(luò),將樞紐基因可視化,發(fā)現(xiàn)了與心肌肥厚顯著相關(guān)的基因模塊內(nèi)的樞紐基因,如APOB、CCND1、SOX9、MYC和APP等。通過進一步在自發(fā)性高血壓大鼠心肌肥厚中驗證樞紐基因,發(fā)現(xiàn)APOB、CCND1、SOX9、MYC和APP mRNA水平均顯著上調(diào)。APOB基因編碼產(chǎn)物是乳糜顆粒和低密度脂蛋白的主要載脂蛋白,并且是低密度脂蛋白受體的配體。以往研究發(fā)現(xiàn),APOB異常與動脈粥樣硬化性心血管疾病包括冠心病、腦卒中和周圍動脈疾病密切相關(guān)。但是近年來有研究報道,降低APOB水平的阿托伐他汀藥物在治療心肌肥厚中也能夠發(fā)揮改善作用[18],再結(jié)合我們當(dāng)前的發(fā)現(xiàn),APOB基因的異?？赡茉谛募》屎竦陌l(fā)生發(fā)展中有重要作用,這有待進一步的研究工作來驗證。

綜上所述,加權(quán)基因共表達網(wǎng)絡(luò)分析方法是一個高效的系統(tǒng)生物學(xué)方法,應(yīng)用本方法我們發(fā)現(xiàn)了心肌肥厚的樞紐基因APOB、CCND1、SOX9、MYC和APP。我們當(dāng)前的工作有助于為心肌肥厚發(fā)生的調(diào)控機制的研究以及治療診斷、靶點選擇提供更多的參考信息。