利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析識(shí)別膽道閉鎖相關(guān)的樞紐基因

周紅，吳梅，李鑫

（1.宜賓市第一人民醫(yī)院兒科，宜賓 644000;2.烏魯木齊兒童醫(yī)院普外科，烏魯木齊 830000）

膽道閉鎖（biliary atresia，BA）是一種嬰幼兒期常見的極為嚴(yán)重的疾病，特點(diǎn)為發(fā)生于生后3 個(gè)月內(nèi)的部分或者全部肝外膽道完全性纖維化梗阻，如果不治療，將發(fā)展為肝硬化、肝衰竭，不經(jīng)治療的平均生存期在1 年左右，主要的治療方式是肝門腸吻合術(shù)（hepatoportoenterostomy，HPE）和肝移植[1-4]。即使在出生后45 d 內(nèi)進(jìn)行HPE，2 歲時(shí)未進(jìn)行肝移植的生存率僅為65.5%[5]。BA 面臨早期診斷非常困難、手術(shù)治療效果不佳、肝移植后并發(fā)癥和免疫抑制管理困難等多種挑戰(zhàn)。目前BA 的發(fā)病機(jī)制、病因仍不清楚[6-7]，近年來專注于鑒定BA 相關(guān)的遺傳和免疫因子作用的研究成為熱點(diǎn)。且BA 是新生兒膽汁淤積最常見的原因，BA 與非BA 的膽汁淤積性疾病在臨床表現(xiàn)和實(shí)驗(yàn)室檢查上非常相似，除手術(shù)探查外的方法難以鑒別[8]，因此尋求一種無創(chuàng)的方法以鑒別這兩種疾病非常有必要。

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）是一種分析多個(gè)樣本基因表達(dá)模式的分析方法，將具有類似表達(dá)模式的基因整合到同一的模塊，并分析模塊與特定性狀或表型之間的關(guān)聯(lián)關(guān)系[9-11]。本研究擬利用WGCNA篩選出與BA 相關(guān)的基因模塊及樞紐基因，以期為BA的發(fā)病機(jī)制及臨床診治帶來新的理解。

1 材料與方法

1.1 數(shù)據(jù)集的收集 從美國(guó)國(guó)立生物信息技術(shù)中心（National Center for Biotechnology Information，NCBI）的基因表達(dá)綜合數(shù)據(jù)庫(kù)（Gene Expression Omnibus，GEO）下載GSE46960數(shù)據(jù)集，芯片測(cè)序平臺(tái)為GPL6244。從該數(shù)據(jù)集中選取64 例診斷為BA 的肝臟樣本（BA 組），14 例非BA 的膽汁淤積性疾病肝臟樣本（Non_BA 組）以及7 名已故無肝臟疾病捐贈(zèng)者的肝臟樣本（Normol 組）數(shù)據(jù)進(jìn)行分析。采用R 軟件進(jìn)行數(shù)據(jù)的處理。

1.2 加權(quán)共表達(dá)基因網(wǎng)絡(luò)的構(gòu)建 以表達(dá)量為篩選條件，選擇平均表達(dá)量最高的5 000 個(gè)基因在R 軟件的WGCNA 包構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。為確保無尺度網(wǎng)絡(luò)分別計(jì)算所有基因之間的成對(duì)Pearson 相關(guān)系數(shù)，設(shè)置閾值進(jìn)行篩選，以將成對(duì)相關(guān)矩陣轉(zhuǎn)換為鄰近相關(guān)矩陣。采用動(dòng)態(tài)混合剪切樹算法標(biāo)準(zhǔn)，設(shè)置每個(gè)基因模塊最少基因數(shù)目為30，并且依次計(jì)算每個(gè)模塊的特征向量值，然后對(duì)模塊進(jìn)行聚類分析，將距離較近的模塊合并為新的模塊。

1.3 模塊特征相關(guān)性的計(jì)算 通過WGCNA 算法計(jì)算模塊基因與疾病分組表型之間的相關(guān)性，相關(guān)性的強(qiáng)度通過熱圖反應(yīng)出來。當(dāng)P＜0.05 時(shí)，認(rèn)為單個(gè)模塊與表型顯著相關(guān)。選擇與BA 相關(guān)系數(shù)最高的模塊作為關(guān)鍵模塊。

