呂曉宇 王新紅 孟 丹

（復(fù)旦大學(xué)基礎(chǔ)醫(yī)學(xué)院生理與病理生理學(xué)系 上海 200032）

胚胎發(fā)育作為一種基本的生物學(xué)過程，具有復(fù)雜的時(shí)空調(diào)控網(wǎng)絡(luò)。利用分子生物技術(shù)和生物信息學(xué)方法，可以揭示多種調(diào)節(jié)因子協(xié)同作用的調(diào)控網(wǎng)絡(luò)［1-5］。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）是利用基因表達(dá)數(shù)據(jù)構(gòu)建無(wú)標(biāo)度網(wǎng)絡(luò)的系統(tǒng)生物學(xué)方法。它能尋找協(xié)同基因模塊，探索基因網(wǎng)絡(luò)和感興趣的表型之間的關(guān)系，以及網(wǎng)絡(luò)中的中樞基因。WGCNA 分析方法目前已成為胚胎發(fā)育研究領(lǐng)域的一種重要生物信息學(xué)分析手段［6-8］。我們?cè)谇捌谘芯恐邪l(fā)現(xiàn)，Bach1 維持人胚胎干細(xì)胞的干細(xì)胞特性，通過募集去泛素化酶Usp7 來(lái)穩(wěn)定多能性因子，通過招募PRC2 復(fù)合體沉默中胚層和內(nèi)胚層基因表達(dá)，并抑制Wnt3 和Nodal信號(hào)通路［9］。全身敲除Bach1的小鼠有亞致死的表型，這提示Bach1是胚胎發(fā)育過程中的一個(gè)重要轉(zhuǎn)錄因子。但是，Bach1在小鼠胚胎發(fā)育中的作用意義并不清楚。

本研究利用小鼠胚胎發(fā)育不同時(shí)期的轉(zhuǎn)錄組數(shù)據(jù)［獲自Gene Expression Omnibus（GEO）數(shù)據(jù)庫(kù)（GSE92634）］，用生物信息學(xué)分析在小鼠胚胎發(fā)育過程中與Bach1共表達(dá)的基因網(wǎng)絡(luò)，分析Bach1共表達(dá)的基因生物功能，以期了解Bach1在胚胎發(fā)育中可能的調(diào)控作用。為深入研究Bach1在胚胎發(fā)育中的作用提供有意義的生物信息學(xué)依據(jù)。

材料和方法

小鼠胚胎樣本RNA-seq 數(shù)據(jù)的表達(dá)分析GEO數(shù)據(jù)庫(kù)由美國(guó)國(guó)立生物技術(shù)信息中心NCBI 建立，數(shù)據(jù)來(lái)源于世界各個(gè)科研機(jī)構(gòu)及組織提交的數(shù)據(jù)，具有較好的生物信息數(shù)據(jù)資源。作為一個(gè)服務(wù)于廣大科研工作者的免費(fèi)數(shù)據(jù)庫(kù)，GEO 資源有大量質(zhì)量較高的高通量測(cè)序數(shù)據(jù)，包括RNA-seq 數(shù)據(jù)、ChIP-seq 數(shù)據(jù)等，數(shù)據(jù)結(jié)果較為可信。GEO 數(shù)據(jù)庫(kù)包含已發(fā)表文獻(xiàn)的數(shù)據(jù)，原始數(shù)據(jù)包括GPL（Platform）、GSM（Sample）和 GSE（Series）。GEO數(shù)據(jù)庫(kù)整理后的數(shù)據(jù)包括數(shù)據(jù)集GDS（DataSets）和表達(dá)譜（Profiles）。GEO 根據(jù)平臺(tái)、數(shù)據(jù)集、系列和樣本等4 種形式組織數(shù)據(jù)。

我們從GEO 數(shù)據(jù)庫(kù)下載小鼠胚胎的RNA 序列數(shù)據(jù)（GSE92634），并利用另外一個(gè)小鼠胚胎的轉(zhuǎn)錄組數(shù)據(jù)庫(kù)（GSE120963）作為驗(yàn)證。原始數(shù)據(jù)集已對(duì)FPKM 數(shù)值進(jìn)行了log2 校正。我們采用校正后的數(shù)據(jù)進(jìn)行mRNA 表達(dá)分析，使用WGCNA 算法評(píng)估基因表達(dá)值，使用R 語(yǔ)言的flashClust 工具對(duì)具有適當(dāng)閾值的樣本進(jìn)行聚類分析。

