BACKGROUND
embryonic stem  cells have two important characteristics: pluripotency, the ability to differentiate into any type of cell in the body, and self-renewal, the ability to replicate indefinitely.
as such, they have tremendous therapeutic potential for regenerative medicine  <cit> .
current work focuses on understanding and extending the network of genes that controls these key characteristics  <cit> .
these efforts identified es cell-specific transcription factors  that are differentially expressed between es cells and differentiated cells .
several studies have identified the targets of these tfs and the mechanism by which they regulate them  <cit> .
highly differentially expressed tfs  have been found capable of reprogramming fibroblasts to a pluripotent state  <cit> .
while standard differential expression analysis techniques have led to remarkable discoveries they ignore the strong correlations that may exist between gene expression profiles.
as a consequence, the user of a standard marginal analysis can drown in information but starve in knowledge.
this is especially true when considering es cells where many genes change expression during differentiation.
for example, in a data set from zhou et al  <dig>  which we consider below, more than  <dig> genes were highly differentially expressed .
it is difficult to further prioritize these genes and to learn the underlying biological pathways.
in contrast, co-expression networks, also referred to as 'association,' 'correlation,' or 'influence' networks  <cit> , realize that genes can be highly correlated and thus can be grouped into large clusters .
for example, our network analysis of the same data organizes the genes into only  <dig> large modules.
next our module-centric analysis focuses on understanding the modules and their key regulators.
since it applies significance testing to the level of modules, co-expression network analysis may greatly alleviate the multiple testing problem that plagues standard gene-centric methods  <cit> .
gene co-expression network methods have been successfully applied in a variety of different settings  <cit> .
in this article, we demonstrate that a co-expression network analysis of stem cell data sets provides novel biological insights that cannot be found using conventional techniques.
using external data , we also contrast the performance of signed and unsigned network construction methods.
we find that signed co-expression network analysis performs best in this stem cell application.
we identify pluripotency and differentiation related co-expression modules and novel es cell regulators.
RESULTS
constructing signed co-expression networks
we first define a gene co-expression similarity measure which is used to define the network.
we denote the gene co-expression similarity measure of a pair of genes i and j by sij.
many co-expression studies use the absolute value of the pearson correlation as an unsigned co-expression similarity measure,
where gene expression profiles xi and xj consist of the expression of genes i and j across multiple microarray samples.
however, using the absolute value of the correlation may obfuscate biologically relevant information, since no distinction is made between gene repression and activation.
in contrast, in signed networks the similarity between genes reflects the sign of the correlation of their expression profiles.
to define a signed co-expression measure between gene expression profiles xi and xj, we use a simple transformation of the correlation:
as the unsigned measure , the signed similarity  takes on a value between  <dig> and  <dig>  note that the unsigned similarity between two oppositely expressed genes  = -1) equals  <dig> while it equals  <dig> for the signed similarity.
similarly, while the unsigned co-expression measure of two genes with zero correlation remains zero, the signed similarity equals  <dig> .
next, an adjacency matrix , a = , is used to quantify how strongly genes are connected to one another.
a is defined by thresholding the co-expression similarity matrix s = .
'hard' thresholding  the similarity measure s results in an unweighted gene co-expression network.
specifically an unweighted network adjacency is defined to be  <dig> if sij > τ and  <dig> otherwise, i.e.
two genes are considered connected if their similarity measure is above a given threshold τ, and are considered separated otherwise.
because hard thresholding encodes gene connections in a binary fashion, it can be sensitive to the choice of the threshold and result in the loss of co-expression information  <cit> .
the continuous nature of the co-expression information can be preserved by employing soft thresholding, which results in a weighted network.
specifically, we use a continuous measure to assess their connection strength:
where the power β is the thresholding parameter.
as a default we use β =  <dig> and β =  <dig> for unsigned and signed networks, respectively.
alternatively, β and be chosen using the scale-free topology criterion  <cit> .
since log = β × log, the weighted network adjacency is linearly related to the co-expression similarity on a logarithmic scale.
figure  <dig> shows the resulting adjacencies after applying the co-expression similarity measures and thresholding.
note that a high power β transforms high similarities into high adjacencies, while pushing low similarities towards  <dig>  since this soft-thresholding procedure leads to weighted adjacency matrix, the ensuing analysis is referred to as weighted gene co-expression network analysis or wgcna  <cit> .
a major step in our module centric analysis is to cluster genes into network modules using a network proximity measure.
roughly speaking, a pair of genes has a high proximity if it is closely interconnected.
we will use the convention that the maximal proximity between two genes is  <dig> and the minimum proximity is  <dig>  specifically, we define the proximity as the topological overlap measure   <cit>  which can also be defined for weighted networks  <cit> .
the tom combines the adjacency of two genes and the connection strengths these two genes share with other "third party" genes .
the tom is a highly robust measure of network interconnectedness .
this proximity is used as input of average linkage hierarchical clustering.
modules are defined as branches of the resulting cluster tree  <cit> .
this module detection procedure has been used in many applications  <cit>  and a comparison to alternative procedures is beyond the scope of this article.
we find it convenient to summarize the gene expression profiles of a given module with the module eigengene, which can be considered as the best summary of the standardized module expression data  <cit> .
the module eigengene of a given module is defined as the first principal component of the standardized expression profiles .
quantifying module membership
to identify possible regulators within a given module, we looked for highly connected intramodular hub genes, i.e.
genes that have strong connections within the module.
in our effort to find these genes, we examined two types of connectivity measures, which can be applied relative to any module q. the first connectivity measure is intramodular connectivity  defined as
where n is the number of genes in the qth module.
in the case of an unweighted network,  simply counts the number of connections to gene i within the qth module.
intramodular connectivity can be interpreted as a measure of module membership: the higher the intramodular connectivity, the more centrally located the gene is in the module and the more certain is its membership with regard to this module.
in signed networks, these highly connected hub genes may up-regulate adjacent genes since they are positively correlated with them, while in unsigned networks they may activate or repress their neighboring genes.
the second connectivity measure is the module eigengene based connectivity,  , defined as
where e is the eigengene of the qth module  and xi is the expression profile of the gene i. we denote modules by colors.
for example,  denotes the module membership measure of the i-th gene with regard to the blue module.
module eigengene based connectivity has several advantages over intramodular connectivity: first, it is naturally scaled to take on values between - <dig> and 1; second, one can use a correlation test to calculate a corresponding p-value for a gene's module membership; third it can be used in signed networks to identify genes that are anti-correlated with a given module eigengene , and fourth, kme can be computed for any gene on the array .
in practice, we found that intramodular and module eigengene based connectivity are highly correlated .
a priori, the connectivity measures defined in equations  <dig> and  <dig> are quite different.
but we show in the methods section that a simple theoretical relationship between them can be derived in the context of a signed co-expression module.
due to its advantages, we used the module eigengene based connectivity  as the measure of module membership in our applications.
signed wgcna identifies pluripotency related modules in ivanova et al  data set
we generated unsigned and signed co-expression networks to analyze over  <dig>  genes measured across  <dig> expression arrays from data published in ivanova et al   <cit> .
this data set contains expression profiles of es cells individually depleted for the transcription factors oct <dig>  nanog, sox <dig>  esrrb, and tbx <dig> by rna interference .
the data set also includes expression profiles for rnai knock downs of tcl <dig> a co-activator of akt kinase, and an est , along with expression profiles of control es cells carrying an empty rnai vector and of es cells differentiated by retinoic acid .
each of these treatments was sampled over approximately eight days.
to compare the performances of unsigned and signed wgcna in identifying gene groups that are important for the regulation of the pluripotent state, we defined gene modules in unsigned and signed networks and assessed module function and importance by determining gene ontology terms associated with each module and examining module membership of genes known to play a role in es cells.
in addition, we analyzed how genes of a given module are bound by chromatin regulators or pluripotency tfs by incorporating independent promoter binding information.
functional enrichment with regard to known es cell related genes
next we used external data to further study the gene modules defined by the networks and reveal their functional roles.
we used two different strategies for this evaluation: first we assigned transcription factors and other regulators with known roles in pluripotency, self-renewal or differentiation to modules  <cit>  and second, we incorporated genome-wide binding data for transcription factors and other regulators implicated in es cell regulation in order to determine if these modules contain genes that are directly controlled by es cell related tfs or differentiation suppressors  <cit> .
many genes known to maintain the pluripotent state of es cells are found in the black module in the signed network.
we defined a measure of gene significance  as the t-statistic from the paired student's t-test of expression in control rnai samples and es cell samples with rnai knock down of oct <dig> .
figure 2c shows gs plotted against its module eigengene based connectivity, kme, in the black and blue modules of the signed network with marker genes labeled.
since the signed module membership kme is defined as the correlation between a gene expression profile and the module eigengene, its values lie between - <dig> and  <dig> with values near  <dig> signifying strong module membership to the corresponding signed module.
figure 2c shows a strong linear relationship between kme and gs in the black module .
as expected, most of the genes whose rnai knock down induced es cell differentiation in ivanova et al  <cit>  belong to the black module .
oct4's high connectivity  makes it a hub gene in the black module, consistent with its known role as a master regulator of the pluripotent state.
furthermore, many genes that are known to be highly expressed in es cells are also in the black module .
klf <dig> is one of the four tfs that can reprogram differentiated cells into a pluripotency state  <cit> .
utf <dig> interacts with oct <dig>  affects chromatin regulation in es cells, and has recently been shown to improve reprogramming efficiency  <cit> .
phc <dig> is a polycomb group  protein.
pcg proteins repress genes that become active upon differentiation of es cells by mediating histone h <dig> lysine  <dig> tri-methylation and histone h2a ubiquitination  <cit> .
the blue module contains gata <dig> and gata <dig>  which are both highly connected .
these tfs are markers of es cell differentiation, particularly into endoderm.
below we provide further evidence that the black and blue modules are related to pluripotency and differentiation respectively.
module enrichment with regard to known es cell regulators
we incorporated genome-wide binding data for tfs  and other regulators  implicated in the maintenance of pluripotency and self-renewal, which were obtained by chromatin immunoprecipitation  and massive parallel sequencing  by chen et al   <cit> .
oct <dig>  sox <dig>  nanog, smad <dig>  and stat <dig> are referred to as the oct <dig> group of tfs, as they have been shown to often co-bind genomic regions; cmyc, nmyc, e2f <dig>  and zfx are referred to as the cmyc group of tfs because they also co-bind genomic regions  <cit> .
together tfs in the oct <dig> and cmyc group are thought to activate expression of genes involved in pluripotency and self-renewal.
suz <dig>  is a subunit of the histone h3k <dig> methyltransferase pcg protein complex, which represses genes that are activated upon differentiation  <cit> .
to determine if binding by the two tf groups and suz <dig> occurs more often in certain modules we computed binding enrichment for each module.
enrichment is defined as the odds ratio, that is the probability of a gene being bound by a particular tf or tf complex for genes in a given module divided by the probability of being bound for genes not located in the module.
table  <dig> shows module enrichment of genes bound by tfs in the oct <dig> group, the cmyc group, and suz <dig> for modules in the unsigned and signed networks.
a gene is called bound by the oct <dig> or cmyc groups if it is bound by at least  <dig> of the  <dig> and  <dig> of the  <dig> tfs in each tf group, respectively .
in agreement with the notion that the black module found in the signed network contains genes that are implicated in pluripotency, this module is strongly enriched with genes bound by tfs in the oct <dig> and cmyc groups and under-enriched for binding by suz <dig>  specifically, the proportion of genes bound by the oct <dig> group in the black module is almost twice the proportion of genes bound in the general population .
similarly genes in the black module are almost twice as likely to be bound by tfs in the cmyc group  and are almost three times less likely to be bound by suz <dig> .
these enrichments further support the idea that the black module is a pluripotency and self-renewal module.
other modules, like the blue, brown, and cyan, are enriched for suz <dig> bound genes and under-enriched for genes bound by the two tf groups.
the blue module is the most significantly enriched for suz <dig> binding further supporting the idea that genes in this module are involved in es cell differentiation.
for modules preserved in both unsigned and signed networks, table  <dig> shows that enrichment values are generally more significant in the signed network.
similar enrichments can be seen for oct <dig> and nanog binding from loh et al  and polycomb group  protein binding from boyer et al  .
the incorporation of binding data suggests that signed wgcna better separates genes into modules based on function and regulation.
specifically, the black module can be considered a pluripotency/self-renewal module, while the blue, brown, and cyan modules can be considered differentiation modules in the signed network.
it is important to reiterate that the pluripotency/self-renewal black module was only found using signed wgcna.
enrichment for binding by the oct <dig> and cmyc tf groups and suz <dig>
each enrichment score  is followed by its corresponding p-value  calculated from the hyper-geometric distribution.
p-values are uncorrected for multiple comparisons.
epigenetic regulation and module membership
recent studies suggest that chromatin structure and epigenetic modifications, like histone modification and dna methylation, play a role in controlling gene expression during es cell self-renewal and differentiation  <cit> .
for example, gene repression by the pcg protein complex via histone h <dig> lysine  <dig> trimethylation  is required for es cell self-renewal and pluripotency  <cit> .
to understand how epigenetic variables contribute to the regulation of es cells we studied the relationship of the pluripotency and differentiation modules with es cell h3k <dig> and h3k <dig> trimethylation, dna methylation, and cpg promoter content from previously published data sets  <cit> .
we related the epigenetic variables to module membership in the black  and blue module .
specifically, we determined what proportion of genes with a given epigenetic mark  are also in the top  <dig> genes with the highest  .
below we show that our findings are highly robust with respect to the number of selected module genes .
h3k4me <dig> is associated with gene activation, whereas h3k27me <dig> is known to silence genes.
figure  <dig> shows that genes with h3k <dig> trimethylation in es cells contain significantly more black module genes than genes without this marker  which may reflect the active role the pluripotency module plays in es cells.
interestingly, genes with h3k <dig> trimethylation or bivalent methylation contain significantly  more blue module genes than other gene classifications.
promoters that are both h3k <dig> and h3k <dig> trimethilated in es cells  are thought to poise key developmental genes for activation upon differentiation  <cit> .
mammalian gene promoters are known to fall into one of at least two major classes: 1) cpg-rich promoters are associated with both ubiquitously expressed 'housekeeping' genes, and genes with more complex expression patterns, particularly those expressed during embryonic development and 2) cpg-poor promoters are generally associated with highly tissue-specific genes.
to understand the role of cpg content in our modules we analyzed three cpg content classifications from mikkelsen et al: high , low , and intermediate .
figure  <dig> shows that hcp genes contain significantly more black module genes  and significantly more blue module genes  than icp or lcp genes.
the lcps are known to have a very different trimethylation pattern than the hcps.
few  of lcps have significant h3k4me <dig> in es cells and virtually none have h3k27me <dig>  hcps and lcps are subject to distinct modes of regulation.
in es cells, all hcps seem to be targets of trithorax group activity, and may therefore drive transcription unless actively repressed by pcg proteins.
in contrast, lcps seem to be inactive by default, independent of repression by pcg proteins, and may instead be selectively activated by cell-type- or tissue-specific factors  <cit> .
variance in kme explained by epigenetic variables
the above results demonstrate highly significant relationships between module membership and epigenetic variables.
in the following, we probe deeper and determine the proportion of variance in module membership that can be explained by the epigenetic variables.
using analysis of variance, we can determine what proportion of the variation in eigengene-based connectivity kme can be explained by the different epigenetic variables.
as can be seen from table  <dig>  the epigenetic variables explain only  <dig> % of the variation in  and  <dig> % of the variation in .
for , histone trimethylation status  and cmyc complex binding  explain most of the variation.
for , histone trimethylation status  and cpg class  explain most of the variation.
in summary, we find highly significant but relatively weak relationships between module membership and epigenetic variables.
this analysis of variance table reports which epigenetic variables and tf binding data have a significant effect on   and  .
for each variable , the columns report the degrees of freedom, the sums of squares, the proportion of total variance explained by the variable and the corresponding f-test p-value.
note that histone trimethylation status is the most significant source of variation for both  and .
signed wgcna identifies a pluripotency module in data from zhou et al
to further investigate wgcna's ability to discover functionally important groups of genes, we turned to an independent data set from zhou et al   <cit> .
in this study, es cells were removed from feeder cells and leukemia inhibitory factor  to induce differentiation.
during the course of differentiation, cells were separated based on expression of an oct <dig> green fluorescent protein  reporter gene.
multiple samples were taken from undifferentiated es cells and cells sorted at days  <dig>   <dig>   <dig>  and  <dig> for high and low oct <dig> expression.
as before, we first identified gene modules via signed and unsigned methods and then related module membership to external data.
in the following we show that a pluripotency/self-renewal and a differentiation module can be found in this new data set.
for consistency between data sets, we have colored these modules black and blue, respectively.
cluster tree comparison of unsigned and signed networks
similar to figure  <dig>  figure  <dig> shows the dendrograms of the unsigned and signed networks with color bands indicating module membership in the unsigned and signed networks.
the heat maps of the expression profiles of genes in the blue and black modules appear the same across the unsigned and signed networks.
the dendrograms, however, reveal network differences.
because module genes roughly stay together, we note that the red, black, blue, and yellow modules are partially preserved in the signed network.
however, the black and blue modules are co-located on the same large branch in the unsigned network  but are separated onto two distinct branches in the signed network , suggesting a change in network topology.
therefore, depending on branch cutting methods used, the blue and black modules could easily be merged into a larger module in the unsigned network.
indeed, the structure of the dendrogram in figure 4a suggests that the black and blue modules  may comprise one large module.
for illustrative purposes, we have separated them into smaller modules.
in the signed network, on the other hand, the black and blue modules will remain separate regardless of the branch cutting technique employed.
module membership in unsigned and signed networks
a comparison of transcription factor binding enrichment in unsigned and signed networks
this conclusion is further supported by our analysis of tf binding patterns in the modules of the zhou et al network.
table  <dig> shows enrichment for binding by the oct <dig> and cmyc groups of tfs and suz <dig> from chip-seq data from chen et al  <cit> .
in the signed network the black, brown, and red modules are enriched for binding by transcription factors in the oct <dig> and cmyc groups and are under-enriched for suz <dig> binding.
this is not surprising since the black, red, and brown modules are co-located on the same large branch in the signed network , which shows that they have similar expression profiles and share many of the same gene connections.
the blue module, on the other hand, is under-enriched for genes bound by the oct <dig> group tfs.
these results were also found when using oct <dig>  nanog, and pcg binding data from two other studies  and boyer et al ) .
the identification of the brown module and the increased enrichment of the black module in the signed network show that signed wgcna is better at grouping genes into modules based on transcriptional regulation.
furthermore, the black and blue modules are well separated in the signed network , while in the unsigned network branch cutting could erroneously group them into one large module .
module enrichments for binding by the oct <dig> and cmyc tf groups amd suz <dig>  each enrichment score  is followed by its corresponding p-value calculated from the hyper-geometric distribution .
functional enrichment analysis of the pluripotency and differentiation modules
given that there was significant overlap between the pluripotency modules and differentiation modules of the ivanova and zhou networks, we focused further analysis on the network constructed from the zhou et al data set as this network was based solely on differentiation induced expression changes.
we determined functional enrichment using the database for annotation visualization and integrated discovery   <cit> .
table  <dig> shows the significantly enriched go terms for genes with the 5% highest .
in agreement with the finding that the blue module is related to differentiation, genes within this module are significantly enriched for a functional group containing organ development, system development, and cell differentiation .
other highly enriched groups are involved in regulating protein localization  and membrane composition .
enriched go terms of genes with the highest 5% black and blue kme's.
highly connected genes in each functional group are given.
all kme values are highly significant .
we also used ingenuity pathway analysis, ipa, to compare functional enrichment in the pluripotency and differentiation modules .
additional file  <dig> shows that functional groups similar to those found using david are enriched in the black and blue module respectively.
cell cycle and dna replication, recombination, and repair are enriched in the black module compared to the blue module and skeletal, muscular, and cardiovascular system development are enriched in the blue module.
comparison to a standard differential expression analysis
here we compare some of our wgcna results with those of a standard differential expression analysis.
in figure 2c and figure  <dig> we showed that for some modules a strong relationship between module membership  and differential expression  can be observed.
in  <cit> , we provide a geometric description of modules for which such a relationship can be observed.
while a close relationship may exist between kme and a student t-test statistic, it does not imply that corresponding gene ranking procedures are equivalent.
here we compare signed wgcna to standard differential expression methods using three different approaches.
first, we show that a gene ranking based on kme is more consistent  than that based on the student t-test in our data.
specifically, we computed two gene rankings for the ivanova et al data set, one ranked by t-statistic and the other by connectivity to a module of interest.
we similarly computed two such rankings for the zhou et al data set and studied the overlap between the two data sets .
of the  <dig> genes most significantly down regulated upon differentiation in each data set  <dig> overlap .
however, when ranking genes by connectivity to each data set's pluripotency model there is an increase in overlap to  <dig> .
this increased consistency is also seen in genes up regulated upon differentiation where  <dig> genes overlap between the two data sets  when ranking by t-statistic and  <dig> genes overlap when ranking by connectivity to the differentiation modules .
a second approach for comparing gene rankings is to use the functional enrichment with regard to known gene ontologies.
to compare the abilities of kme versus the conventional t-statistic in identifying functionally interesting groups of genes we consider functional enrichment of genes found by one ranking method but not the other .
of the  <dig> genes most strongly connected to the ivanova et al pluripotency module ,  <dig> overlap with the  <dig> genes most significantly down regulated upon oct <dig> rnai in the same data set.
figure  <dig> shows the functional enrichment from ingenuity pathway analysis, ipa, of the  <dig> genes in each group that do not overlap.
genes found only by signed wgcna are significantly enriched for functions important to es cells like dna replication, recombination, and repair, cell cycle, cancer, protein synthesis, and rna post-transcriptional modification, which have previously been found using network methods  <cit> .
generally, genes identified by signed wgcna exhibit more significant enrichment of functional classifications than those found by standard differential analysis.
we similarly compare the  <dig> most highly connected genes in the differentiation  module and the  <dig> genes that are most significantly up regulated upon oct <dig> rnai in the ivanova network.
interestingly, only five genes overlap .
by examining those genes that do not overlap we see that ranking by connectivity yields greater significant enrichment for many functional groups important in es cell differentiation including organ development, tissue development, cell morphology etc.
similar analysis of pluripotency genes in zhou et al yields consistent results with dna replication, recombination, and repair being more enriched when ranked by connectivity  while analysis of highly connected genes in the differentiation module shows that differential analysis moderately out performs ranking by connectivity.
the differences in functional enrichment in the zhou et al data set are subtle given that there is more overlap between the two rankings .
this large overlap is likely due to the simplicity of the expression array samples which are filtered into only two groups, those that exhibit oct <dig> expression and those that do not.
meanwhile, signed wgcna is especially useful in ivanova et al where smaller overlap is caused by the complexity of the expression samples which are made of many different rnai treatments.
a third approach for comparing gene rankings is to use the enrichment with regard to epigenetic and transcriptional regulators.
in additional file  <dig> we relate different gene rankings to enrichment significance with regard to the following variables  histone h3k <dig> alone versus all others,  bivalent h3k <dig> & h3k <dig> versus all others  <cit> ,  high cpg class versus all others ,  promoter cpg methylation status  <cit> ,  oct  <dig> complex binding status,  cmyc complex binding status.
we report results for  <dig> different gene rankings using the ivanova data: the black and blue curve represent gene rankings according to  and , respectively.
the grey curve represents ranking according to a student t-test of differential expression.
additional file  <dig> shows that black and blue module genes can have very different enrichment results that tend to be very different from those of a standard analysis.
this analysis illustrates how module membership provides important complementary variables along the student t-test for understanding differences between genes.
the increased functional enrichment and improved consistency between data sets suggest that signed wgcna is a complementary method to standard differential analysis.
in practice, we recommend to use both kme and the student t-test to find highly differentially expressed intramodular hub genes.
pluripotency module genes involved in transcriptional regulation and chromatin structure
to gain a better understanding of the regulatory network involved in maintaining es cell pluripotency and self-renewal, we first looked at genes with high black  and go terms related to transcriptional regulation or chromatin structure.
figure  <dig> lists such genes and contains information about how these genes are bound by tfs in the oct <dig> and cmyc groups, and suz <dig> plus klf <dig>  esrrb, and tcfcp2l <dig>  the sign of black  allows us to distinguish genes that promote pluripotency from those that repress it.
as expected, many of the genes with high positive module membership measure () are known to participate in es cell regulation   <cit> .
importantly, four tfs capable of reprogramming differentiated cells into a es cell like state are in this list )  <cit> , confirming that the black module captures the known core transcriptional regulatory network responsible for maintaining a cell's stemness.
genes with high positive values of black module membership () but have not been implicated in maintaining pluripotency and self-renewal should therefore be strong candidates for further functional study.
such genes include msh <dig> , rbpj , zfp <dig> , and nrf <dig> .
similarly, many of the highly negatively connected genes in the black module are known to play a role in es cell differentiation .
because  and  are highly negatively correlated , genes with high positive  have negative  and vice versa.
highly negatively connected genes that are not known to be involved in es cell differentiation are also candidates for functional investigation.
some genes like maged <dig> are known to play a role in cell differentiation but have not been shown to be important in es cell differentiation specifically  <cit> .
other genes like lass <dig> have little known about their role in cell differentiation.
figure  <dig> also shows that genes positively connected to the black module tend to be bound by more es cell related tfs compared to the negatively connected genes, supporting the idea that these tfs bind and activate pluripotency and self-renewal genes  <cit> .
pluripotency module genes not involved in transcriptional regulation or chromatin structure
genes that are not involved in transcriptional regulation or chromatin modification  but have high average  are also of interest.
figure  <dig> lists such genes along with tf and suz <dig> binding information, connectivity, and fold change.
once again the importance of some genes has been validated, while others should be candidates for further research.
genes like dppa <dig>  dppa <dig>  and tcl <dig> are markers of pluripotency and nup <dig>  a nuclear pore complex subunit, has recently been shown to be necessary in the maintenance of pluripotency  <cit> .
nup <dig> highlights the usefulness of signed wgcna.
using the t-statistic from standard differential analysis nup <dig> is ranked 222th most significantly down regulated upon differentiation while using connectivity its rank moves to 28th.
other candidate genes include sh3gl <dig>  which binds lipids and proteins  <cit> , mrpl <dig>  a mitochondrial ribosomal protein, and ppif, involved in mitochondrial function and oxidative stress-induced cell death  <cit> .
of the negatively connected genes in figure  <dig>  ctsl  has recently been shown to cleave the histone h <dig> n-terminus during es cell differentiation  <cit> , while little is known about ctsz, also a cathepsin, while gnas and ctgf are differentiation genes  <cit> .
the high negative connectivity of uqcrh, a mitochondrial inner membrane protein  <cit> , along with the high positive connectivity of mrpl <dig> and ppif, confirms that mitochondrial regulation may be distinct in es cells  <cit>  and suggests that they may be important regulators of mitochondrial function in es cells.
pluripotency module genes that lack binding by known pluripotency tfs
to further investigate the role of these genes, we used motif scanning methods described in zhou et al   <cit>  to determine if the binding sites of these genes are contained in regions co-bound by tfs in the oct <dig> group or cmyc group in chip-seq data from chen et al .
we concentrated solely on lrh <dig>  and elk <dig> since their motifs have known position specific weight matrices while the other genes lack known motifs.
table  <dig> shows the enrichment and significance of motifs scanned.
sox <dig> and oct <dig> bind to the composite soxoct motif besides their own, stat <dig> binds the stat <dig> motif, and cmyc binds the ebox motif.
both the oct <dig> and cmyc groups' expected motifs are enriched.
for example, the soxoct motif has over three fold enrichment in regions bound by tfs in the oct <dig> group.
interestingly, the lrh <dig> motif is more enriched than the nanog motif in sequences bound by the oct <dig> group, which contains nanog binding by definition.
this reinforces the hypothesis that lrh <dig> co-binds regions bound by tfs in the oct <dig> group  <cit> .
furthermore, lrh <dig> sites are found in the promoter regions of pou5f <dig> , klf <dig>  dppa <dig>  and suz <dig> with pou5f <dig> having three separate sites.
these motif sites and lrh1's known importance in es cells, suggest that it may be an upstream regulator of these pluripotency factors and as such is a candidate for experimental validation  <cit> .
the elk <dig> motif is also significantly enriched in sequences bound by the cmyc group, thus ekl <dig> may co-regulate genes bound by tfs in this group.
sequences co-bound by all tfs of the oct <dig> group  or all tfs of the cmyc group  in the chen et al data set were scanned for relevant motifs, plus the lrh <dig>  and elk <dig> motifs.
there were  <dig> genes bound by all tfs in the oct <dig> group and  <dig> genes bound by all tfs in the cmyc group.
enrichments, computed by comparing against motif scans in control sequences, and p-values are shown below.
a geometric interpretation of signed wgcna modules
to understand how signed wgcna is better able to separate genes into functional modules in the ivanova data set, we plotted genes in the signed black or turquoise module relative to the unsigned turquoise module eigengene .
note that genes located in the black and turquoise modules in the signed network are clearly separated into two clusters.
because a module eigengene is defined as the first principle component of its module, it describes the main direction in which the module's gene expressions vary.
note that the signed module eigengenes are oriented in the direction of their clusters.
the direction of the unsigned turquoise module eigengene is more difficult to interpret.
because the turquoise module in the unsigned network contains two distinct signed modules , its module eigengene describes the variance between these two sub-modules and the variance within the larger sub-module, the signed turquoise.
as such, the unsigned turquoise module eigengene fails to quantify the true importance of highly connected genes in the signed black module.
for example, oct4's  in the unsigned turquoise module is - <dig>  while it is  <dig>  in the signed black module.
thus, oct <dig> is not identified as a hub gene in the unsigned network while it is clearly a hub gene in the signed network.
CONCLUSIONS
we show that a systems biology approach, which utilizes gene expression, transcription factor binding, genomic, epigenetic and gene ontology data, can be improved by accounting for the sign of co-expression relationships.
we also show that signed wgcna has advantages over standard differential expression methods.
specifically, signed wgcna has more consistent gene rankings between data sets , is better able to identify functionally enriched groups of genes , and its focus on module eigengenes circumvents the multiple testing problems that plague standard gene-based expression analysis.
below, we highlight several novel stem cell related genes that would not have been found using a standard differential expression analysis.
signed wgcna provides novel insight into murine es cell biology, which unsigned wgcna is unable to provide.
applying these signed methods to previously published data, we identified pluripotency and differentiation gene modules not found in unsigned networks or differential analysis.
the results of signed wgcna are robust as it identifies similar modules in independently published data sets.
we show that module eigengene based connectivity kme is valuable for annotating genes with regard to module membership and for identifying genes related to pluripotency and differentiation.
as a resource, we provide a module membership annotation for each gene with regard to the signed modules .
many current studies focus on the role transcriptional regulators play in es cell maintenance.
as expected, the pluripotency module is enriched with genes active in transcriptional regulation, e.g.
oct <dig>  sox <dig>  klf <dig>  nanog, jarid1b, jarid <dig>  nodal, tgif <dig>  and esrrb, and contains other genes expected to play a role in es cell function, such as dppa <dig> and dppa <dig>  the module also contains genes that have recently been shown to be necessary for maintaining the pluripotent state, nup <dig> and utf <dig>  <cit> .
interestingly, the pluripotency module contains genes with roles in two other pathways, dna repair and mitochondrial function, which are not found by standard differential analysis.
the enrichment for genes that respond to dna damage is not surprising given that es cells spend a larger portion of their cell cycle in s phase and have a shorter g <dig> phase than differentiated cells  <cit> .
an emphasis on accurate dna replication is expected since it helps es cells maintain a stable genome and prevents errors from being inherited by differentiated cells.
mitochondria in es cells may assist in the prevention of dna damage  <cit> .
during aerobic production of adenosine triphosphate , mitochondria leak superoxides leading to the creation of reactive oxygen species , which damage dna.
es cells, however, produce atp anaerobically and thus minimize the amount of dna damaging ros  <cit> .
es cells also have fewer mitochondria than differentiated cells and their mitochondria are smaller, have fewer cristae, lack dense matrices, and are perinuclearly located  <cit> .
our use of signed wgcna reveals that in addition to genes involved in transcriptional regulation, genes that prevent or repair dna damage are key to maintaining pluriotency and self-renewal.
using module eigengene based connectivity  we find that many known differentiation related genes are highly connected in the differentiation  module, cited <dig>  gata <dig>  and gata <dig>  along with ctsl, which has recently been shown to be active in differentiation  <cit> .
we also find that uqcrh, a gene involved in the electron transport chain, is highly connected in this module, lending support to the argument that es cell mitochondria differ from those in differentiated cells.
module eigengene based connectivity enabled us to identify novel candidate genes in the differentiation module, like uqcrh, that warrant experimental validation .
for the pluripotency module interesting candidate genes are msh <dig>  ppif, sh3gl <dig>  rbpj, elk <dig>  nrf <dig>  nup <dig>  mrpl <dig>  and zfp <dig> .
these genes lack significant fold change but are highly connected and thus would not be found using standard differential analysis.
using sequence data with motif analysis we confirm the importance of two genes, nr5a <dig> and elk <dig>  computationally.
we use gene ontology information and literature results to provide strong statistical evidence that these candidate genes are very promising and justify further biological study.
our article provides a resource in form of module based gene annotation tables that could form the starting point of future biological validation studies.
depending on their function, these candidate genes can be tested by rnai knock down, viral infection in order to increase the efficiency of reprogramming, or, if they bind dna, analyzing their binding sites.
our article demonstrates that signed wgcna not only identifies many well known es cell regulators; it also yields novel insights regarding es cell function.

 embryonic stem cells have two important characteristics: pluripotency, the ability to differentiate into any type of cell in the body, and self-renewal, the ability to replicate indefinitely. as such, they have tremendous therapeutic potential for regenerative medicine. current work focuses on understanding and extending the network of genes that controls these key characteristics. these efforts identified es cell-specific transcription factors that are differentially expressed between es cells and differentiated cells. several studies have identified the targets of these tfs and the mechanism by which they regulate them. highly differentially expressed tfs have been found capable of reprogramming fibroblasts to a pluripotent state. while standard differential expression analysis techniques have led to remarkable discoveries they ignore the strong correlations that may exist between gene expression profiles. as a consequence, the user of a standard marginal analysis can drown in information but starve in knowledge. this is especially true when considering es cells where many genes change expression during differentiation. for example, in a data set from zhou et al which we consider below, more than genes were highly differentially expressed. it is difficult to further prioritize these genes and to learn the underlying biological pathways. in contrast, co-expression networks, also referred to as'association,''correlation,' or'influence' networks, realize that genes can be highly correlated and thus can be grouped into large clusters. for example, our network analysis of the same data organizes the genes into only large modules. next our module-centric analysis focuses on understanding the modules and their key regulators. since it applies significance testing to the level of modules, co-expression network analysis may greatly alleviate the multiple testing problem that plagues standard gene-centric methods. gene co-expression network methods have been successfully applied in a variety of different settings. in this article, we demonstrate that a co-expression network analysis of stem cell data sets provides novel biological insights that cannot be found using conventional techniques. using external data, we also contrast the performance of signed and unsigned network construction methods. we find that signed co-expression network analysis performs best in this stem cell application. we identify pluripotency and differentiation related co-expression modules and novel es cell regulators.  constructing signed co-expression networks we first define a gene co-expression similarity measure which is used to define the network. we denote the gene co-expression similarity measure of a pair of genes i and j by sij. many co-expression studies use the absolute value of the pearson correlation as an unsigned co-expression similarity measure, where gene expression profiles xi and xj consist of the expression of genes i and j across multiple microarray samples. however, using the absolute value of the correlation may obfuscate biologically relevant information, since no distinction is made between gene repression and activation. in contrast, in signed networks the similarity between genes reflects the sign of the correlation of their expression profiles. to define a signed co-expression measure between gene expression profiles xi and xj, we use a simple transformation of the correlation: as the unsigned measure, the signed similarity takes on a value between and note that the unsigned similarity between two oppositely expressed genes = -1) equals while it equals for the signed similarity. similarly, while the unsigned co-expression measure of two genes with zero correlation remains zero, the signed similarity equals. next, an adjacency matrix, a =, is used to quantify how strongly genes are connected to one another. a is defined by thresholding the co-expression similarity matrix s =.'hard' thresholding the similarity measure s results in an unweighted gene co-expression network. specifically an unweighted network adjacency is defined to be if sij > τ and otherwise, i.e. two genes are considered connected if their similarity measure is above a given threshold τ, and are considered separated otherwise. because hard thresholding encodes gene connections in a binary fashion, it can be sensitive to the choice of the threshold and result in the loss of co-expression information. the continuous nature of the co-expression information can be preserved by employing soft thresholding, which results in a weighted network. specifically, we use a continuous measure to assess their connection strength: where the power β is the thresholding parameter. as a default we use β = and β = for unsigned and signed networks, respectively. alternatively, β and be chosen using the scale-free topology criterion. since log = β × log, the weighted network adjacency is linearly related to the co-expression similarity on a logarithmic scale. figure shows the resulting adjacencies after applying the co-expression similarity measures and thresholding. note that a high power β transforms high similarities into high adjacencies, while pushing low similarities towards since this soft-thresholding procedure leads to weighted adjacency matrix, the ensuing analysis is referred to as weighted gene co-expression network analysis or wgcna. a major step in our module centric analysis is to cluster genes into network modules using a network proximity measure. roughly speaking, a pair of genes has a high proximity if it is closely interconnected. we will use the convention that the maximal proximity between two genes is and the minimum proximity is specifically, we define the proximity as the topological overlap measure which can also be defined for weighted networks. the tom combines the adjacency of two genes and the connection strengths these two genes share with other "third party" genes. the tom is a highly robust measure of network interconnectedness. this proximity is used as input of average linkage hierarchical clustering. modules are defined as branches of the resulting cluster tree. this module detection procedure has been used in many applications and a comparison to alternative procedures is beyond the scope of this article. we find it convenient to summarize the gene expression profiles of a given module with the module eigengene, which can be considered as the best summary of the standardized module expression data. the module eigengene of a given module is defined as the first principal component of the standardized expression profiles. quantifying module membership to identify possible regulators within a given module, we looked for highly connected intramodular hub genes, i.e. genes that have strong connections within the module. in our effort to find these genes, we examined two types of connectivity measures, which can be applied relative to any module q. the first connectivity measure is intramodular connectivity defined as where n is the number of genes in the qth module. in the case of an unweighted network, simply counts the number of connections to gene i within the qth module. intramodular connectivity can be interpreted as a measure of module membership: the higher the intramodular connectivity, the more centrally located the gene is in the module and the more certain is its membership with regard to this module. in signed networks, these highly connected hub genes may up-regulate adjacent genes since they are positively correlated with them, while in unsigned networks they may activate or repress their neighboring genes. the second connectivity measure is the module eigengene based connectivity,, defined as where e is the eigengene of the qth module and xi is the expression profile of the gene i. we denote modules by colors. for example, denotes the module membership measure of the i-th gene with regard to the blue module. module eigengene based connectivity has several advantages over intramodular connectivity: first, it is naturally scaled to take on values between - and 1; second, one can use a correlation test to calculate a corresponding p-value for a gene's module membership; third it can be used in signed networks to identify genes that are anti-correlated with a given module eigengene, and fourth, kme can be computed for any gene on the array. in practice, we found that intramodular and module eigengene based connectivity are highly correlated. a priori, the connectivity measures defined in equations and are quite different. but we show in the methods section that a simple theoretical relationship between them can be derived in the context of a signed co-expression module. due to its advantages, we used the module eigengene based connectivity as the measure of module membership in our applications. signed wgcna identifies pluripotency related modules in ivanova et al data set we generated unsigned and signed co-expression networks to analyze over genes measured across expression arrays from data published in ivanova et al. this data set contains expression profiles of es cells individually depleted for the transcription factors oct nanog, sox esrrb, and tbx by rna interference. the data set also includes expression profiles for rnai knock downs of tcl a co-activator of akt kinase, and an est, along with expression profiles of control es cells carrying an empty rnai vector and of es cells differentiated by retinoic acid. each of these treatments was sampled over approximately eight days. to compare the performances of unsigned and signed wgcna in identifying gene groups that are important for the regulation of the pluripotent state, we defined gene modules in unsigned and signed networks and assessed module function and importance by determining gene ontology terms associated with each module and examining module membership of genes known to play a role in es cells. in addition, we analyzed how genes of a given module are bound by chromatin regulators or pluripotency tfs by incorporating independent promoter binding information. functional enrichment with regard to known es cell related genes next we used external data to further study the gene modules defined by the networks and reveal their functional roles. we used two different strategies for this evaluation: first we assigned transcription factors and other regulators with known roles in pluripotency, self-renewal or differentiation to modules and second, we incorporated genome-wide binding data for transcription factors and other regulators implicated in es cell regulation in order to determine if these modules contain genes that are directly controlled by es cell related tfs or differentiation suppressors. many genes known to maintain the pluripotent state of es cells are found in the black module in the signed network. we defined a measure of gene significance as the t-statistic from the paired student's t-test of expression in control rnai samples and es cell samples with rnai knock down of oct. figure 2c shows gs plotted against its module eigengene based connectivity, kme, in the black and blue modules of the signed network with marker genes labeled. since the signed module membership kme is defined as the correlation between a gene expression profile and the module eigengene, its values lie between - and with values near signifying strong module membership to the corresponding signed module. figure 2c shows a strong linear relationship between kme and gs in the black module. as expected, most of the genes whose rnai knock down induced es cell differentiation in ivanova et al belong to the black module. oct4's high connectivity makes it a hub gene in the black module, consistent with its known role as a master regulator of the pluripotent state. furthermore, many genes that are known to be highly expressed in es cells are also in the black module. klf is one of the four tfs that can reprogram differentiated cells into a pluripotency state. utf interacts with oct affects chromatin regulation in es cells, and has recently been shown to improve reprogramming efficiency. phc is a polycomb group protein. pcg proteins repress genes that become active upon differentiation of es cells by mediating histone h lysine tri-methylation and histone h2a ubiquitination. the blue module contains gata and gata which are both highly connected. these tfs are markers of es cell differentiation, particularly into endoderm. below we provide further evidence that the black and blue modules are related to pluripotency and differentiation respectively. module enrichment with regard to known es cell regulators we incorporated genome-wide binding data for tfs and other regulators implicated in the maintenance of pluripotency and self-renewal, which were obtained by chromatin immunoprecipitation and massive parallel sequencing by chen et al. oct sox nanog, smad and stat are referred to as the oct group of tfs, as they have been shown to often co-bind genomic regions; cmyc, nmyc, e2f and zfx are referred to as the cmyc group of tfs because they also co-bind genomic regions. together tfs in the oct and cmyc group are thought to activate expression of genes involved in pluripotency and self-renewal. suz is a subunit of the histone h3k methyltransferase pcg protein complex, which represses genes that are activated upon differentiation. to determine if binding by the two tf groups and suz occurs more often in certain modules we computed binding enrichment for each module. enrichment is defined as the odds ratio, that is the probability of a gene being bound by a particular tf or tf complex for genes in a given module divided by the probability of being bound for genes not located in the module. table shows module enrichment of genes bound by tfs in the oct group, the cmyc group, and suz for modules in the unsigned and signed networks. a gene is called bound by the oct or cmyc groups if it is bound by at least of the and of the tfs in each tf group, respectively. in agreement with the notion that the black module found in the signed network contains genes that are implicated in pluripotency, this module is strongly enriched with genes bound by tfs in the oct and cmyc groups and under-enriched for binding by suz specifically, the proportion of genes bound by the oct group in the black module is almost twice the proportion of genes bound in the general population. similarly genes in the black module are almost twice as likely to be bound by tfs in the cmyc group and are almost three times less likely to be bound by suz. these enrichments further support the idea that the black module is a pluripotency and self-renewal module. other modules, like the blue, brown, and cyan, are enriched for suz bound genes and under-enriched for genes bound by the two tf groups. the blue module is the most significantly enriched for suz binding further supporting the idea that genes in this module are involved in es cell differentiation. for modules preserved in both unsigned and signed networks, table shows that enrichment values are generally more significant in the signed network. similar enrichments can be seen for oct and nanog binding from loh et al and polycomb group protein binding from boyer et al. the incorporation of binding data suggests that signed wgcna better separates genes into modules based on function and regulation. specifically, the black module can be considered a pluripotency/self-renewal module, while the blue, brown, and cyan modules can be considered differentiation modules in the signed network. it is important to reiterate that the pluripotency/self-renewal black module was only found using signed wgcna. enrichment for binding by the oct and cmyc tf groups and suz each enrichment score is followed by its corresponding p-value calculated from the hyper-geometric distribution. p-values are uncorrected for multiple comparisons. epigenetic regulation and module membership recent studies suggest that chromatin structure and epigenetic modifications, like histone modification and dna methylation, play a role in controlling gene expression during es cell self-renewal and differentiation. for example, gene repression by the pcg protein complex via histone h lysine trimethylation is required for es cell self-renewal and pluripotency. to understand how epigenetic variables contribute to the regulation of es cells we studied the relationship of the pluripotency and differentiation modules with es cell h3k and h3k trimethylation, dna methylation, and cpg promoter content from previously published data sets. we related the epigenetic variables to module membership in the black and blue module. specifically, we determined what proportion of genes with a given epigenetic mark are also in the top genes with the highest. below we show that our findings are highly robust with respect to the number of selected module genes. h3k4me is associated with gene activation, whereas h3k27me is known to silence genes. figure shows that genes with h3k trimethylation in es cells contain significantly more black module genes than genes without this marker which may reflect the active role the pluripotency module plays in es cells. interestingly, genes with h3k trimethylation or bivalent methylation contain significantly more blue module genes than other gene classifications. promoters that are both h3k and h3k trimethilated in es cells are thought to poise key developmental genes for activation upon differentiation. mammalian gene promoters are known to fall into one of at least two major classes: 1) cpg-rich promoters are associated with both ubiquitously expressed'housekeeping' genes, and genes with more complex expression patterns, particularly those expressed during embryonic development and 2) cpg-poor promoters are generally associated with highly tissue-specific genes. to understand the role of cpg content in our modules we analyzed three cpg content classifications from mikkelsen et al: high, low, and intermediate. figure shows that hcp genes contain significantly more black module genes and significantly more blue module genes than icp or lcp genes. the lcps are known to have a very different trimethylation pattern than the hcps. few of lcps have significant h3k4me in es cells and virtually none have h3k27me hcps and lcps are subject to distinct modes of regulation. in es cells, all hcps seem to be targets of trithorax group activity, and may therefore drive transcription unless actively repressed by pcg proteins. in contrast, lcps seem to be inactive by default, independent of repression by pcg proteins, and may instead be selectively activated by cell-type- or tissue-specific factors. variance in kme explained by epigenetic variables the above results demonstrate highly significant relationships between module membership and epigenetic variables. in the following, we probe deeper and determine the proportion of variance in module membership that can be explained by the epigenetic variables. using analysis of variance, we can determine what proportion of the variation in eigengene-based connectivity kme can be explained by the different epigenetic variables. as can be seen from table the epigenetic variables explain only % of the variation in and % of the variation in. for, histone trimethylation status and cmyc complex binding explain most of the variation. for, histone trimethylation status and cpg class explain most of the variation. in summary, we find highly significant but relatively weak relationships between module membership and epigenetic variables. this analysis of variance table reports which epigenetic variables and tf binding data have a significant effect on and. for each variable, the columns report the degrees of freedom, the sums of squares, the proportion of total variance explained by the variable and the corresponding f-test p-value. note that histone trimethylation status is the most significant source of variation for both and. signed wgcna identifies a pluripotency module in data from zhou et al to further investigate wgcna's ability to discover functionally important groups of genes, we turned to an independent data set from zhou et al. in this study, es cells were removed from feeder cells and leukemia inhibitory factor to induce differentiation. during the course of differentiation, cells were separated based on expression of an oct green fluorescent protein reporter gene. multiple samples were taken from undifferentiated es cells and cells sorted at days and for high and low oct expression. as before, we first identified gene modules via signed and unsigned methods and then related module membership to external data. in the following we show that a pluripotency/self-renewal and a differentiation module can be found in this new data set. for consistency between data sets, we have colored these modules black and blue, respectively. cluster tree comparison of unsigned and signed networks similar to figure figure shows the dendrograms of the unsigned and signed networks with color bands indicating module membership in the unsigned and signed networks. the heat maps of the expression profiles of genes in the blue and black modules appear the same across the unsigned and signed networks. the dendrograms, however, reveal network differences. because module genes roughly stay together, we note that the red, black, blue, and yellow modules are partially preserved in the signed network. however, the black and blue modules are co-located on the same large branch in the unsigned network but are separated onto two distinct branches in the signed network, suggesting a change in network topology. therefore, depending on branch cutting methods used, the blue and black modules could easily be merged into a larger module in the unsigned network. indeed, the structure of the dendrogram in figure 4a suggests that the black and blue modules may comprise one large module. for illustrative purposes, we have separated them into smaller modules. in the signed network, on the other hand, the black and blue modules will remain separate regardless of the branch cutting technique employed. module membership in unsigned and signed networks a comparison of transcription factor binding enrichment in unsigned and signed networks this conclusion is further supported by our analysis of tf binding patterns in the modules of the zhou et al network. table shows enrichment for binding by the oct and cmyc groups of tfs and suz from chip-seq data from chen et al. in the signed network the black, brown, and red modules are enriched for binding by transcription factors in the oct and cmyc groups and are under-enriched for suz binding. this is not surprising since the black, red, and brown modules are co-located on the same large branch in the signed network, which shows that they have similar expression profiles and share many of the same gene connections. the blue module, on the other hand, is under-enriched for genes bound by the oct group tfs. these results were also found when using oct nanog, and pcg binding data from two other studies and boyer et al ). the identification of the brown module and the increased enrichment of the black module in the signed network show that signed wgcna is better at grouping genes into modules based on transcriptional regulation. furthermore, the black and blue modules are well separated in the signed network, while in the unsigned network branch cutting could erroneously group them into one large module. module enrichments for binding by the oct and cmyc tf groups amd suz each enrichment score is followed by its corresponding p-value calculated from the hyper-geometric distribution. functional enrichment analysis of the pluripotency and differentiation modules given that there was significant overlap between the pluripotency modules and differentiation modules of the ivanova and zhou networks, we focused further analysis on the network constructed from the zhou et al data set as this network was based solely on differentiation induced expression changes. we determined functional enrichment using the database for annotation visualization and integrated discovery. table shows the significantly enriched go terms for genes with the 5% highest. in agreement with the finding that the blue module is related to differentiation, genes within this module are significantly enriched for a functional group containing organ development, system development, and cell differentiation. other highly enriched groups are involved in regulating protein localization and membrane composition. enriched go terms of genes with the highest 5% black and blue kme's. highly connected genes in each functional group are given. all kme values are highly significant. we also used ingenuity pathway analysis, ipa, to compare functional enrichment in the pluripotency and differentiation modules. additional file shows that functional groups similar to those found using david are enriched in the black and blue module respectively. cell cycle and dna replication, recombination, and repair are enriched in the black module compared to the blue module and skeletal, muscular, and cardiovascular system development are enriched in the blue module. comparison to a standard differential expression analysis here we compare some of our wgcna results with those of a standard differential expression analysis. in figure 2c and figure we showed that for some modules a strong relationship between module membership and differential expression can be observed. in, we provide a geometric description of modules for which such a relationship can be observed. while a close relationship may exist between kme and a student t-test statistic, it does not imply that corresponding gene ranking procedures are equivalent. here we compare signed wgcna to standard differential expression methods using three different approaches. first, we show that a gene ranking based on kme is more consistent than that based on the student t-test in our data. specifically, we computed two gene rankings for the ivanova et al data set, one ranked by t-statistic and the other by connectivity to a module of interest. we similarly computed two such rankings for the zhou et al data set and studied the overlap between the two data sets. of the genes most significantly down regulated upon differentiation in each data set overlap. however, when ranking genes by connectivity to each data set's pluripotency model there is an increase in overlap to. this increased consistency is also seen in genes up regulated upon differentiation where genes overlap between the two data sets when ranking by t-statistic and genes overlap when ranking by connectivity to the differentiation modules. a second approach for comparing gene rankings is to use the functional enrichment with regard to known gene ontologies. to compare the abilities of kme versus the conventional t-statistic in identifying functionally interesting groups of genes we consider functional enrichment of genes found by one ranking method but not the other. of the genes most strongly connected to the ivanova et al pluripotency module, overlap with the genes most significantly down regulated upon oct rnai in the same data set. figure shows the functional enrichment from ingenuity pathway analysis, ipa, of the genes in each group that do not overlap. genes found only by signed wgcna are significantly enriched for functions important to es cells like dna replication, recombination, and repair, cell cycle, cancer, protein synthesis, and rna post-transcriptional modification, which have previously been found using network methods. generally, genes identified by signed wgcna exhibit more significant enrichment of functional classifications than those found by standard differential analysis. we similarly compare the most highly connected genes in the differentiation module and the genes that are most significantly up regulated upon oct rnai in the ivanova network. interestingly, only five genes overlap. by examining those genes that do not overlap we see that ranking by connectivity yields greater significant enrichment for many functional groups important in es cell differentiation including organ development, tissue development, cell morphology etc. similar analysis of pluripotency genes in zhou et al yields consistent results with dna replication, recombination, and repair being more enriched when ranked by connectivity while analysis of highly connected genes in the differentiation module shows that differential analysis moderately out performs ranking by connectivity. the differences in functional enrichment in the zhou et al data set are subtle given that there is more overlap between the two rankings. this large overlap is likely due to the simplicity of the expression array samples which are filtered into only two groups, those that exhibit oct expression and those that do not. meanwhile, signed wgcna is especially useful in ivanova et al where smaller overlap is caused by the complexity of the expression samples which are made of many different rnai treatments. a third approach for comparing gene rankings is to use the enrichment with regard to epigenetic and transcriptional regulators. in additional file we relate different gene rankings to enrichment significance with regard to the following variables histone h3k alone versus all others, bivalent h3k & h3k versus all others, high cpg class versus all others, promoter cpg methylation status, oct complex binding status, cmyc complex binding status. we report results for different gene rankings using the ivanova data: the black and blue curve represent gene rankings according to and, respectively. the grey curve represents ranking according to a student t-test of differential expression. additional file shows that black and blue module genes can have very different enrichment results that tend to be very different from those of a standard analysis. this analysis illustrates how module membership provides important complementary variables along the student t-test for understanding differences between genes. the increased functional enrichment and improved consistency between data sets suggest that signed wgcna is a complementary method to standard differential analysis. in practice, we recommend to use both kme and the student t-test to find highly differentially expressed intramodular hub genes. pluripotency module genes involved in transcriptional regulation and chromatin structure to gain a better understanding of the regulatory network involved in maintaining es cell pluripotency and self-renewal, we first looked at genes with high black and go terms related to transcriptional regulation or chromatin structure. figure lists such genes and contains information about how these genes are bound by tfs in the oct and cmyc groups, and suz plus klf esrrb, and tcfcp2l the sign of black allows us to distinguish genes that promote pluripotency from those that repress it. as expected, many of the genes with high positive module membership measure () are known to participate in es cell regulation. importantly, four tfs capable of reprogramming differentiated cells into a es cell like state are in this list ), confirming that the black module captures the known core transcriptional regulatory network responsible for maintaining a cell's stemness. genes with high positive values of black module membership () but have not been implicated in maintaining pluripotency and self-renewal should therefore be strong candidates for further functional study. such genes include msh, rbpj, zfp, and nrf. similarly, many of the highly negatively connected genes in the black module are known to play a role in es cell differentiation. because and are highly negatively correlated, genes with high positive have negative and vice versa. highly negatively connected genes that are not known to be involved in es cell differentiation are also candidates for functional investigation. some genes like maged are known to play a role in cell differentiation but have not been shown to be important in es cell differentiation specifically. other genes like lass have little known about their role in cell differentiation. figure also shows that genes positively connected to the black module tend to be bound by more es cell related tfs compared to the negatively connected genes, supporting the idea that these tfs bind and activate pluripotency and self-renewal genes. pluripotency module genes not involved in transcriptional regulation or chromatin structure genes that are not involved in transcriptional regulation or chromatin modification but have high average are also of interest. figure lists such genes along with tf and suz binding information, connectivity, and fold change. once again the importance of some genes has been validated, while others should be candidates for further research. genes like dppa dppa and tcl are markers of pluripotency and nup a nuclear pore complex subunit, has recently been shown to be necessary in the maintenance of pluripotency. nup highlights the usefulness of signed wgcna. using the t-statistic from standard differential analysis nup is ranked 222th most significantly down regulated upon differentiation while using connectivity its rank moves to 28th. other candidate genes include sh3gl which binds lipids and proteins, mrpl a mitochondrial ribosomal protein, and ppif, involved in mitochondrial function and oxidative stress-induced cell death. of the negatively connected genes in figure ctsl has recently been shown to cleave the histone h n-terminus during es cell differentiation, while little is known about ctsz, also a cathepsin, while gnas and ctgf are differentiation genes. the high negative connectivity of uqcrh, a mitochondrial inner membrane protein, along with the high positive connectivity of mrpl and ppif, confirms that mitochondrial regulation may be distinct in es cells and suggests that they may be important regulators of mitochondrial function in es cells. pluripotency module genes that lack binding by known pluripotency tfs to further investigate the role of these genes, we used motif scanning methods described in zhou et al to determine if the binding sites of these genes are contained in regions co-bound by tfs in the oct group or cmyc group in chip-seq data from chen et al. we concentrated solely on lrh and elk since their motifs have known position specific weight matrices while the other genes lack known motifs. table shows the enrichment and significance of motifs scanned. sox and oct bind to the composite soxoct motif besides their own, stat binds the stat motif, and cmyc binds the ebox motif. both the oct and cmyc groups' expected motifs are enriched. for example, the soxoct motif has over three fold enrichment in regions bound by tfs in the oct group. interestingly, the lrh motif is more enriched than the nanog motif in sequences bound by the oct group, which contains nanog binding by definition. this reinforces the hypothesis that lrh co-binds regions bound by tfs in the oct group. furthermore, lrh sites are found in the promoter regions of pou5f, klf dppa and suz with pou5f having three separate sites. these motif sites and lrh1's known importance in es cells, suggest that it may be an upstream regulator of these pluripotency factors and as such is a candidate for experimental validation. the elk motif is also significantly enriched in sequences bound by the cmyc group, thus ekl may co-regulate genes bound by tfs in this group. sequences co-bound by all tfs of the oct group or all tfs of the cmyc group in the chen et al data set were scanned for relevant motifs, plus the lrh and elk motifs. there were genes bound by all tfs in the oct group and genes bound by all tfs in the cmyc group. enrichments, computed by comparing against motif scans in control sequences, and p-values are shown below. a geometric interpretation of signed wgcna modules to understand how signed wgcna is better able to separate genes into functional modules in the ivanova data set, we plotted genes in the signed black or turquoise module relative to the unsigned turquoise module eigengene. note that genes located in the black and turquoise modules in the signed network are clearly separated into two clusters. because a module eigengene is defined as the first principle component of its module, it describes the main direction in which the module's gene expressions vary. note that the signed module eigengenes are oriented in the direction of their clusters. the direction of the unsigned turquoise module eigengene is more difficult to interpret. because the turquoise module in the unsigned network contains two distinct signed modules, its module eigengene describes the variance between these two sub-modules and the variance within the larger sub-module, the signed turquoise. as such, the unsigned turquoise module eigengene fails to quantify the true importance of highly connected genes in the signed black module. for example, oct4's in the unsigned turquoise module is - while it is in the signed black module. thus, oct is not identified as a hub gene in the unsigned network while it is clearly a hub gene in the signed network.  we show that a systems biology approach, which utilizes gene expression, transcription factor binding, genomic, epigenetic and gene ontology data, can be improved by accounting for the sign of co-expression relationships. we also show that signed wgcna has advantages over standard differential expression methods. specifically, signed wgcna has more consistent gene rankings between data sets, is better able to identify functionally enriched groups of genes, and its focus on module eigengenes circumvents the multiple testing problems that plague standard gene-based expression analysis. below, we highlight several novel stem cell related genes that would not have been found using a standard differential expression analysis. signed wgcna provides novel insight into murine es cell biology, which unsigned wgcna is unable to provide. applying these signed methods to previously published data, we identified pluripotency and differentiation gene modules not found in unsigned networks or differential analysis. the results of signed wgcna are robust as it identifies similar modules in independently published data sets. we show that module eigengene based connectivity kme is valuable for annotating genes with regard to module membership and for identifying genes related to pluripotency and differentiation. as a resource, we provide a module membership annotation for each gene with regard to the signed modules. many current studies focus on the role transcriptional regulators play in es cell maintenance. as expected, the pluripotency module is enriched with genes active in transcriptional regulation, e.g. oct sox klf nanog, jarid1b, jarid nodal, tgif and esrrb, and contains other genes expected to play a role in es cell function, such as dppa and dppa the module also contains genes that have recently been shown to be necessary for maintaining the pluripotent state, nup and utf. interestingly, the pluripotency module contains genes with roles in two other pathways, dna repair and mitochondrial function, which are not found by standard differential analysis. the enrichment for genes that respond to dna damage is not surprising given that es cells spend a larger portion of their cell cycle in s phase and have a shorter g phase than differentiated cells. an emphasis on accurate dna replication is expected since it helps es cells maintain a stable genome and prevents errors from being inherited by differentiated cells. mitochondria in es cells may assist in the prevention of dna damage. during aerobic production of adenosine triphosphate, mitochondria leak superoxides leading to the creation of reactive oxygen species, which damage dna. es cells, however, produce atp anaerobically and thus minimize the amount of dna damaging ros. es cells also have fewer mitochondria than differentiated cells and their mitochondria are smaller, have fewer cristae, lack dense matrices, and are perinuclearly located. our use of signed wgcna reveals that in addition to genes involved in transcriptional regulation, genes that prevent or repair dna damage are key to maintaining pluriotency and self-renewal. using module eigengene based connectivity we find that many known differentiation related genes are highly connected in the differentiation module, cited gata and gata along with ctsl, which has recently been shown to be active in differentiation. we also find that uqcrh, a gene involved in the electron transport chain, is highly connected in this module, lending support to the argument that es cell mitochondria differ from those in differentiated cells. module eigengene based connectivity enabled us to identify novel candidate genes in the differentiation module, like uqcrh, that warrant experimental validation. for the pluripotency module interesting candidate genes are msh ppif, sh3gl rbpj, elk nrf nup mrpl and zfp. these genes lack significant fold change but are highly connected and thus would not be found using standard differential analysis. using sequence data with motif analysis we confirm the importance of two genes, nr5a and elk computationally. we use gene ontology information and literature results to provide strong statistical evidence that these candidate genes are very promising and justify further biological study. our article provides a resource in form of module based gene annotation tables that could form the starting point of future biological validation studies. depending on their function, these candidate genes can be tested by rnai knock down, viral infection in order to increase the efficiency of reprogramming, or, if they bind dna, analyzing their binding sites. our article demonstrates that signed wgcna not only identifies many well known es cell regulators; it also yields novel insights regarding es cell function. 
using a connectivity measure of module membership, we not only identify known regulators of es cells but also show that mrpl msh nrf nup ppif, rbpj, sh3gl and zfp 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, which are known to play a role in controlling gene expression during es cell self-renewal and differentiation. interestingly, we find that the pluripotency module is enriched with genes related to dna damage repair and mitochondrial function in addition to transcriptional regulation. specifically, using signed weighted gene co-expression network analysis, we found a pluripotency module and a differentiation module, which are not identified in unsigned networks. 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. recent work has revealed that a core group of transcription factors regulates the key characteristics of embryonic stem cells: pluripotency and self-renewal. to that end, we investigated the use of unsigned and signed network analysis to identify pluripotency and differentiation related genes. 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. 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.
 recent work has revealed that a core group of transcription factors regulates the key characteristics of embryonic stem 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.  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, 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 mrpl msh nrf nup ppif, rbpj, sh3gl and zfp 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, which are known to play a role in controlling gene expression during es cell self-renewal and differentiation.  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. 