1.4 樞紐基因的篩選 計(jì)算每個(gè)共表達(dá)模塊與基因特征值的皮爾森相關(guān)系數(shù)。選取模塊身份（module membership，MM）＞0.8 且基因顯著性（gene significance，GS）＞0.65 的基因?yàn)闃屑~基因。

1.5 關(guān)鍵模塊的富集分析 提取關(guān)鍵模塊中的基因，利用R 軟件ClusterProfiler 包進(jìn)行基因本體論（Gene Ontology，GO）功能注釋和《京都基因和基因組百科全書》（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析。P＜0.05 為差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 樣本數(shù)據(jù)篩選 經(jīng)過歸一化和基因ID 轉(zhuǎn)換生成18 762 個(gè)基因和85 個(gè)樣本的表達(dá)矩陣。選取平均表達(dá)量最高的5 000 個(gè)基因構(gòu)建基因共表達(dá)模塊。數(shù)據(jù)集處理后進(jìn)行離群值檢測(cè)，分析結(jié)果示該數(shù)據(jù)集中無明顯離群值，故直接進(jìn)行下一步分析，將基因聚類模塊與臨床分組表型進(jìn)行分析（圖1）。

圖1 基于歐幾里得距離的樣本聚類樹圖Fig 1 The sample cluster tree diagram based on Euclidean distance

2.2 軟閾值計(jì)算 利用WGCNA 構(gòu)建共表達(dá)網(wǎng)絡(luò)需先確定軟閾值（soft-thresholding power），本研究設(shè)定1~20，對(duì)每個(gè)軟閾值計(jì)算無標(biāo)度拓?fù)鋽M合指數(shù)R2（scale free topology fitting index R2）和平均連通性（mean connectivity）。軟閾值對(duì)應(yīng)的R2＞0.8 且平均連通性接近零，說明此網(wǎng)絡(luò)符合無尺度條件。因此，本研究選取軟閾值β = 9（圖2）。

圖2 軟閾值參數(shù)的確定Fig 2 Determination of soft threshold parameters

2.3 構(gòu)建基因共表達(dá)模塊 使用動(dòng)態(tài)混合切割構(gòu)建層次聚類樹，樹上的每一片葉子代表一個(gè)基因，具有相似表達(dá)數(shù)據(jù)的基因靠在一起，形成樹的一個(gè)分支，代表一個(gè)基因模塊，生成了20 個(gè)模塊（圖3）。

圖3 模塊構(gòu)建Fig 3 Construction of modules

2.4 共表達(dá)模塊之間的相關(guān)性分析 將層次聚類得到的基因模塊進(jìn)行相關(guān)性分析，觀察模塊與模塊之間的關(guān)聯(lián)，顏色越深則表示模塊之間相關(guān)性越高（圖4）。

圖4 特征基因樹圖和特征基因鄰接圖Fig 4 Characteristic gene tree diagram and characteristic gene adjacency diagram

2.5 共表達(dá)模塊的關(guān)聯(lián)分析 根據(jù)拓?fù)渲丿B程度，選擇全部的基因制作熱圖（圖5）。

圖5 共表達(dá)網(wǎng)絡(luò)的熱圖可視化Fig 5 Thermal graph visualization of co-expression network

2.6 計(jì)算BA 疾病相關(guān)的模塊相關(guān)性 對(duì)于每一個(gè)模塊，計(jì)算基因表達(dá)與不同疾病組之間的相關(guān)性。與BA 相關(guān)性最大的模塊是黃色模塊，相關(guān)系數(shù)為0.69（P＜0.05），同時(shí)黃色模塊與非BA 的膽汁淤積性疾病及正常組是負(fù)相關(guān)關(guān)系（圖6），故將黃色模塊識(shí)別為關(guān)鍵模塊。

圖6 模塊基因與臨床表型相關(guān)性圖Fig 6 Correlation between modular genes and clinical phenotypes