我們使用加利福尼亞大學(xué)圣克魯茲分校（University of California Santa Cruz，UCSC）數(shù)據(jù)庫(kù)中人胚胎干細(xì)胞Bach1 的染色質(zhì)免疫沉淀和測(cè)序（chromatin immunoprecipitation and sequencing，ChIP-seq）數(shù)據(jù)庫(kù)（GSE31477），使用 UCSC 基因組瀏覽器進(jìn)行可視化。MACS2 用于鑒定Bach1的富集峰（默認(rèn)參數(shù)為“銳”峰）。

小鼠胚胎共表達(dá)模塊構(gòu)建分析用WGCNA算法在模塊構(gòu)建中篩選出功率值。梯度法用于測(cè)試不同功率值（范圍：1～20）的模塊的獨(dú)立性和平均連通度。當(dāng)獨(dú)立程度為0.9 時(shí)確定適當(dāng)?shù)墓β手?，用WGCNA 算法繼續(xù)構(gòu)造模塊，提取出每個(gè)模塊的相應(yīng)基因信息，將最小基因數(shù)設(shè)定為30。應(yīng)用WGCNA 算法識(shí)別共表達(dá)模塊，在R 軟件包（http：//www.r-project.org/）中實(shí)現(xiàn)，并繪制熱圖。

構(gòu)建小鼠胚胎的模塊-特征關(guān)系使用模塊eigengene 和表型之間的相關(guān)性來(lái)估計(jì)模塊-性狀關(guān)聯(lián)，進(jìn)一步鑒定與表型高度相關(guān)的模塊。對(duì)于每種表達(dá)譜，用基因顯著性（gene significane，GS）計(jì)算表達(dá)譜與每種性狀之間相關(guān)性的絕對(duì)值；絕對(duì)值＞0的基因被聚類到模塊中進(jìn)行特征基因與表型的相關(guān)性分析。

共表達(dá)模塊的功能富集分析構(gòu)建的模塊的數(shù)量由基因的數(shù)量決定，然后對(duì)這些模塊中的基因進(jìn)行功能富集分析。將模塊的基因信息輸入到Metascape（用于注釋、可視化和集成發(fā)現(xiàn)的數(shù)據(jù)庫(kù)）數(shù)據(jù)集（http：//metascape.org/）進(jìn)行功能富集［10］。提取GO（Gene Ontology）分析結(jié)果，對(duì)感興趣的模塊用Cytoscape 軟件進(jìn)行可視化分析。

統(tǒng)計(jì)分析數(shù)據(jù)用±s表示，多組間比較用單因素方差分析（One-way ANOVA），P＜0.05 為差異有統(tǒng)計(jì)學(xué)意義。用GraphPad Prism 軟件進(jìn)行分析。

結(jié) 果

小鼠胚胎發(fā)育表達(dá)矩陣的WGCNA 模型構(gòu)建從表達(dá)譜數(shù)據(jù)和表型數(shù)據(jù)矩陣中獲得總共10 個(gè)樣品和17 360 個(gè)基因。 根據(jù)表達(dá)譜計(jì)算每個(gè)樣品中每個(gè)基因的方差，選擇標(biāo)準(zhǔn)偏差大于1.2 的基因并進(jìn)一步聚類所有樣品（圖1A）。所有樣品均符合WGCNA 的選擇標(biāo)準(zhǔn)，未排除樣品。具有樣本特征的聚類結(jié)果如圖1B 所示。圖中的紅色表示表型表中標(biāo)記為非零的樣品。我們得到一個(gè)新的數(shù)據(jù)表達(dá)譜，其中包含10 個(gè)樣本和8 933 個(gè)基因（圖1B）。表達(dá)矩陣被轉(zhuǎn)換為鄰接矩陣，鄰接矩陣被轉(zhuǎn)換為拓?fù)渚仃??；赥OM 矩陣，使用平均連鎖層次聚類方法來(lái)聚類基因。根據(jù)混合動(dòng)態(tài)切割樹的標(biāo)準(zhǔn)，將每個(gè)基因網(wǎng)絡(luò)模塊的最小基數(shù)設(shè)置為30。通過動(dòng)態(tài)剪切方法確定基因模塊后，依次計(jì)算每個(gè)模塊的特征向量，然后聚類模塊并將更近的模塊合并到新模塊中?？偣搏@得17 個(gè)模塊，其中灰色模塊是不能聚合到其他模塊中的基因的集合。Bach1被富集到黃綠模塊中（圖1C）。

