轉錄因子Bach1在小鼠胚胎發育時期的生物信息學分析

呂曉宇 王新紅 孟 丹

（復旦大學基礎醫學院生理與病理生理學系 上海 200032）

胚胎發育作為一種基本的生物學過程，具有復雜的時空調控網絡。利用分子生物技術和生物信息學方法，可以揭示多種調節因子協同作用的調控網絡［1-5］。加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）是利用基因表達數據構建無標度網絡的系統生物學方法。它能尋找協同基因模塊，探索基因網絡和感興趣的表型之間的關系，以及網絡中的中樞基因。WGCNA 分析方法目前已成為胚胎發育研究領域的一種重要生物信息學分析手段［6-8］。我們在前期研究中發現，Bach1 維持人胚胎干細胞的干細胞特性，通過募集去泛素化酶Usp7 來穩定多能性因子，通過招募PRC2 復合體沉默中胚層和內胚層基因表達，并抑制Wnt3 和Nodal信號通路［9］。全身敲除Bach1的小鼠有亞致死的表型，這提示Bach1是胚胎發育過程中的一個重要轉錄因子。但是，Bach1在小鼠胚胎發育中的作用意義并不清楚。

本研究利用小鼠胚胎發育不同時期的轉錄組數據［獲自Gene Expression Omnibus（GEO）數據庫（GSE92634）］，用生物信息學分析在小鼠胚胎發育過程中與Bach1共表達的基因網絡，分析Bach1共表達的基因生物功能，以期了解Bach1在胚胎發育中可能的調控作用。為深入研究Bach1在胚胎發育中的作用提供有意義的生物信息學依據。

材料和方法

小鼠胚胎樣本RNA-seq 數據的表達分析GEO數據庫由美國國立生物技術信息中心NCBI 建立，數據來源于世界各個科研機構及組織提交的數據，具有較好的生物信息數據資源。作為一個服務于廣大科研工作者的免費數據庫，GEO 資源有大量質量較高的高通量測序數據，包括RNA-seq 數據、ChIP-seq 數據等，數據結果較為可信。GEO 數據庫包含已發表文獻的數據，原始數據包括GPL（Platform）、GSM（Sample）和 GSE（Series）。GEO數據庫整理后的數據包括數據集GDS（DataSets）和表達譜（Profiles）。GEO 根據平臺、數據集、系列和樣本等4 種形式組織數據。

我們從GEO 數據庫下載小鼠胚胎的RNA 序列數據（GSE92634），并利用另外一個小鼠胚胎的轉錄組數據庫（GSE120963）作為驗證。原始數據集已對FPKM 數值進行了log2 校正。我們采用校正后的數據進行mRNA 表達分析，使用WGCNA 算法評估基因表達值，使用R 語言的flashClust 工具對具有適當閾值的樣本進行聚類分析。

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

小鼠胚胎共表達模塊構建分析用WGCNA算法在模塊構建中篩選出功率值。梯度法用于測試不同功率值（范圍：1～20）的模塊的獨立性和平均連通度。當獨立程度為0.9 時確定適當的功率值，用WGCNA 算法繼續構造模塊，提取出每個模塊的相應基因信息，將最小基因數設定為30。應用WGCNA 算法識別共表達模塊，在R 軟件包（http：//www.r-project.org/）中實現，并繪制熱圖。

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

共表達模塊的功能富集分析構建的模塊的數量由基因的數量決定，然后對這些模塊中的基因進行功能富集分析。將模塊的基因信息輸入到Metascape（用于注釋、可視化和集成發現的數據庫）數據集（http：//metascape.org/）進行功能富集［10］。提取GO（Gene Ontology）分析結果，對感興趣的模塊用Cytoscape 軟件進行可視化分析。

統計分析數據用±s表示，多組間比較用單因素方差分析（One-way ANOVA），P＜0.05 為差異有統計學意義。用GraphPad Prism 軟件進行分析。

結 果

小鼠胚胎發育表達矩陣的WGCNA 模型構建從表達譜數據和表型數據矩陣中獲得總共10 個樣品和17 360 個基因。 根據表達譜計算每個樣品中每個基因的方差，選擇標準偏差大于1.2 的基因并進一步聚類所有樣品（圖1A）。所有樣品均符合WGCNA 的選擇標準，未排除樣品。具有樣本特征的聚類結果如圖1B 所示。圖中的紅色表示表型表中標記為非零的樣品。我們得到一個新的數據表達譜，其中包含10 個樣本和8 933 個基因（圖1B）。表達矩陣被轉換為鄰接矩陣，鄰接矩陣被轉換為拓撲矩陣。基于TOM 矩陣，使用平均連鎖層次聚類方法來聚類基因。根據混合動態切割樹的標準，將每個基因網絡模塊的最小基數設置為30。通過動態剪切方法確定基因模塊后，依次計算每個模塊的特征向量，然后聚類模塊并將更近的模塊合并到新模塊中。總共獲得17 個模塊，其中灰色模塊是不能聚合到其他模塊中的基因的集合。Bach1被富集到黃綠模塊中（圖1C）。

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

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

