Signed weighted gene co-expression network analysis of transcriptional regulation in murine embryonic stem cells

Abstract

Background: Recent work has revealed that a core group of transcription factors (TFs) regulates the key characteristics of embryonic stem (ES) cells: pluripotency and self-renewal. Current efforts focus on identifying genes that play important roles in maintaining pluripotency and self-renewal in ES cells and aim to understand the interactions among these genes. To that end, we investigated the use of unsigned and signed network analysis to identify pluripotency and differentiation related genes.

Results: We show that signed networks provide a better systems level understanding of the regulatory mechanisms of ES cells than unsigned networks, using two independent murine ES cell expression data sets. Specifically, using signed weighted gene co-expression network analysis (WGCNA), we found a pluripotency module and a differentiation module, which are not identified in unsigned networks. We confirmed the importance of these modules by incorporating genome-wide TF binding data for key ES cell regulators. Interestingly, we find that the pluripotency module is enriched with genes related to DNA damage repair and mitochondrial function in addition to transcriptional regulation. Using a connectivity measure of module membership, we not only identify known regulators of ES cells but also show that Mrpl15, Msh6, Nrf1, Nup133, Ppif, Rbpj, Sh3gl2, and Zfp39, among other genes, have important roles in maintaining ES cell pluripotency and self-renewal. We also report highly significant relationships between module membership and epigenetic modifications (histone modifications and promoter CpG methylation status), which are known to play a role in controlling gene expression during ES cell self-renewal and differentiation.

Conclusion: Our systems biologic re-analysis of gene expression, transcription factor binding, epigenetic and gene ontology data provides a novel integrative view of ES cell biology.

Figure 1

Network Connection Strength Versus Expression Correlation. Network adjacency (y-axis) versus correlation (x-axis) for an unweighted network (black step function with τ = 0.8) and weighted networks (dashed lines corresponding to different powers, β) in an unsigned network (a) and a signed network (b). Note that cor(x_i, x_j) = -1 leads to adjacency = 0 in the signed network. The weighted network preserves the continuous nature of the co-expression information while an unweighted network dichotomizes the correlation.

Figure 2

Unsigned and Signed Mouse ES cell Networks in Ivanova et al. a, left, Dendrogram of the unsigned network of the Ivanova et al (2006) data set with color bands below indicating module membership for the unsigned network (U) and the signed network (S). a, right, heat map for visualizing standardized gene expressions (rows) across samples (columns) for genes in the turquoise module in the unsigned network. b, left, Dendrogram of the signed network of the Ivanova et al (2006) data set with color bands below indicating module membership for the signed network (S) and the unsigned network (U). b, right, heat map of expression profiles across samples for genes in the turquoise, black, and blue modules in the signed network. Note, modules are not scaled to reflect the number of genes in each module. c, scatter plot of module membership, k_ME, (x-axis) plotted against gene significance, GS, (y-axis) for the black and blue modules in the signed network with known ES cell regulators and differentiation genes labelled.

Figure 3

Relating Module Membership to Epigenetic Regulation. The top 1000 genes with highest module membership in the black module (top row) and blue module (bottom row) are related to 3 epigenetic variables (corresponding to the 3 columns). The y-axis reports the proportion of top 1000 genes that are known to belong to the group of genes defined on the x-axis. Histone H3K4me3 trimethylation status is abbreviated K4, H3K27me3 trimethylation statys is abbreviated by K27. Promoters that are both H3K4 and H3K27 trimethilated in ES cells (denoted K4&K27) are thought to poise key developmental genes for activation upon differentiation [50,51]. Note that genes with promoter CpG methylation are significantly (p = 2.0 × 10^-14) under-enriched with respect to the top 1000 black module genes.

Figure 4

Unsigned and Signed Networks of the Zhou et al ES Expression Data. a, left, Dendrogram of the unsigned network of the Zhou et al (2007) data with color bands below indicating module membership for the unsigned network (U) and the signed network (S). a, right, A heat map shows microarray expression profiles accross samples for genes in the blue and black modules in the unsigned network. b, left, Dendrogram of the signed network of the Zhou et al (2007) data with color bands below indicating module membership for the signed network (S) and the unsigned network (U). b, right, heatmap of expression profiles (rows) across samples (columns) for genes in the blue and black modules in the signed network.

Figure 5

Expression Changes Versus Module Membership in the Black and Blue Modules (Zhou et al). Module membership, k_ME, is plotted against log₂ expression fold change (FC) for the black (a) and blue (b) modules in the unsigned network of the Zhou et al (2007) data. FC is the ratio between the average expression in Oct4 positive samples and Oct4 negative microarray samples. Known ES cell regulators are labeled. Genes are colored by module membership in the signed network. c and d are analogous to a and b but module membership is with regard to the signed black and blue modules.

Figure 6

Comparison of Genes Ranked by Network Connectivity and Differential Expression in the Ivanova et al data set. Ingenuity Pathway Analysis of functional enrichments in the set of genes ranked within the top 1000 by Student's t-test and k_ME and yet do not overlap with each other. Venn diagrams show the amount of gene overlap between the top 1000 black (pluripotency) module genes and the top 1000 genes most significantly down-regulated upon Oct4 RNAi (left); gene overlap between the top 1000 blue (differentiation) module genes and the 1000 genes most significantly up-regulated with Oct4 RNAi (right). Significance of differential expression was determined using Student's t-statistic. p-values have been corrected for multiple hypothesis tests (Benjamini-Hochberg). Only significantly enriched functional groups are shown.

Figure 7

Transcriptional Regulators related to Pluripotency and Differentiation in the Zhou Network. TF and Suz12 binding in the promoter regions of highly connected genes related to ES cell pluripotency and self-renewal with GO terms of transcriptional regulation or chromatin structure. Genes are listed by black k_ME (positive, left and negative, right) along with their corresponding significance level (log₁₀ of the Bonferroni corrected p-value generated by a correlation test). Binding data from Chen et al (2008), Boyer et al (2006), and Loh et al (2006), are marked in blue (bound) and beige (unbound). For Oct4, Sox2, and Suz12, where binding is given by two studies, binding will be blue if it is found in both studies and light blue if found in only one.

Figure 8

Non-Transcriptional Regulators Related to Pluripotency and Differentiation. TFs and Suz12 binding of highly connected genes related to ES cell pluripotency and self-renewal lacking GO terms for transcriptional regulation or chromatin structure. Genes are tablulated in the same format as Figure 7.

Figure 9

Pluripotency Transcriptional Regulators that are not Bound by the Core TF Machinery in ES Cells. Genes related to ES cell pluripotency and self-renewal with GO terms of transcriptional regulation or chromatin structure and little pluripotency TF binding. Genes are listed by black k_ME along with their corresponding significance level (log₁₀ of the Bonferroni corrected p-value generated by a correlation test). Binding data from Chen et al, Boyer et al and Loh et al, are marked in blue (bound) and beige (unbound). For Suz12, light blue indicates that a genes is called bound in Chen et al or Boyer et al but not in both.