小鼠胚胎轉(zhuǎn)錄組的WGCNA 聚集模塊的Bach1 GO 富集分析根據(jù)每個(gè)模塊的特征向量計(jì)算這些模塊和每個(gè)表型之間的相關(guān)性（圖2）。對(duì)富集到的模塊進(jìn)行GS 分布預(yù)測(cè)。GS=0 表示該基因與表型無(wú)關(guān)。每個(gè)模塊中每種表型的GS 分布可以顯示表型和基因之間的整體相關(guān)性（圖2A）。根據(jù)每個(gè)模塊的特征向量，計(jì)算這些模塊之間的相關(guān)性。模塊聚集在同一個(gè)分支中，結(jié)果發(fā)現(xiàn)黃綠模塊和magenta 聚集到一起（圖2B）。對(duì)黃綠模塊關(guān)鍵基因Bach1進(jìn)行檢驗(yàn) ，結(jié)果 R 值為 0.820，P值為 0.003。所以選擇Bach1為樞紐基因。

圖1 小鼠胚胎（E5.5-E7.5）基因表達(dá)矩陣的WGCNA 模型構(gòu)建Fig 1 Construction of WGCNA model of mouse embryo（E5.5-E7.5）developmental expression matrix

圖2 經(jīng)GO 富集分析提取的Bach1 共表達(dá)基因相關(guān)模塊Fig 2 Extraction of the Bach1 co-expressing genes related module by GO enrichment analysis

WGCNA 富集Bach1 模塊的共表達(dá)網(wǎng)絡(luò)構(gòu)建根據(jù)每個(gè)模塊中基因的共表達(dá)權(quán)重，權(quán)重的閾值為0.02。提取每個(gè)模塊的共表達(dá)網(wǎng)絡(luò)文件，并將其導(dǎo)入Cytoscape 可視化。將Bach1 所在模塊構(gòu)建PPI網(wǎng)絡(luò)后發(fā)現(xiàn)，與Bach1在小鼠胚胎中共表達(dá)的基因包括VEGFb、Fam105b等與血管發(fā)育密切相關(guān)的基因，以及Tcf15、Znf622等與胚胎發(fā)育密切相關(guān)的基因（圖3）。

圖3 WGCNA 富集Bach1 模塊的PPI 網(wǎng)絡(luò)構(gòu)建Fig 3 Construction of PPI Network for WGCNA enrichment Bach1 module

小鼠胚胎發(fā)育過程中Bach1 與Vegfb 基因的表達(dá)變化通過分析小鼠胚胎發(fā)育E5.5Epi、E6.0Epi、E6.5P、E7.0P、E7.5P 等不同時(shí)期的數(shù)據(jù)，我們發(fā)現(xiàn)Bach1與Vegfb基因的表達(dá)趨勢(shì)非常相似，都是在E5.5Epi 到 E6.0Epi 期 表 達(dá) 增 加 ，E6.5P 到 E7.5P 期表達(dá)降低（圖4），二者具有很好的共表達(dá)相關(guān)性。為了驗(yàn)證這一結(jié)果，我們用小鼠胚胎發(fā)育的另一個(gè)轉(zhuǎn)錄組數(shù)據(jù)庫(kù)（GSE120963）分析了胚胎發(fā)育不同時(shí)期Bach1與Vegfb基因的表達(dá)，發(fā)現(xiàn)二者從E5.5A到E7.0A 期表達(dá)趨勢(shì)也基本相似。與E5.5A 期相比，Bach1的表達(dá)在E6.0A 期增加，差異有統(tǒng)計(jì)學(xué)意義（P＜0.05，圖5）。