圖2 經GO 富集分析提取的Bach1 共表達基因相關模塊Fig 2 Extraction of the Bach1 co-expressing genes related module by GO enrichment analysis

WGCNA 富集Bach1 模塊的共表達網絡構建根據每個模塊中基因的共表達權重，權重的閾值為0.02。提取每個模塊的共表達網絡文件，并將其導入Cytoscape 可視化。將Bach1 所在模塊構建PPI網絡后發現，與Bach1在小鼠胚胎中共表達的基因包括VEGFb、Fam105b等與血管發育密切相關的基因，以及Tcf15、Znf622等與胚胎發育密切相關的基因（圖3）。

圖3 WGCNA 富集Bach1 模塊的PPI 網絡構建Fig 3 Construction of PPI Network for WGCNA enrichment Bach1 module

小鼠胚胎發育過程中Bach1 與Vegfb 基因的表達變化通過分析小鼠胚胎發育E5.5Epi、E6.0Epi、E6.5P、E7.0P、E7.5P 等不同時期的數據，我們發現Bach1與Vegfb基因的表達趨勢非常相似，都是在E5.5Epi 到 E6.0Epi 期 表 達 增 加 ，E6.5P 到 E7.5P 期表達降低（圖4），二者具有很好的共表達相關性。為了驗證這一結果，我們用小鼠胚胎發育的另一個轉錄組數據庫（GSE120963）分析了胚胎發育不同時期Bach1與Vegfb基因的表達，發現二者從E5.5A到E7.0A 期表達趨勢也基本相似。與E5.5A 期相比，Bach1的表達在E6.0A 期增加，差異有統計學意義（P＜0.05，圖5）。

Bach1 富集在一些共表達基因啟動子或增強子區使用UCSC 數據庫中人胚胎干細胞Bach1的ChIP-seq 數據（GSE31477）［11］，分析與Bach1基因具有共表達關系的基因是否受Bach1的調控。結果發現，Bach1 蛋白在一些與其共表達的基因，如Tcf15、Znf622和Fam105b基因啟動子或增強子區具有較高的富集（圖6），已知這些基因與胚胎發育、血管生成有著密切的關系［12-14］。這些結果表明Bach1 可能直接調控這些基因的轉錄，影響胚胎發育和血管生成。

圖4 轉錄組數據庫（GSE92634）分析Bach1 和Vegfb 基因在小鼠胚胎發育E5.5Epi 到E7.5P 的表達Fig 4 mRNA expression of Bach1 and Vegfb from E5.5Epi to E7.5P in mouse embryos from transcription group data（GSE92634）

圖5 轉錄組數據庫（GSE120963）分析Bach1 和Vegfb 基因在小鼠胚胎發育E5.5A 到E7.0A 的表達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 富集分析對WGCNA 建立的共表達模型中Bach1 聚集的基因模塊進行GO 通路富集分析，結果發現多個與Wnt信號通路、蛋白修飾與翻譯、染色質重塑、DNA 損傷反應、細胞周期調節、泛素化修飾等功能密切相關的基因模塊。我們以往研究已證實Bach1 通過乙酰化調控Wnt信號通路，Bach1 通過招募去泛素化酶增加多能性基因蛋白穩定性，這些都與我們預測的Bach1 聚集的基因功能具有高度的一致性，提示該生物信息學分析具有可信性（圖7）。

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

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

討 論

胚胎發育不同時間和空間的分子調控機制研究對再生醫學具有重要意義。本研究分析了胚胎E5.5 天到E7.5 天的小鼠胚胎發育過程中不同天數和不同部位的RNA-seq 數據，通過對數據的提取和標準化，進行WGCNA 分析，旨在了解Bach1 在胚胎發育過程中的時空調控規律。

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

綜上所述，我們通過生物信息學分析，發現Bach1 在小鼠胚胎發育中可能起著重要作用，其與血管發育密切相關。這些信息對闡明胚胎發育過程中的調控網絡具有重要的參考價值。