转录因子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 在小鼠胚胎发育中可能起着重要作用，其与血管发育密切相关。这些信息对阐明胚胎发育过程中的调控网络具有重要的参考价值。