2020, Vol. 47
胚胎发育作为一种基本的生物学过程，具有复杂的时空调控网络。利用分子生物技术和生物信息学方法，可以揭示多种调节因子协同作用的调控网络[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软件进行可视化分析。
统计分析 数据用x±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）。
|
| A:Correlation clustering of expression matrix gene expression data in mouse embryos of E5.5-E7.5;B:A cluster heat map was drawn for mouse embryo samples of E5.5-E7.5;C:The model of WGCNA. 图 1 小鼠胚胎（E5.5-E7.5）基因表达矩阵的WGCNA模型构建 Fig 1 Construction of WGCNA model of mouse embryo (E5.5-E7.5) developmental expression matrix |
小鼠胚胎转录组的WGCNA聚集模块的Bach1 GO富集分析 根据每个模块的特征向量计算这些模块和每个表型之间的相关性(图 2)。对富集到的模块进行GS分布预测。GS=0表示该基因与表型无关。每个模块中每种表型的GS分布可以显示表型和基因之间的整体相关性（图 2A）。根据每个模块的特征向量，计算这些模块之间的相关性。模块聚集在同一个分支中，结果发现黄绿模块和magenta聚集到一起（图 2B）。对黄绿模块关键基因Bach1进行检验，结果R值为0.820，P值为0.003。所以选择Bach1为枢纽基因。
|
| A:The construction of GS distribution of module gene and phenotype; B:The correlation clustering of WGCNA module. 图 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）。
|
| mRNA level was calculated using the log 2 FPKM.Epi:Eepiblast; P:Posterior. 图 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) |
|
| mRNA level was calculated using the log 2 FPKM.A:Anterior.(1)P < 0.05 compared with E5.5A;One-way analysis of variance. 图 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) |
Bach1富集在一些共表达基因启动子或增强子区 使用UCSC数据库中人胚胎干细胞Bach1的ChIP-seq数据（GSE31477）[11]，分析与Bach1基因具有共表达关系的基因是否受Bach1的调控。结果发现，Bach1蛋白在一些与其共表达的基因，如Tcf15、Znf622和Fam105b基因启动子或增强子区具有较高的富集（图 6），已知这些基因与胚胎发育、血管生成有着密切的关系[12-14]。这些结果表明Bach1可能直接调控这些基因的转录，影响胚胎发育和血管生成。
|
| 图 6 胚胎干细胞中Bach1在共表达基因Tcf15、Znf622和Fam105b上的信号富集 Fig 6 Bach1 signal tracks for representative loci Tcf15, Znf622, and Fam105b in embryonic stem cells |
WGCNA聚集Bach1模块的GO富集分析 对WGCNA建立的共表达模型中Bach1聚集的基因模块进行GO通路富集分析，结果发现多个与Wnt信号通路、蛋白修饰与翻译、染色质重塑、DNA损伤反应、细胞周期调节、泛素化修饰等功能密切相关的基因模块。我们以往研究已证实Bach1通过乙酰化调控Wnt信号通路，Bach1通过招募去泛素化酶增加多能性基因蛋白稳定性，这些都与我们预测的Bach1聚集的基因功能具有高度的一致性，提示该生物信息学分析具有可信性（图 7）。
|
| 图 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在小鼠胚胎发育中可能起着重要作用，其与血管发育密切相关。这些信息对阐明胚胎发育过程中的调控网络具有重要的参考价值。
| [1] |
MOREY L, SANTANACH A, BLANCO E, et al. Polycomb regulates mesoderm cell fate-specification in embryonic stem cells through activation and repression mechanisms[J]. Cell Stem Cell, 2015, 17(3): 300-315.
[DOI]
|
| [2] |
BOWER NI, VOGRIN AJ, LE GUEN L, et al. Vegfd modulates both angiogenesis and lymphangiogenesis during zebrafish embryonic development[J]. Development, 2017, 144(3): 507-518.
[DOI]
|
| [3] |
YUAN S, SCHUSTER A, TANG C, et al. Sperm-borne miRNAs and endo-siRNAs are important for fertilization and preimplantation embryonic development[J]. Development, 2016, 143(4): 635-647.
[DOI]
|
| [4] |
GRIJZENHOUT A, GODWIN J, KOSEKI H, et al. Functional analysis of AEBP2, a PRC2 Polycomb protein, reveals a Trithorax phenotype in embryonic development and in ESCs[J]. Development, 2016, 143(15): 2716-2723.
[DOI]
|
| [5] |
LIU C, WANG R, HE Z, et al. Suppressing nodal signaling activity predisposes ectodermal differentiation of epiblast stem cells[J]. Stem Cell Reports, 2018, 11(1): 43-57.
[URI]
|
| [6] |
WAN Q, TANG J, HAN Y, et al. Co-expression modules construction by WGCNA and identify potential prognostic markers of uveal melanoma[J]. Exp Eye Res, 2018, 166: 13-20.
[DOI]
|
| [7] |
DING M, LI F, WANG B, et al. A comprehensive analysis of WGCNA and serum metabolomics manifests the lung cancer-associated disordered glucose metabolism[J]. J Cell Biochem, 2018, 120(6): 10855-10863.
[URI]
|
| [8] |
LANGFELDER P, HORVATH S. WGCNA:an R package for weighted correlation network analysis[J]. BMC Bioinformatics, 2008, 9: 559.
[DOI]
|
| [9] |
WEI X, GUO J, LI Q, et al. Bach1 regulates self-renewal and impedes mesendodermal differentiation of human embryonic stem cells[J]. Sci Adv, 2019, 5(3): u7887.
[DOI]
|
| [10] |
ZHOU Y, ZHOU B, PACHE L, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets[J]. Nat Commun, 2019, 10(1): 1523.
[DOI]
|
| [11] |
POPE BD, RYBA T, DILEEP V, et al. Topologically associating domains are stable units of replication-timing regulation[J]. Nature, 2014, 515(7527): 402-405.
[DOI]
|
| [12] |
AKAGI T, KUURE S, URANISHI K, et al. ETS-related transcription factors ETV4 and ETV5 are involved in proliferation and induction of differentiation-associated genes in embryonic stem (ES) cells[J]. J Biol Chem, 2015, 290(37): 22460-22473.
[DOI]
|
| [13] |
RIVKIN E, ALMEIDA SM, CECCARELLI DF, et al. The linear ubiquitin-specific deubiquitinase gumby regulates angiogenesis[J]. Nature, 2013, 498(7454): 318-324.
[DOI]
|
| [14] |
SEONG H, JUNG H, MANOHARAN R, et al. Positive regulation of apoptosis signal-regulating kinase 1 signaling by ZPR9 protein, a zinc finger protein[J]. J Biol Chem, 2011, 286(36): 31123-31135.
[DOI]
|
| [15] |
JENSEN LD, NAKAMURA M, BRÄUTIGAM L, et al. VEGF-B-Neuropilin-1 signaling is spatiotemporally indispensable for vascular and neuronal development in zebrafish[J]. Proc Natl Acad Sci U S A, 2015, 112(44): E5944-E5953.
[DOI]
|
| [16] |
JIANG L, YIN M, XU J, et al. The transcription factor Bach1 suppresses the developmental angiogenesis of zebrafish[J]. Oxid Med Cell Longev, 2017, 2017: 2143875.
[URI]
|
| [17] |
JIANG L, YIN M, WEI X, et al. Bach1 represses Wnt/β-catenin signaling and angiogenesis[J]. Circ Res, 2015, 117(4): 364-375.
[DOI]
|