Bach1 富集在一些共表達(dá)基因啟動(dòng)子或增強(qiáng)子區(qū)使用UCSC 數(shù)據(jù)庫(kù)中人胚胎干細(xì)胞Bach1的ChIP-seq 數(shù)據(jù)（GSE31477）［11］，分析與Bach1基因具有共表達(dá)關(guān)系的基因是否受Bach1的調(diào)控。結(jié)果發(fā)現(xiàn)，Bach1 蛋白在一些與其共表達(dá)的基因，如Tcf15、Znf622和Fam105b基因啟動(dòng)子或增強(qiáng)子區(qū)具有較高的富集（圖6），已知這些基因與胚胎發(fā)育、血管生成有著密切的關(guān)系［12-14］。這些結(jié)果表明Bach1 可能直接調(diào)控這些基因的轉(zhuǎn)錄，影響胚胎發(fā)育和血管生成。

圖4 轉(zhuǎn)錄組數(shù)據(jù)庫(kù)（GSE92634）分析Bach1 和Vegfb 基因在小鼠胚胎發(fā)育E5.5Epi 到E7.5P 的表達(dá)Fig 4 mRNA expression of Bach1 and Vegfb from E5.5Epi to E7.5P in mouse embryos from transcription group data（GSE92634）

圖5 轉(zhuǎn)錄組數(shù)據(jù)庫(kù)（GSE120963）分析Bach1 和Vegfb 基因在小鼠胚胎發(fā)育E5.5A 到E7.0A 的表達(dá)Fig 5 mRNA expression of Bach1 and Vegfb from E5.5A to E7.0A of mouse embryos from transcription group data（GSE120963）