2.7 識(shí)別關(guān)鍵模塊內(nèi)的樞紐基因 分析黃色模塊內(nèi)基因連通性（MM）及基因顯著性（GS）之間的相關(guān)性，二者相關(guān)性良好，呈明顯線性相關(guān)（圖7）。以模塊身份MM＞0.8 且基因顯著性GS＞0.65 作為篩選條件，共篩選出16 個(gè)樞紐基因，包括基質(zhì)金屬蛋白酶7（matrix metallopeptidase 7，MMP7）、重組人分泌型磷蛋白1（secreted phosphoprotein 1，SPP1）、多功能蛋白聚糖（versican，VCAN）、含V-Set 域T 細(xì)胞激活抑制因子1（V-set domain containing T cell activation inhibitor，VTCN1）、重組人光蛋白聚酶（lumican，LUM）、G-蛋白信號(hào)轉(zhuǎn)導(dǎo)調(diào)節(jié)子4（regulator of G protein signaling 4，RGS4）、上皮細(xì)胞黏附分子（epithelial cell adhesion molecule，EPCAM）、超長(zhǎng)鏈脂肪酸延伸酶7（elongation of very long fatty acid elongase 7，ELOVL7）、MET 轉(zhuǎn)錄調(diào)節(jié)因子MACC1（MET transcriptional regulator MACC1，MACC1）、層黏連蛋白γ2（laminin subunit gamma 2，LAMC2）、染色體12開放閱讀框75（chromosome 12 open reading frame 75，C12orf75）、IV 型膠原蛋白基因α1（collagen type IV alpha 1 chain，COL4A1）、波形纖維蛋白（vimentin，VIM）、跨膜蛋白（transmembraneprotein156，TMEM156）、胸腺素β10（thymosin beta 10，TMSB10）、肌球蛋白重鏈9（myosin heavy chain 9，MYH9）（表1）。

表1 黃色模塊中的樞紐基因Tab1 Hub gene in yellow module

圖7 黃色模塊基因相關(guān)性散點(diǎn)圖Fig 7 Scatter plot of module eigengenes in the yellow module

2.8 模塊的富集分析 黃色模塊GO 富集結(jié)果表明，生物學(xué)過程（biological process，BP）主要富集在細(xì)胞外基質(zhì)組織和細(xì)胞外結(jié)構(gòu)組織，細(xì)胞成分（cellular component，CC）主要富集于含膠原的細(xì)胞外基質(zhì)，分子功能（molecular function，MF）主要富集于細(xì)胞黏附分子。黃色模塊KEGG 富集結(jié)果示，BA 的發(fā)生可能與色氨酸代謝、甘油磷脂代謝、過氧化物酶體增殖物激活受體（peroxisome proliferator-activated receptor，PPAR）信號(hào)通路等密切相關(guān)。見圖8。

圖8 GO 和KEGG 富集氣泡圖Fig 8 Bubble diagrams of GO and KEGG enrichment

3 討論

BA 是發(fā)生在嬰幼兒期的一種嚴(yán)重的先天性疾病，由于肝門附近的膽道系統(tǒng)狹窄、閉鎖或缺如，結(jié)合膽紅素排入腸道受阻，最終導(dǎo)致肝臟膽汁性硬化、肝衰竭，甚至死亡，生后1 個(gè)月內(nèi)行HPE 治療肝臟存活率明顯增高，但診斷該疾病需排除非BA 的膽汁淤積性疾病，如果不能區(qū)分BA 與非BA 的膽汁淤積性疾病可能導(dǎo)致手術(shù)治療的延遲[1-2]。既往的研究通過差異分析的方法鑒定BA 的特征性基因，可能忽略一些具有重要功能但差異表達(dá)不明顯的基因，也可能混雜與非BA 的膽汁淤積性疾病的相關(guān)基因。WGCNA 是一種無監(jiān)督的層次聚類方法，能夠識(shí)別與表型相關(guān)的“基因模塊”，即可以整合在基礎(chǔ)生物學(xué)途徑中的高度互連的基因。與基因差異分析相比，WGCNA 劃分的基因模塊具有明顯的生物學(xué)意義。結(jié)合臨床信息分析，可以確定與疾病發(fā)病機(jī)制相關(guān)的重要模塊和潛在樞紐基因[12]。為進(jìn)一步分析BA 的發(fā)生、發(fā)展機(jī)制，并將非BA 的膽汁淤積性疾病加以鑒別，本研究采用WGCNA 分析BA、非BA 的膽汁淤積性疾病和無肝臟疾病肝臟轉(zhuǎn)錄組數(shù)據(jù)。

