基于加權基因共表達網絡分析識別心肌肥厚中的樞紐基因

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

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

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

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

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

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

1.3 基因富集分析

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

1.4 蛋白-蛋白互作網絡構建及樞紐基因識別

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

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

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

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

在取材之前測量大鼠體質量(BW),動物麻醉后,迅速取出大鼠心臟,濾紙吸干,沿房室分隔處分離心房及右心室,保留室間隔及左心室,稱左心室質量(LVW)和全心質量(HW),計算左心室質量指數(shù)(LVWI)=左心室質量(LVW)/體質量(BW)和全心質量指數(shù)(HWI)=全心質量(HW)/體質量(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Ⅱ試劑盒,應用實時定量PCR擴增儀(ABI StepOne Plus)檢測每個擴增循環(huán)后SYBR Green熒光信號水平。PCR擴增條件:94 ℃ 30 s,58 ℃ 15 s,72 ℃ 10 s,共計40個循環(huán)。收集熒光信號,建立熔解曲線,運用熔解曲線檢測擴增產物純度。采用相對定量2-ΔΔCt法對獲得的Ct值進行分析。用Student’st檢驗分析樞紐基因在心肌肥厚與正常對照之間的差異[13]。由上海生工生物工程公司合成引物,序列見表1。

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

1.8 統(tǒng)計學方法

2 結果

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

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

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

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

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

2.3 閾值選擇與網絡拓撲結構分析

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

圖2 軟閾值參數(shù)網絡拓撲結構的分析Figure 2 Analysis of network topology for soft thresholding parameters

2.4 TOM聚類鑒定基因模塊

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

每個分支代表一個基因,下面的每種顏色代表一個網絡模塊

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

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

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

圖B中每一列代表一個組別,每一行對應一個模塊特征基因,顏色表示相關性大小

2.6 關鍵模塊的功能分析與信號通路分析

經過GO富集分析,我們得到了Blue模塊和Turquoise模塊的詳細功能(見圖6A和圖6B),發(fā)現(xiàn)Blue模塊中的基因主要參與肌細胞發(fā)育、表皮細胞增殖、氨基多糖結合、膠原結合、TGF-βg受體結合等分子功能和生物學過程及信號通路,而Turquoise模塊中的基因主要參與免疫系統(tǒng)調節(jié)、刺激反應、炎癥反應、趨化因子和細胞因子介導的炎癥反應信號通路等分子功能和生物學過程,這可能與心肌肥厚的發(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 樞紐基因識別與共表達網絡可視化

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

對Blue模塊和Turquoise模塊內連通性最大的30個基因構建PPI網絡,應用Network Analyzer計算節(jié)點度(degree)

2.8 樞紐基因的驗證

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

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

3 討論

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

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

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

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