WGCNA 聚集Bach1 模塊的GO 富集分析對(duì)WGCNA 建立的共表達(dá)模型中Bach1 聚集的基因模塊進(jìn)行GO 通路富集分析，結(jié)果發(fā)現(xiàn)多個(gè)與Wnt信號(hào)通路、蛋白修飾與翻譯、染色質(zhì)重塑、DNA 損傷反應(yīng)、細(xì)胞周期調(diào)節(jié)、泛素化修飾等功能密切相關(guān)的基因模塊。我們以往研究已證實(shí)Bach1 通過乙?；{(diào)控Wnt信號(hào)通路，Bach1 通過招募去泛素化酶增加多能性基因蛋白穩(wěn)定性，這些都與我們預(yù)測(cè)的Bach1 聚集的基因功能具有高度的一致性，提示該生物信息學(xué)分析具有可信性（圖7）。

圖6 胚胎干細(xì)胞中Bach1 在共表達(dá)基因Tcf15、Znf622 和Fam105b 上的信號(hào)富集Fig 6 Bach1 signal tracks for representative loci Tcf15，Znf622，and Fam105b in embryonic stem cells

圖7 小鼠胚胎轉(zhuǎn)錄組的WGCNA 聚集模塊的Bach1 GO 富集分析Fig 7 Bach1 GO enrichment analysis of WGCNA aggregation module of mouse embryonic transcriptome

討 論

胚胎發(fā)育不同時(shí)間和空間的分子調(diào)控機(jī)制研究對(duì)再生醫(yī)學(xué)具有重要意義。本研究分析了胚胎E5.5 天到E7.5 天的小鼠胚胎發(fā)育過程中不同天數(shù)和不同部位的RNA-seq 數(shù)據(jù)，通過對(duì)數(shù)據(jù)的提取和標(biāo)準(zhǔn)化，進(jìn)行WGCNA 分析，旨在了解Bach1 在胚胎發(fā)育過程中的時(shí)空調(diào)控規(guī)律。

WGCNA 分析側(cè)重于共表達(dá)模塊和表型特征之間的關(guān)聯(lián)，因此與其他方法相比，分析結(jié)果具有更高的可靠性和生物學(xué)意義。在本研究中，我們通過WGCNA 方法構(gòu)建共表達(dá)模塊，分析了Bach1 的胚胎發(fā)育E5.5 天到E7.5 天的表達(dá)數(shù)據(jù)，對(duì)來(lái)自10個(gè)小鼠胚胎樣品的8 933 個(gè)基因構(gòu)建了總共17 個(gè)共表達(dá)模塊，應(yīng)用于研究模塊和表型的基礎(chǔ)關(guān)系。通過分析，確定了兩個(gè)與小鼠胚胎發(fā)育顯著相關(guān)的共表達(dá)模塊，用于檢測(cè)小鼠胚胎轉(zhuǎn)錄組與小鼠胚胎發(fā)生的時(shí)空轉(zhuǎn)錄調(diào)節(jié)之間的關(guān)系。進(jìn)一步對(duì)特定模塊的這些共表達(dá)基因進(jìn)行功能富集分析，發(fā)現(xiàn)轉(zhuǎn)錄因子Bach1被WGCNA 共表達(dá)網(wǎng)絡(luò)鑒定為胚胎發(fā)育重要的時(shí)空監(jiān)管中樞基因。我們發(fā)現(xiàn)Bach1與Vegfb基因的表達(dá)趨勢(shì)基本一致，二者具有很好的相關(guān)性。由于分析的小鼠胚胎發(fā)育不同時(shí)期RNA的數(shù)據(jù)集樣本數(shù)量偏少，為了驗(yàn)證這一結(jié)果的可信性，我們對(duì)小鼠早期胚胎發(fā)育的另一個(gè)轉(zhuǎn)錄組數(shù)據(jù)集進(jìn)行了分析，發(fā)現(xiàn)Bach1與Vegfb基因的表達(dá)趨勢(shì)也非常相似。已知Vegfb在斑馬魚的血管發(fā)育中起重要作用，因此我們推測(cè)Bach1 與血管發(fā)育密切相關(guān)［15］。兩個(gè)數(shù)據(jù)集分析發(fā)現(xiàn) Bach1 從 E5.5 天到E6.0 天表達(dá)均增加，隨后E6.5 天開始降低。我們最近的研究證實(shí)，Bach1 抑制人胚胎干細(xì)胞向中內(nèi)胚層細(xì)胞的分化［9］，因此推測(cè)可能在胚胎發(fā)育早期Bach1 短暫表達(dá)增高抑制中內(nèi)胚層分化，但隨著胚胎發(fā)育的進(jìn)行，Bach1 表達(dá)降低，這種抑制作用減弱。此外，我們還發(fā)現(xiàn)Bach1 蛋白存在一些與其共表達(dá)的基因，如Tcf15、Znf622和Fam105b基因啟動(dòng)子區(qū)具有較高的富集。已有文獻(xiàn)報(bào)道Tcf15、Znf622和Fam105b與胚胎發(fā)育和血管生成有著密切的關(guān)系，提示在小鼠胚胎發(fā)育過程中，Bach1 可能直接調(diào)控與胚胎發(fā)育、血管生成相關(guān)的基因表達(dá)。分析與Bach1 共表達(dá)基因的 GO 富集發(fā)現(xiàn)，Bach1 與Wnt信號(hào)通路、蛋白修飾與翻譯、染色質(zhì)重塑、DNA 損傷反應(yīng)、細(xì)胞周期調(diào)節(jié)、泛素化修飾等密切相關(guān)［16-17］。全身敲除Bach1 的小鼠胚胎有亞致死的表型，但Bach1 在胚胎發(fā)育中的作用仍不清楚。我們的研究表明Bach1通過促進(jìn)干細(xì)胞維持自我更新，抑制中內(nèi)胚層分化。分析小鼠胚胎發(fā)育不同時(shí)期轉(zhuǎn)錄組數(shù)據(jù)后，我們發(fā)現(xiàn)Bach1在小鼠胚胎E5.5 天到E7.5天，與Wnt信號(hào)通路、染色質(zhì)重塑、DNA 損傷反應(yīng)、細(xì)胞周期調(diào)節(jié)、泛素化修飾等密切相關(guān)，這與我們之前的研究相符合，提示本次生物信息學(xué)分析結(jié)果具有可信度。后續(xù)我們將在Bach1內(nèi)皮細(xì)胞特異敲除小鼠上進(jìn)一步驗(yàn)證本次分析結(jié)果，研究Bach1 對(duì)小鼠胚胎血管發(fā)育的影響。

綜上所述，我們通過生物信息學(xué)分析，發(fā)現(xiàn)Bach1 在小鼠胚胎發(fā)育中可能起著重要作用，其與血管發(fā)育密切相關(guān)。這些信息對(duì)闡明胚胎發(fā)育過程中的調(diào)控網(wǎng)絡(luò)具有重要的參考價(jià)值。