本研究構(gòu)建了20 個(gè)基因共表達(dá)模塊，鑒定了1 個(gè)與BA 高度相關(guān)的基因模塊即黃色模塊，其與非BA 的膽汁淤積性疾病和無肝臟疾病無明顯相關(guān)性。本研究對(duì)關(guān)鍵模塊進(jìn)行GO 及KEGG 富集，發(fā)現(xiàn)與BA 顯著相關(guān)的基因主要富集在細(xì)胞外基質(zhì)（extracellular matrix，ECM）、細(xì)胞間黏附分子、色氨酸代謝、甘油磷脂代謝、PPAR 信號(hào)通路。ECM 含有大量信號(hào)分子，參與調(diào)節(jié)細(xì)胞的遷移、增殖、凋亡或分化，細(xì)胞不斷通過合成、降解、重新組裝和化學(xué)改造來重建和改造ECM，以保持組織平衡，不良的ECM 可能加劇疾病進(jìn)展[13]。PPAR 受體激動(dòng)劑能減少腎小球硬化和腎小管間質(zhì)損傷[14]。

本研究根據(jù)模塊身份和基因重要性識(shí)別出16 個(gè)BA 相關(guān)的樞紐基因，包括MMP7、SPP1、VCAN、VTCN1、LUM、RGS4、EpCAM、ELOVL7、MACC1、LAMC2、C12orf75、COL4A1、VIM、TMEM156、TMSB10和MYH9。有研究顯示MMP7 是一種肝星形細(xì)胞分泌的蛋白酶，與BA 的肝纖維化與組織重塑相關(guān)[15]。血清蛋白組學(xué)研究提示MMP7 測(cè)定具有高靈敏度和特異性，可區(qū)別BA 與非BA 的膽汁淤積性疾病，可作為BA的可靠的生物標(biāo)志[16-18]。這些提示MMP7 直接參與了BA 發(fā)病，但MMP7 對(duì)于診斷BA 的截?cái)嘀瞪形唇y(tǒng)一（1.43~52.85 ng/mL）[18]，因此MMP7 的臨床應(yīng)用仍需要更大規(guī)模的研究。SPP1 在膽汁異常患者中被過度表達(dá)，但它的滅活并沒有影響表型的發(fā)展[19]，這表明SPP1 可能涉及肝臟對(duì)膽汁損傷的一般生理反應(yīng)，但不一定與BA 的發(fā)病機(jī)制有關(guān)。在BA 患者及BA 動(dòng)物模型的肝臟中均檢測(cè)到未成熟的EpCAM，且EpCAM 陽(yáng)性的不成熟膽道上皮細(xì)胞在短時(shí)間內(nèi)迅速增生，EpCAM 可能是BA 患兒肝纖維化進(jìn)展迅速的主要原因[20]。文獻(xiàn)報(bào)道ELOVL7 能產(chǎn)生對(duì)于構(gòu)建正常功能的人巨細(xì)胞病毒包膜必需的飽和超長(zhǎng)鏈脂肪酸[21]，這與BA 的病因可能是病毒感染的說法一致。已有研究顯示16 個(gè)樞紐基因中的部分基因在BA 發(fā)病中可能發(fā)揮重要作用，表明本研究采用的鑒定樞紐基因的方法精確可靠。VCAN、VTCN1、LUM、ELOVL7、MACC1、LAMC2、C12orf75、COL4A1、VIM、TMEM156、TMSB10 及MYH9 在BA發(fā)病機(jī)制中的研究，樞紐基因及富集分析所示功能是否能減緩甚至逆轉(zhuǎn)BA 的進(jìn)展，有待進(jìn)一步深入研究。

綜上所述，基于WGCNA 識(shí)別了與BA 發(fā)病相關(guān)的1 個(gè)關(guān)鍵模塊及16 個(gè)樞紐基因，新發(fā)現(xiàn)了12 個(gè)可能和BA 的發(fā)病密切相關(guān)的基因，值得進(jìn)一步實(shí)驗(yàn)驗(yàn)證。