BACKGROUND
relationships between genes and gene pathways can be inferred based on gene expression profiles across conditions, treatments, or samples obtained from microarray experiments  <cit> .
results from these studies can aid in the confirmation of previously known pathways and motivate the study of newly uncovered relationships among genes.
a bayesian network approach is well-suited to detect relationships between genes using an acyclic direct graph  <cit> .
needham et al.
provided an in depth primer on learning bayesian networks for computational biology  <cit> .
this approach has a solid theoretical foundation, offers a probabilistic framework to describe the variation typically observed in microarray data, accommodates missing data, and incorporates prior knowledge on gene relationships.
several implementations of bayesian networks to infer gene networks have been reported  <cit>  and two potential weaknesses have been identified.
first, some applications require the transformation of continuous gene expression data into discrete input data, and this can influence the resulting network or resulting in loss of information.
thus, analysis of data that has not been discretized is favored.
second, most bayesian network implementations assume standard binomial, multinomial probability or single gaussian probability density functions of gene expression across a wide range of conditions.
these distributions may fail to accommodate multimodal or skewed distributions associated with condition-dependent networks that exhibit changes in gene expression or gene relationships across conditions.
this is because gene network inference is typically based on gene expression measurements across a wide range of conditions.
mixture models can be used to address this limitation because these models can describe potentially complex distributions of gene expression across a wide range of conditions.
newman and leicht explored the capability of probabilistic mixture models to detect a very broad range of network structures without prior knowledge  <cit> .
ko et al.
presented preliminary results from the application of bayesian mixture algorithm to infer pathways  <cit> .
they presented a histone pathway, did not use parameter estimates to identify condition-dependent gene relationships, and did not use cross-validation to assess the adequacy of the inferred network.
in this study, gene networks were inferred using a mixture bayesian network approach.
the superiority of this approach was evidenced on the new insights into general and condition-dependent relationships between genes and gene expression profiles gained while addressing some of the limitations encountered in previous applications of bayesian network to infer gene networks.
we demonstrate that the integration of mixture models and bayesian network is well-suited to infer the structure of gene networks from continuous gene expression measurements across a wide range of conditions.
the estimated parameters were easy to interpret, and aided in the characterization of gene expression and co-expression profiles.
furthermore, functions of these parameter estimates allowed the identification of condition-dependent networks.
the mixture bayesian network approach was used to infer three pathways on three independent data sets.
the first prediction of the circadian rhythm pathway in honey bees and the first prediction of a cell communication pathway  in mouse embryos were obtained.
in addition, the mixture bayesian network approach was applied to a yeast cell-cycle data set commonly used in many gene network studies.
comparisons of the predicted and known pathways and benchmark tests confirmed the outstanding performance of the mixture bayesian network approach.
novel insights into general and condition-specific gene relationships and expression patterns were obtained, thus demonstrating the strength of the proposed approach to infer gene networks.
methods
mixture bayesian network model
bayesian networks can be described as a directed acyclic graph with nodes representing random variables or genes and directed edges representing the relationships between the nodes  <cit> .
given a set of genes {g <dig>  g <dig> ..., gn}, the corresponding bayesian network g is represented as the joint probability distribution over all genes in g or p. applying the narrow sense form of the markov property for a stochastic process, each gene is independent of non-descendant genes in the network, given the parent genes or nodes.
this conditional independence property allows the factorization of the joint probability distribution as the product of conditional probabilities.
the overall network is conceived as a set of gene subnetworks, each corresponding to a given gene node gj and its associated parent gene nodes a for j =  <dig> to n.
models that accommodate for potential changes in gene expression and co-expression patterns across conditions are necessary when gene expression data across multiple conditions is used to infer gene networks.
mixture models are a flexible and effective option to describe changes in gene co-expression patterns across conditions because the joint probability functions of sub-networks are modeled with a combination of mixture components.
the concept of integrating mixtures of gaussian densities into bayesian networks was introduced by davies and moore  <cit> , and applied to word and social networks by newman and leicht  <cit> .
under the bayesian network framework,
where n is the total number of genes  in the network  and a is the set of parent genes of child gene gj in the j-th sub-network.
the conditional probability density function of the j-th gene p), given the set of parent genes is expressed as the ratio between the joint probability density function of the parent and child genes in the j-th sub-network )), and the marginal probability density function of the parent genes in the j-th sub-network )).
the joint and marginal probability density functions are assumed to follow a multivariate gaussian distribution and a mixture model is used to describe possible sub-network structures.
the probability density function  for the j-th sub-network is represented with a mixture of kj multivariate gaussian distributions  each with a weight αjk.
here, xji denotes a pj dimensional vector including the i-th expression observation  of the child and parent genes in the j-th sub-network.
for example, if gj has three parents then pj =  <dig>  therefore, the joint probability density function for the j-th sub-network is
and
where .
for the joint probability density function, each component k of the mixture is described with a mean vector μjk of dimension pj and a variance-covariance matrix Σjk of dimension pj × pj.
the marginal probability density function of the  parent genes a is described in the same way as the joint probability density function but the dimensions of the i-th vector of observations , the mean vector  and the variance-covariance matrix  are , , and  × , respectively.
therefore, the overall likelihood of the data across all genes in the network g is:
an example of a bayesian gene sub-network is provided in additional file  <dig>  the relationship between the child and parent genes in each sub-network can be described by the weighted sum of the correlations between the parent geneand child gene estimated for each mixture component.
the weights correspond to the mixture component weights.
this computation was described by bland and altman  <cit>  and is presented in additional file  <dig>
overall network learning algorithm
the learning algorithm integrated two major components, a) the inference of each individual gene sub-network structure and b) the inference of the overall network structure.
the first component consisted of the identification of the sub-network for each gene , the estimation of the parameters of the mixture models using the expectation-maximization algorithm  <cit> , and the identification of the number of components of the mixture model best supported by the data.
the second component consisted of the combination of the individual gene sub-networks into an overall network , and the removal of cyclic relationships between genes in the overall network.
a detailed description of the network structure learning approach is presented in additional file  <dig>  the equations in the expectation-maximization steps corresponding to a gene sub-network follow the derivations of bilmes  <cit>  and are summarized in additional file  <dig>  the bayesian information criterion  was used to identify the mixture model best supported by the data  for each gene sub-network, and to evaluate the overall network  <cit> .
bayesian information criterion offers a good compromise between model adequacy and parsimony, minimizing the possibility of over-fitting the data with highly parameterized networks.
this feature of the bic is advantageous in gene network inference because of the large number of relationships and conditions that could potentially be considered.
newman and leicht used probabilistic mixture models to infer a wide range of network structures and noted that the convergence of the algorithm to the global maxima is not guaranteed  <cit> .
to address this potential problem, we followed the approach implemented by newman and leicht, and networks were inferred using a range of starting points for each mixture parameter  <cit> .
each network presented in this study had the highest likelihood over all sets of initial conditions or runs evaluated.
furthermore, davies and moore demonstrated  the suitability of mixtures of gaussian densities estimated with the em algorithm to infer bayesian networks using data from an extensive astronomical survey and a high-throughput biological cell assay  <cit> .
extensive evaluation of the performance of the mixture bayesian network approach to infer networks was presented by newman and leicht and davies and moore  <cit> .
in this study, we concentrated on the application of this methodology to infer gene networks using microarray data, and demonstrated how the interpretation of parameter estimates and function thereof offered new insights into the relationships between genes.
a unique advantage of the integration of the mixture models and bayesian network is the ability to uncover general and condition-specific gene networks.
the parameter estimates of each mixture component of a gene sub-network were used as indicators of the gene expression level , variation , and co-variation between genes .
in addition, the mixture parameter estimates were used to compute the odds that the observation corresponded to a particular mixture component relative to the other components.
each sample had a non-zero probability  to pertain to each of the mixture components.
for ease of interpretation and to overcome the limitation of using a single probability threshold across sub-networks with different number of mixture components, no threshold was used to assign samples to mixtures.
instead, the samples were assigned to the mixture component that had the highest odds.
the assignment of samples to the most likely mixture component allowed the identification of samples  that shared the same expression and co-expression patterns for a gene sub-network.
data
gene networks corresponding to three gene pathways  were inferred from three microarray gene expression data sets  using a mixture bayesian network approach.
the goal of these applications was to gain insights into general and condition-dependent relationships between genes and gene expression profiles, and demonstrate the flexibility and adequacy of the mixture bayesian network approach to characterize networks.
honey bee data set and circadian rhythm pathway
a detailed description of the honey bee microarray data set is presented by whitfield et al.
and rodriguez-zas et al.
<cit> .
briefly, gene expressions during behavioral maturation in the brains of honey bees  from two subspecies  rose in one of two host colonies  representing two different environments were measured.
within combinations of bee and host race , nurse bees were sampled on days  <dig>   <dig>   <dig>   <dig> and  <dig> after adult emergence  and forager bees  were sampled at day  <dig> after emergence.
the experimental design allowed the division of the honey bee data set into four independent bee-host subspecies-combination sub-data set .
thus, the gene network was inferred on the complete data set and also on each of the four sub-data set.
processing of the data from  <dig> spotted cdna microarrays included background subtraction and log2-transformation of fluorescence intensities, and removal of flagged spots or spots that did not surpass a minimum intensity threshold of  <dig>  the log <dig> intensity values were normalized using a loess transformation, and global dye and microarray effects were removed using the approach implemented beehive .
six genes assigned to the fruit fly circadian rhythm pathway in the kegg database  were present on the microarray platform.
the kegg identifier for the circadian rhythm pathway in the fruit fly is dme <dig> .
mouse embryo data set and adherens junction pathway
data from nine microarray experiments that evaluated the effect of multiple toxic agents on mouse  embryo gene expression levels were used to reconstruct the adherens junction pathway.
a detailed description of this data set is provided by rodriguez-zas et al.
<cit> .
briefly, the experimental design consisted of two types of direct comparisons with reverse labeling; comparisons between samples that received different treatments , and comparisons between treated and control samples.
a total of  <dig> microarrays from the same platform were available for analysis, and data processing and normalization followed the same protocols applied to the honey bee data set  <cit>  and implemented in beehive.
seven genes assigned to the mouse adherens junction pathway in the kegg database were present on the microarray platform.
the kegg identifier for the adherens junction pathway in the mouse is mmu <dig> .
yeast synchronization data set and cell cycle pathway
the mixture bayesian network approach was used to predict the yeast  cell-cycle network using gene expression data reported by  <cit> .
the data set consisted of gene expression measurements from  <dig> microarrays across six experiments.
the yeast experiments evaluated yeast cell synchronization via the arrest of cells by one of four conditions, followed by activation or release of the cells from the arresting condition.
the synchronization conditions were cdc15-, cdc28-, α factor- and size or elutrition- based synchronization, and samples were obtained at various time points after the removal of the last four arresting conditions.
in a separate study, the effects of inducing either the g <dig> cyclin cln3p  or the b-type cyclin clb2p  were examined, each in two experiments.
the normalized gene expression data is available at the stanford microarray database  and the experiment identifiers are:  <dig>   <dig>   <dig>   <dig>   <dig>   <dig> to  <dig>   <dig> and  <dig>  a network including  <dig> genes analyzed by friedman et al.
using multinomial and single gaussian distribution was inferred and results were compared  <cit> .
the kegg identifier for the yeast cell cycle pathway is sce <dig> .
validation of the inferred gene networks
the relationships between genes inferred by the mixture bayesian network approach were compared to relationships reported in the kegg database, biogrid  database, and saccharomyces genome database, or sgd , and in the literatures.
these databases were accessed on august  <dig> and provide complementary information.
the kegg database includes cause-effect or directional and non-directional relationships between genes depicted in pathway diagrams.
in directional relationships, the child and parent genes are identified, meanwhile in non-directional relationships, genes are linked but there is neither a parent  nor a child  gene.
the biogrid database lists all known interactions or undirected relationships between genes.
due to differences in the curation and annotation protocols of each database, biogrid can encompass more relationships than kegg, but these relationships do not have associated direction.
two strategies were used to confirm the reliability of the gene networks predicted by the mixture bayesian network approach.
first, a k-fold cross-validation approach was implemented.
the design of the circadian rhythm experiment allowed to infer the gene network in each of four independent data sets corresponding to the four bee-host race combinations that encompassed the same five behavioral maturation stages or ages.
assessment of the mixture bayesian network approach consisted on the comparison of the four inferred networks.
second, a permutation-based resampling approach was used to generate  <dig> randomized data sets with the same structure as the original circadian rhythm data set.
these data sets were analyzed using the mixture bayesian network approach and the predicted networks were compared to the observed network and the known pathway in the database.
not all genes present in the pathway databases had reporters in the microarray platforms studied and thus two types of gene relationships, direct and indirect  were inferred.
direct relationships were detected between gene nodes that were present in the microarray platform and directly linked on the pathway database used to evaluate the results .
indirect relationships were detected when the database indicated relationships between genes that were related through intermediate genes not present on the microarray platform.
the adequacy of the approach was demonstrated in the three independent data sets using a bilateral strategy.
first, our approach encompasses a wide range of network models, from no relationships between genes and single gaussian distributions to networks with variable number of relationships across sub-networks and mixtures of multiple gaussian distributions.
our approach identified the gene relationships and number of mixture components best supported by the data for each sub-network, accounting for model complexity and data available using the bic score.
multiple mixture components were favored over single components in the vast majority of the sub-networks inferred in the three data sets studied suggesting that adequacy of our approach to describe gene networks.
second, the inferred networks were compared to known network and the vast majority of the relationships were confirmed.
RESULTS
honey bee data set and circadian rhythm pathway
the mixture bayesian network approach provided the first published description of the relationship between genes in the circadian rhythm pathway of honey bees.
in addition, this approach offered novel information on general and condition-specific gene relationships.
figure  <dig> summarizes the circadian rhythm gene network predicted by the mixture bayesian network approach on each of the four bee-host sub-data sets used for cross-validation.
solid and dash-dotted arrows  denote direct and indirect relationships, respectively, confirmed in the circadian rhythm pathway for the fruit fly in the kegg database.
the number beside each edge denotes the fraction of data sets that supports that edge.
for example, 4/ <dig> denotes that the gene relationship was detected on all four sub-data sets.
all direct and indirect relationships detected by the mixture bayesian network on the honey bee genes were consistent with the relationships depicted for the fruit fly circadian rhythm pathway on the kegg database.
in addition, there were no predicted relationships that were not present on the kegg fruit fly pathway.
previous reports indicated that maturation age was the highest source of differentially expressed genes in the honey bee microarray study under consideration, followed by the subspecies of the bee sampled and the species of the host colony  <cit> .
in agreement with these studies, the mixture bayesian network consistently detected most edges of the circadian rhythm pathway across multiple bee-host sub-data sets .
however, not all relationships were detected in all four sub-data sets.
four edges were only detected on half of the sub-data sets , and one edge was only detected on one sub-data set .
the inability to consistently predict the same edges on all four sub-data sets may be due to two reasons.
first, previous work noted that the subspecies of the honey bee sampled was another source of differentially expressed genes, albeit less important than maturation  <cit> .
differences on edge prediction among sub-data sets may be due to differences in the honey bee subspecies used in each sub-data set.
second, the size of each bee-host sub-data set may have resulted in insufficient information to predict some edges in some sub-data sets.
to test this hypothesis, the mixture bayesian network approach was applied to the complete data set.
all the edges depicted in figure  <dig> were also predicted on the complete data set.
this result suggests that most or all data sets contained consistent information on the edges, although the information within data set may have not permitted an accurate prediction of some edges.
alternatively, these gene relationships were present on one or two sub-data sets, but the strength of the relationship on these sub-data sets was sufficient to be detected across all four sub-data sets.
all the edges in figure  <dig> were also inferred when the complete honey bee data set was analyzed and thus, results based on the complete honey bee data set are discussed.
the reliability of the gene network predicted by the mixture bayesian network approach was assessed using permutation-based resampling.
of the networks predicted from the randomized data sets, 15% and 24% did not share any or only shared one relationship with the network inferred from the original data and confirmed against the kegg database, respectively.
the maximum number of gene relationships in common between the randomized and original data sets was five.
only 1% and 7% of the networks predicted from the randomized data sets shared five and four relationships out of the  <dig> gene relationships with the network originally inferred, respectively.
similarly, 32% of the networks predicted from randomized data sets had more than eight false positive relationships and only 5% of networks had four or less false positives relationships.
the hegemony of these findings further advocates the capability of the mixture bayesian network approach to uncover evidence on gene relationships contained in the data and minimize the probability of networks predicted by chance alone.
the insights into condition-dependent gene relationships offered by the mixture bayesian network were evident on the sub-network of gene cyc, including parent genes per, dbt, and vri.
using the complete honey bee data set, the cyc sub-network was best described by a mixture of two components.
the parameter estimates corresponding to each mixture component of a gene sub-network offered information on the expression and co-expression patterns.
mixture components  <dig> and  <dig> had weights equal to  <dig>  and  <dig> , respectively.
the most notable differences between the two mixture components were the level of gene cyc  and the correlations between gene cyc and parent genes vri and sgg.
the correlations of cyc with both parent genes were negative in mixture component  <dig> and positive in mixture component  <dig> and the difference between the correlations of both mixture components was significant at p-value <  <dig>  and <  <dig>  × 10- <dig> for cyc and vri and cyc and sgg, respectively  <cit> .
bee-host subspecies combinations: ll, lm, ml, and mm.
maturation ages: h <dig>  h <dig>  h <dig>  h <dig>  h <dig>  f4
mix <dig> and mix2: mixture components  <dig> and  <dig>  respectively.
the mixture bayesian network approach was able to identify different models  of association between genes and of gene expression profiles depending on the conditions under consideration, and also provided a general model that is the weighted average of the condition-specific components.
cyran et al.
noted that the fruit fly circadian clock encompasses two interlocked transcriptional feedback loops  <cit> .
in one loop, dclock/cyc activates the expression of per, and the per protein in turn inhibits the activity of dclock/cyc.
in addition, vri and pdp <dig> encode related transcription factors whose expressions are directly activated by dclock/cyc.
in the second loop, the vri and pdp <dig> proteins feed back and directly regulate the expression of dclock.
results from the mixture bayesian network approach  confirm these findings.
the identification and interpretation of general and condition-dependent gene co-regulation patterns, using the odds that observations pertain to each mixture component, and the mixture component parameter estimates presented for the gene cyc sub-network can be repeated for all other genes in the network.
although the optimal number of mixture components varied across gene sub-networks , similar general and condition-dependent conclusions were obtained for other gene sub-networks in the circadian rhythm pathway and thus are not presented here.
this situation may be due to the high degree of interconnectedness among the gene nodes in this network.
for networks encompassing less related sub-networks, parameter estimates from each sub-network module should be evaluated.
rubin et al.
postulated, based on phylogenetic and correlation analyses of the function and domains of two proteins , that the circadian rhythm pathway of the honey bees resembles more the mammalian than the fruit fly counterpart  <cit> .
the similarity between the honey bee and fruit fly circadian rhythm networks found in the present study may be due to absence of reporters associated with genes that differ across insect species, such as the genes considered by rubin et al., on the microarray study under consideration  <cit> .
several indirect relationships among genes reported in the fruit fly kegg pathway were predicted as direct ones due to the absence of microarray reporters corresponding to the intermediate genes.
although the identification of relationships  among genes constitutes a major step towards a comprehensive description of biological pathways, precise characterization of the relationships is also important.
differences between the predicted and reported  gene relationships may be due to differences among the experiments considered  and insufficient or inaccurate information .
additional studies are necessary to ascertain the true nature  of these gene relationships.
mouse embryo data set and adherens junction pathway
yokoyama et al.
studied the relationship among various proteins in the adherens junction pathway  <cit> .
using the mixture bayesian network approach, the first depiction of general and condition-specific gene profiles and interactions in the adherens junction pathway of mouse embryos was obtained.
figure  <dig> depicts the adherens junction gene network predicted by the mixture bayesian network approach.
this approach predicted directional relationships in which the parent and child gene nodes are identified.
however, the adherens-junction mouse pathway in the kegg database included some non-directional gene relationships.
because the inferred network was validated against the kegg pathway, non-directional gene relationships are denoted with lines instead of arrows.
all direct and indirect relationships between genes predicted by the mixture bayesian network approach were present in the kegg database.
consistency between the predicted and reported adherens junction networks was also observed for two additional pathways  predicted using the same mouse microarray data set  <cit> .
the parameter estimates and odds that observations pertained to different mixture model components enhanced the understanding of the differences in gene profiles and relationships between mouse embryo treatments.
for example, a mixture of three components was best suited to describe the sub-network of gene α-actinin.
mixture components one, two and three had weights  <dig> ,  <dig>  and  <dig>  and modeled low, intermediate and high levels of expression of all genes in the α-actinin sub-network, respectively.
in addition, gene α-actinin had a negative correlation with parent gene actin in the first mixture component meanwhile this correlation was positive in the other two mixture components.
the statistical significance of the difference between correlations from mixture components one and two, and one and three were p-value <  <dig>  × 10- <dig> and p-value <  <dig>  × 10- <dig>  respectively meanwhile the correlations from mixture components two and three were non-significantly different .
studies: gene expression omnibus series  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>
treatments: control or cnt, methylmercury or mhg, ethanol or eth, 2-chloro-2'deoxyadenosine or cda, oxygen or oxg, pk <dig> or pk, and ro5- <dig> or ro.
mix <dig> and mix2: mixture components  <dig> and  <dig> respectively.
genes in sub-network: actin, rhoa, afadin, α-catenin, vinculin, α-actinin, and zo- <dig>
consideration of the assignment of treatment observations, regardless of study, to the mixture components for the actin sub-network revealed new and confirmed known biological interpretations .
for example, the only three conditions receiving 2-chloro-2'deoxyadenosine or cda that were assigned to the mixture component mix <dig> correspond to sampling intervals before  and after  the critical event of p <dig> protein induction expressed in cda-treated embryos between  <dig> - and  <dig>  h post-exposure  and to embryos also treated with ro5- <dig> , a presumed agonist that is weakly teratogenic .
all the observations during the critical event  were assigned to mix <dig>  these results suggested that sampling time and treatments in addition to cda could have substantial impact on the actin sub-network.
the odds that observations pertain to each mixture component also provided insights into the effects of two pharmacologic agents, pk <dig>  and ro5- <dig> .
the pk agent can rescue embryos from the effects of toxic agents meanwhile ro is a presumed agonist that can have weak teratogenic effects.
the fairly even assignment of pk and ro treatments to both mixture components suggested that the impact of these treatments on the actin sub-network was dependent on the other conditions of the experiments.
similarly, the evaluation of the assignment of control  samples to mixture components offered information on the behavior of seemingly similar samples.
out of ten control samples , seven were better described using mix <dig>  among the three control samples assigned to mix <dig>  one corresponded to the mouse sub-strain c57bl/6j  and was consistent with the assignment of all other b6j samples that received treatment to mix <dig>  this finding suggested the potentially high effect of the genotype of the embryo on the relationship among genes in the actin sub-network.
the other two control samples assigned to mix <dig> corresponded to series gse  <dig> and  <dig>  and these assignments were consistent with the assignment of all  or the vast majority  of the samples to mix <dig>  these assignments are consistent with the strong series effect previously discussed.
yeast synchronization data set and cell cycle pathway
of the  <dig> gene relationships predicted by the mixture bayesian approach, nine were confirmed in friedman et al., seven and three were listed in the biogrid and kegg databases respectively, and three and one were reported in the literature  <cit>  and sgd database, respectively.
several direct relationsh ips predicted by the mixture bayesian network approach were not reported in the biogrid, kegg, literature, sgd databases or in friedman et al.
<cit> .
likewise, only three out of the  <dig> directed relationships and four out of the  <dig> genes in the pathway predicted by friedman et al.
are present in the kegg pathway  <cit> .
although experimental studies are necessary to confirm or disprove these findings, additional gene relationships and literature review imply that some are likely to be true  relationships.
for example, through text mining of biomedical literature, liu et al.
found that genes mnn <dig>  cln <dig>  and cln <dig> were all influenced by the same activators  <cit> .
this finding supports the novel relationships between mnn <dig> and either cln <dig> or cln <dig> detected by the mixture bayesian network approach.
had microarray data on the common activators been available, the relationship between these genes could have been further resolved.
additional studies are necessary to determine if the unconfirmed relationships are true positives instead of false positives.
the mixture bayesian network approach was not able to detect relationships between genes svs <dig> and alk <dig> and the rest of the genes in the network.
these results are consistent with the absence of evidence relating genes svs <dig> and alk <dig> in the literature or databases.
assessment of the mixture model parameter estimates was used to gain insights into the gene sub-networks.
for example, a mixture model with two components receiving equal weight was found to describe the sub-network of gene cdc <dig> depicted in figure  <dig>  the mean expression of the genes in mixture component one was between a half and a third of the mean expression of the genes in mixture component two.
in addition, genes cdc <dig> and clb <dig> had a positive correlation  in mixture component two, meanwhile the correlation was weaker and negative  in mixture component one and these correlations were different at p-value <  <dig> ).
the assignment of the samples to the components of the mixture model describing gene sub-networks helped in the understanding of the changes in the network across the different treatments or conditions studied.
for example, consideration of the probability that a sample pertained to either one of the two components of the mixture model used to describe the sub-network of gene cdc <dig> confirmed some expected assignments and suggested unexpected commonalities between samples.
the assignment of elutrition samples collected every  <dig> minutes to mixture components provided a good example of a slow transition of the network behavior across time.
the first six samples collected between  <dig> and  <dig> minutes were assigned to one mixture component and six samples obtained from  <dig> to  <dig> minutes were assigned to the other mixture component.
samples collected at intermediate time points  had a similar odds ratio to pertain to either mixture component.
a discussion of the assignment of samples in the α-factor, cdc <dig> and cdc <dig> synchronization is provided in additional file  <dig>
overall, the mixture bayesian network approach detected more confirmed relationships than the approach presented by friedman et al.
<cit> .
however, friedman et al.
detected five relationships  that were not predicted by the mixture bayesian approach and were not present in the biogrid database.
the higher number of relationships predicted by the mixture bayesian network approach relative to the bayesian network implementations of friedman et al.
may be due to the flexibility of the mixture model in the mixture bayesian network approach to accommodate fluctuations on the gene patterns and relationships across the multiple conditions evaluated  <cit> .
single, fixed-number, and optimal mixture densities
the mixture bayesian network approach to infer gene networks implemented in this study used the em algorithm to obtain parameter estimates and the bic score to identify the optimal number of mixture components supported by the data.
this framework offered great flexibility to accommodate a wide range of models and probability density functions across gene sub-networks and within gene sub-network, across the conditions under study.
a series of benchmarking tests were conducted to evaluate the benefits of the more flexible and data-driven mixture bayesian network approach for different scenarios.
the description of the sub-networks supported by the data based on bic score favored mixture over single distributions.
this indicates that for the networks and data sets considered in our study, mixtures offer a better description of the relationship between genes than single distributions.
this result was further validated when the assignment of samples from different treatments or conditions to mixture components was investigated.
the robustness of the mixture bayesian network learning process granted by the use of mixture distributions and potentially variable numbers of mixture components across the network was investigated.
figure  <dig> presents the bic scores of the networks predicted using the mixture bayesian network approach with the number of mixture components fixed to single or mixture of two or three gaussian densities , or estimated from the data  for the circadian rhythm and adherens junction networks previously presented.
the bic score of the optimal mixture implementation was higher than the approaches that used a fixed number of mixtures across all subnetworks.
the superiority of the mixture bayesian network approach that estimates the optimal number of mixture components for each gene sub-network from the data over approaches with fixed number of components was evident in both networks.
in addition, the maximum number of parents in a gene sub-network was set to five in this study, based on a review of pathways reported in the kegg database and associated to the pathways studied.
this maximum was supported by the fact that the optimal bic score was always found with one to four parent genes.
computational evaluation
the computational time required for the prediction of a network across a wide range of number of genes was evaluated.
the mixture bayesian network approach was implemented using the c++ programming language.
figure  <dig> shows the computational time  required to reconstruct networks of different size, ranging from five to  <dig> genes sampled at random from the circadian rhythm  and adherens junction  microarray data sets.
five data sets were obtained from the complete microarray data sets for each network size considered, and a network was estimated for each data set and network size.
the computing time to predict a network depends in part on the degree of connectivity between the genes studied.
because the gene sets were drawn at random from the microarray data sets, the degree of relationship between these genes  was expected to vary across data sets within network size and across network sizes studied.
thus, both the number of edges predicted in each network and the computational time required to predict the network were recorded.
figure  <dig> and figure  <dig> present the average computational time and number of edges across the five data sets within network size evaluated for two cases, when the number of mixture components was estimated from the data  or fixed at  <dig> components .
these comparisons allowed the assessment of the impact of the mixture component of the algorithm on the computational time and inferred edges.
as expected, the higher the number of genes in the network, the higher the computational time required for estimation of the complete network .
the lower computational time to infer the adherens junction network compared to the circadian rhythm network may be due to the fewer relationships  detected in the first network .
in the adherens junction network, the trends of computational time and number of edges across number of genes were similar.
in the circadian rhythm network, the computing time increased at a slightly or moderately higher rate than the number of edges across the number of genes.
the disparity between the trends of time and number of edges was higher in the adherens junction network with a fixed number of mixture components than when the number of components was estimated from the data.
this result may be due to the potential inadequacy of the fixed number of mixture components to describe some or all gene sub-networks.
thus, more iterations and more time were required to estimate the model parameters and infer the network structure.
furthermore, the nearly constant number of edges  inferred in the adherens junction across the number of gene nodes studied  suggested that the mixture bayesian network algorithm was not biased towards inferring higher number of edges with higher number of gene nodes network .
algorithmic extensions
the identification of a bayesian network structure is an np-hard problem because the number of possible candidate structures super-exponentially increases with the number of gene nodes.
two structure learning methods that have complementary advantages were evaluated in this study, the sparse candidate algorithm  <cit>  and the simulated annealing method  <cit> .
the sparse candidate algorithm achieves fast structure learning by restricting the search space of potential parent gene nodes.
intuitively we would like to restrict the search space as less as possible.
in this study, the algorithm was modified to minimize the impact of restricting the search space, while also maintaining efficiency by evaluating additional parent genes for each gene sub-network with less than the maximum number of parent genes.
the network structure was also inferred using a simulated annealing algorithm to overcome the potential risk of inferring networks at local maxima instead of global maximum  <cit> .
at each iteration, the simulated annealing algorithm randomly jumped from the current network topology to a nearby topology formed from the addition, deletion, or reversion of an edge in the current topology.
simulated annealing can escape from potential local minima by occasionally sampling from less suitable networks and the probability of jumping to another status is empirically adjusted.
this procedure is repeated until convergence.
the implementation of these two algorithms allowed the evaluation of the mixture bayesian network approach to be independent of the actual algorithm used to learn the topology of the network.
the vast majority of the gene relationships was predicted by both the sparse candidate and simulated annealing algorithms.
therefore, only results from the sparse candidate algorithm are presented here.
the mixture bayesian approach presented here can be extended to include other gene  selection criteria, consider other distributional assumptions and prior information on the relationship between genes in a network.
these extensions are discussed in additional file  <dig>
CONCLUSIONS
a comprehensive characterization of the circadian rhythm pathway in honey bees and of the adherens junction pathway in mouse embryos was obtained using a mixture bayesian network approach.
all the gene relationships detected by the mixture bayesian network based on honey bee gene expression information were consistent with known relationships in the fruit fly circadian rhythm pathway.
the mixture bayesian network approach was also able to uncover changes in gene relationships and profiles in the cyc gene sub-network associated with changes on the maturation age, subspecies of honey bee sampled, and subspecies of honey bees in the host hive.
likewise, all the gene relationships detected by the mixture bayesian network using mouse-embryo gene expression information were consistent with known relationships in the general mouse adherens junction pathway.
all or the vast majority of the samples within experiment were assigned to the same mixture component model in seven out of the nine experiments considered.
this assignment of samples to mixture components suggested that experiment had a major impact on the gene relationships and expression profiles in the actin gene sub-network.
the application of the mixture bayesian network approach to the well-studied yeast cell cycle confirmed many reported relationships.
in addition, the assignment of samples to the components of the mixture model for the cdc <dig> gene sub-network provided novel insights into gene associations.
for example, the assignment of elutrition samples to the mixture components suggested a slow transition of the network behavior across time.
on the other hand, the assignment of cdc15-synchronized samples to mixture components suggested that the networks were instable across most of the time-points considered.
this study demonstrated that the mixture bayesian network approach is well-suited to infer gene pathways based on microarray gene expression data.
the estimation of the number of mixture components from the data allowed the simultaneous modeling of gene sub-networks with either constant or variable behavior across conditions.
this modeling approach has two potential advantages, a) it can increase the accuracy of the inferred relationships and precision of the parameter estimates, and b) it can help identify groups or clusters of observations that support particular  models or mixture components.
the superior accuracy of the mixture bayesian network approach was demonstrated by the better bic score of the network with the number of mixture components estimated from the data, compared to the bic score of approaches with a fixed number of mixture components.
the mixture of gaussian distributions demonstrated in this study naturally accommodated the continuous and typically normally-distributed gene expression measurements.
the proposed approach can be extended to use mixtures of other discrete or continuous distributions, or incorporate prior information on the genes in the network.
the proposed mixture bayesian network approach allows the detection of general and condition-specific relationships among genes, and the prediction of an overall network structure within the framework of probabilistic inference.
the interpretation of parameter estimates and functions  can offer a more comprehensive understanding of gene patterns and relationships.
insights into the gene patterns, and the relationships between them, can be used to design effective follow-up experiments and investigate the novel relationships predicted by the mixture bayesian network approach.
authors' contributions
yk implemented the mixture bayesian network algorithm, inferred the three gene networks, contributed to the interpretation of results, and drafted the manuscript.
c-xz helped interpret the algorithm and reviewed the manuscript.
srz obtained funding for the study, participated in its conception, coordination, interpretation of results, and helped write the manuscript.
all authors have read and approved the final version of this manuscript.
supplementary material
additional file 1
representation of a mixture model for the bayesian sub-network corresponding to child gene genej and parent genes: gene <dig>  gene <dig> ..., genep- <dig>  the figures represent the bayesian sub-network of the jth gene with p- <dig> parent genes.
click here for file
additional file 2
illustration of the network-structure learning process.
diagram of the steps of the algorithm to infer the topology of the network.
click here for file
additional file 3
estimation of mixture parameters.
formulae to estimate the parameters using the expectation-maximization algorithm.
click here for file
additional file 4
assignment of yeast cell cycle conditions to mixture components.
distribution of yeast cell cycle treatments across the mixture components of the gene sub-network.
click here for file
additional file 5
additional algorithmic extensions.
potential modifications of the proposed algorithm are discussed.
click here for file
acknowledgements
we would like to thank heather adams, jaebum kim and bruce southey for useful comments.
the support of nih/nigms , nsf/itr  and, nih nida  are greatly appreciated.

 relationships between genes and gene pathways can be inferred based on gene expression profiles across conditions, treatments, or samples obtained from microarray experiments. results from these studies can aid in the confirmation of previously known pathways and motivate the study of newly uncovered relationships among genes. a bayesian network approach is well-suited to detect relationships between genes using an acyclic direct graph. needham et al. provided an in depth primer on learning bayesian networks for computational biology. this approach has a solid theoretical foundation, offers a probabilistic framework to describe the variation typically observed in microarray data, accommodates missing data, and incorporates prior knowledge on gene relationships. several implementations of bayesian networks to infer gene networks have been reported and two potential weaknesses have been identified. first, some applications require the transformation of continuous gene expression data into discrete input data, and this can influence the resulting network or resulting in loss of information. thus, analysis of data that has not been discretized is favored. second, most bayesian network implementations assume standard binomial, multinomial probability or single gaussian probability density functions of gene expression across a wide range of conditions. these distributions may fail to accommodate multimodal or skewed distributions associated with condition-dependent networks that exhibit changes in gene expression or gene relationships across conditions. this is because gene network inference is typically based on gene expression measurements across a wide range of conditions. mixture models can be used to address this limitation because these models can describe potentially complex distributions of gene expression across a wide range of conditions. newman and leicht explored the capability of probabilistic mixture models to detect a very broad range of network structures without prior knowledge. ko et al. presented preliminary results from the application of bayesian mixture algorithm to infer pathways. they presented a histone pathway, did not use parameter estimates to identify condition-dependent gene relationships, and did not use cross-validation to assess the adequacy of the inferred network. in this study, gene networks were inferred using a mixture bayesian network approach. the superiority of this approach was evidenced on the new insights into general and condition-dependent relationships between genes and gene expression profiles gained while addressing some of the limitations encountered in previous applications of bayesian network to infer gene networks. we demonstrate that the integration of mixture models and bayesian network is well-suited to infer the structure of gene networks from continuous gene expression measurements across a wide range of conditions. the estimated parameters were easy to interpret, and aided in the characterization of gene expression and co-expression profiles. furthermore, functions of these parameter estimates allowed the identification of condition-dependent networks. the mixture bayesian network approach was used to infer three pathways on three independent data sets. the first prediction of the circadian rhythm pathway in honey bees and the first prediction of a cell communication pathway in mouse embryos were obtained. in addition, the mixture bayesian network approach was applied to a yeast cell-cycle data set commonly used in many gene network studies. comparisons of the predicted and known pathways and benchmark tests confirmed the outstanding performance of the mixture bayesian network approach. novel insights into general and condition-specific gene relationships and expression patterns were obtained, thus demonstrating the strength of the proposed approach to infer gene networks. methods mixture bayesian network model bayesian networks can be described as a directed acyclic graph with nodes representing random variables or genes and directed edges representing the relationships between the nodes. given a set of genes {g g..., gn}, the corresponding bayesian network g is represented as the joint probability distribution over all genes in g or p. applying the narrow sense form of the markov property for a stochastic process, each gene is independent of non-descendant genes in the network, given the parent genes or nodes. this conditional independence property allows the factorization of the joint probability distribution as the product of conditional probabilities. the overall network is conceived as a set of gene subnetworks, each corresponding to a given gene node gj and its associated parent gene nodes a for j = to n. models that accommodate for potential changes in gene expression and co-expression patterns across conditions are necessary when gene expression data across multiple conditions is used to infer gene networks. mixture models are a flexible and effective option to describe changes in gene co-expression patterns across conditions because the joint probability functions of sub-networks are modeled with a combination of mixture components. the concept of integrating mixtures of gaussian densities into bayesian networks was introduced by davies and moore, and applied to word and social networks by newman and leicht. under the bayesian network framework, where n is the total number of genes in the network and a is the set of parent genes of child gene gj in the j-th sub-network. the conditional probability density function of the j-th gene p), given the set of parent genes is expressed as the ratio between the joint probability density function of the parent and child genes in the j-th sub-network )), and the marginal probability density function of the parent genes in the j-th sub-network )). the joint and marginal probability density functions are assumed to follow a multivariate gaussian distribution and a mixture model is used to describe possible sub-network structures. the probability density function for the j-th sub-network is represented with a mixture of kj multivariate gaussian distributions each with a weight αjk. here, xji denotes a pj dimensional vector including the i-th expression observation of the child and parent genes in the j-th sub-network. for example, if gj has three parents then pj = therefore, the joint probability density function for the j-th sub-network is and where. for the joint probability density function, each component k of the mixture is described with a mean vector μjk of dimension pj and a variance-covariance matrix Σjk of dimension pj × pj. the marginal probability density function of the parent genes a is described in the same way as the joint probability density function but the dimensions of the i-th vector of observations, the mean vector and the variance-covariance matrix are,, and ×, respectively. therefore, the overall likelihood of the data across all genes in the network g is: an example of a bayesian gene sub-network is provided in additional file the relationship between the child and parent genes in each sub-network can be described by the weighted sum of the correlations between the parent geneand child gene estimated for each mixture component. the weights correspond to the mixture component weights. this computation was described by bland and altman and is presented in additional file overall network learning algorithm the learning algorithm integrated two major components, a) the inference of each individual gene sub-network structure and b) the inference of the overall network structure. the first component consisted of the identification of the sub-network for each gene, the estimation of the parameters of the mixture models using the expectation-maximization algorithm, and the identification of the number of components of the mixture model best supported by the data. the second component consisted of the combination of the individual gene sub-networks into an overall network, and the removal of cyclic relationships between genes in the overall network. a detailed description of the network structure learning approach is presented in additional file the equations in the expectation-maximization steps corresponding to a gene sub-network follow the derivations of bilmes and are summarized in additional file the bayesian information criterion was used to identify the mixture model best supported by the data for each gene sub-network, and to evaluate the overall network. bayesian information criterion offers a good compromise between model adequacy and parsimony, minimizing the possibility of over-fitting the data with highly parameterized networks. this feature of the bic is advantageous in gene network inference because of the large number of relationships and conditions that could potentially be considered. newman and leicht used probabilistic mixture models to infer a wide range of network structures and noted that the convergence of the algorithm to the global maxima is not guaranteed. to address this potential problem, we followed the approach implemented by newman and leicht, and networks were inferred using a range of starting points for each mixture parameter. each network presented in this study had the highest likelihood over all sets of initial conditions or runs evaluated. furthermore, davies and moore demonstrated the suitability of mixtures of gaussian densities estimated with the em algorithm to infer bayesian networks using data from an extensive astronomical survey and a high-throughput biological cell assay. extensive evaluation of the performance of the mixture bayesian network approach to infer networks was presented by newman and leicht and davies and moore. in this study, we concentrated on the application of this methodology to infer gene networks using microarray data, and demonstrated how the interpretation of parameter estimates and function thereof offered new insights into the relationships between genes. a unique advantage of the integration of the mixture models and bayesian network is the ability to uncover general and condition-specific gene networks. the parameter estimates of each mixture component of a gene sub-network were used as indicators of the gene expression level, variation, and co-variation between genes. in addition, the mixture parameter estimates were used to compute the odds that the observation corresponded to a particular mixture component relative to the other components. each sample had a non-zero probability to pertain to each of the mixture components. for ease of interpretation and to overcome the limitation of using a single probability threshold across sub-networks with different number of mixture components, no threshold was used to assign samples to mixtures. instead, the samples were assigned to the mixture component that had the highest odds. the assignment of samples to the most likely mixture component allowed the identification of samples that shared the same expression and co-expression patterns for a gene sub-network. data gene networks corresponding to three gene pathways were inferred from three microarray gene expression data sets using a mixture bayesian network approach. the goal of these applications was to gain insights into general and condition-dependent relationships between genes and gene expression profiles, and demonstrate the flexibility and adequacy of the mixture bayesian network approach to characterize networks. honey bee data set and circadian rhythm pathway a detailed description of the honey bee microarray data set is presented by whitfield et al. and rodriguez-zas et al.. briefly, gene expressions during behavioral maturation in the brains of honey bees from two subspecies rose in one of two host colonies representing two different environments were measured. within combinations of bee and host race, nurse bees were sampled on days and after adult emergence and forager bees were sampled at day after emergence. the experimental design allowed the division of the honey bee data set into four independent bee-host subspecies-combination sub-data set. thus, the gene network was inferred on the complete data set and also on each of the four sub-data set. processing of the data from spotted cdna microarrays included background subtraction and log2-transformation of fluorescence intensities, and removal of flagged spots or spots that did not surpass a minimum intensity threshold of the log intensity values were normalized using a loess transformation, and global dye and microarray effects were removed using the approach implemented beehive. six genes assigned to the fruit fly circadian rhythm pathway in the kegg database were present on the microarray platform. the kegg identifier for the circadian rhythm pathway in the fruit fly is dme. mouse embryo data set and adherens junction pathway data from nine microarray experiments that evaluated the effect of multiple toxic agents on mouse embryo gene expression levels were used to reconstruct the adherens junction pathway. a detailed description of this data set is provided by rodriguez-zas et al.. briefly, the experimental design consisted of two types of direct comparisons with reverse labeling; comparisons between samples that received different treatments, and comparisons between treated and control samples. a total of microarrays from the same platform were available for analysis, and data processing and normalization followed the same protocols applied to the honey bee data set and implemented in beehive. seven genes assigned to the mouse adherens junction pathway in the kegg database were present on the microarray platform. the kegg identifier for the adherens junction pathway in the mouse is mmu. yeast synchronization data set and cell cycle pathway the mixture bayesian network approach was used to predict the yeast cell-cycle network using gene expression data reported by. the data set consisted of gene expression measurements from microarrays across six experiments. the yeast experiments evaluated yeast cell synchronization via the arrest of cells by one of four conditions, followed by activation or release of the cells from the arresting condition. the synchronization conditions were cdc15-, cdc28-, α factor- and size or elutrition- based synchronization, and samples were obtained at various time points after the removal of the last four arresting conditions. in a separate study, the effects of inducing either the g cyclin cln3p or the b-type cyclin clb2p were examined, each in two experiments. the normalized gene expression data is available at the stanford microarray database and the experiment identifiers are: to and a network including genes analyzed by friedman et al. using multinomial and single gaussian distribution was inferred and results were compared. the kegg identifier for the yeast cell cycle pathway is sce. validation of the inferred gene networks the relationships between genes inferred by the mixture bayesian network approach were compared to relationships reported in the kegg database, biogrid database, and saccharomyces genome database, or sgd, and in the literatures. these databases were accessed on august and provide complementary information. the kegg database includes cause-effect or directional and non-directional relationships between genes depicted in pathway diagrams. in directional relationships, the child and parent genes are identified, meanwhile in non-directional relationships, genes are linked but there is neither a parent nor a child gene. the biogrid database lists all known interactions or undirected relationships between genes. due to differences in the curation and annotation protocols of each database, biogrid can encompass more relationships than kegg, but these relationships do not have associated direction. two strategies were used to confirm the reliability of the gene networks predicted by the mixture bayesian network approach. first, a k-fold cross-validation approach was implemented. the design of the circadian rhythm experiment allowed to infer the gene network in each of four independent data sets corresponding to the four bee-host race combinations that encompassed the same five behavioral maturation stages or ages. assessment of the mixture bayesian network approach consisted on the comparison of the four inferred networks. second, a permutation-based resampling approach was used to generate randomized data sets with the same structure as the original circadian rhythm data set. these data sets were analyzed using the mixture bayesian network approach and the predicted networks were compared to the observed network and the known pathway in the database. not all genes present in the pathway databases had reporters in the microarray platforms studied and thus two types of gene relationships, direct and indirect were inferred. direct relationships were detected between gene nodes that were present in the microarray platform and directly linked on the pathway database used to evaluate the results. indirect relationships were detected when the database indicated relationships between genes that were related through intermediate genes not present on the microarray platform. the adequacy of the approach was demonstrated in the three independent data sets using a bilateral strategy. first, our approach encompasses a wide range of network models, from no relationships between genes and single gaussian distributions to networks with variable number of relationships across sub-networks and mixtures of multiple gaussian distributions. our approach identified the gene relationships and number of mixture components best supported by the data for each sub-network, accounting for model complexity and data available using the bic score. multiple mixture components were favored over single components in the vast majority of the sub-networks inferred in the three data sets studied suggesting that adequacy of our approach to describe gene networks. second, the inferred networks were compared to known network and the vast majority of the relationships were confirmed.  honey bee data set and circadian rhythm pathway the mixture bayesian network approach provided the first published description of the relationship between genes in the circadian rhythm pathway of honey bees. in addition, this approach offered novel information on general and condition-specific gene relationships. figure summarizes the circadian rhythm gene network predicted by the mixture bayesian network approach on each of the four bee-host sub-data sets used for cross-validation. solid and dash-dotted arrows denote direct and indirect relationships, respectively, confirmed in the circadian rhythm pathway for the fruit fly in the kegg database. the number beside each edge denotes the fraction of data sets that supports that edge. for example, 4/ denotes that the gene relationship was detected on all four sub-data sets. all direct and indirect relationships detected by the mixture bayesian network on the honey bee genes were consistent with the relationships depicted for the fruit fly circadian rhythm pathway on the kegg database. in addition, there were no predicted relationships that were not present on the kegg fruit fly pathway. previous reports indicated that maturation age was the highest source of differentially expressed genes in the honey bee microarray study under consideration, followed by the subspecies of the bee sampled and the species of the host colony. in agreement with these studies, the mixture bayesian network consistently detected most edges of the circadian rhythm pathway across multiple bee-host sub-data sets. however, not all relationships were detected in all four sub-data sets. four edges were only detected on half of the sub-data sets, and one edge was only detected on one sub-data set. the inability to consistently predict the same edges on all four sub-data sets may be due to two reasons. first, previous work noted that the subspecies of the honey bee sampled was another source of differentially expressed genes, albeit less important than maturation. differences on edge prediction among sub-data sets may be due to differences in the honey bee subspecies used in each sub-data set. second, the size of each bee-host sub-data set may have resulted in insufficient information to predict some edges in some sub-data sets. to test this hypothesis, the mixture bayesian network approach was applied to the complete data set. all the edges depicted in figure were also predicted on the complete data set. this result suggests that most or all data sets contained consistent information on the edges, although the information within data set may have not permitted an accurate prediction of some edges. alternatively, these gene relationships were present on one or two sub-data sets, but the strength of the relationship on these sub-data sets was sufficient to be detected across all four sub-data sets. all the edges in figure were also inferred when the complete honey bee data set was analyzed and thus, results based on the complete honey bee data set are discussed. the reliability of the gene network predicted by the mixture bayesian network approach was assessed using permutation-based resampling. of the networks predicted from the randomized data sets, 15% and 24% did not share any or only shared one relationship with the network inferred from the original data and confirmed against the kegg database, respectively. the maximum number of gene relationships in common between the randomized and original data sets was five. only 1% and 7% of the networks predicted from the randomized data sets shared five and four relationships out of the gene relationships with the network originally inferred, respectively. similarly, 32% of the networks predicted from randomized data sets had more than eight false positive relationships and only 5% of networks had four or less false positives relationships. the hegemony of these findings further advocates the capability of the mixture bayesian network approach to uncover evidence on gene relationships contained in the data and minimize the probability of networks predicted by chance alone. the insights into condition-dependent gene relationships offered by the mixture bayesian network were evident on the sub-network of gene cyc, including parent genes per, dbt, and vri. using the complete honey bee data set, the cyc sub-network was best described by a mixture of two components. the parameter estimates corresponding to each mixture component of a gene sub-network offered information on the expression and co-expression patterns. mixture components and had weights equal to and, respectively. the most notable differences between the two mixture components were the level of gene cyc and the correlations between gene cyc and parent genes vri and sgg. the correlations of cyc with both parent genes were negative in mixture component and positive in mixture component and the difference between the correlations of both mixture components was significant at p-value < and < × 10- for cyc and vri and cyc and sgg, respectively. bee-host subspecies combinations: ll, lm, ml, and mm. maturation ages: h h h h h f4 mix and mix2: mixture components and respectively. the mixture bayesian network approach was able to identify different models of association between genes and of gene expression profiles depending on the conditions under consideration, and also provided a general model that is the weighted average of the condition-specific components. cyran et al. noted that the fruit fly circadian clock encompasses two interlocked transcriptional feedback loops. in one loop, dclock/cyc activates the expression of per, and the per protein in turn inhibits the activity of dclock/cyc. in addition, vri and pdp encode related transcription factors whose expressions are directly activated by dclock/cyc. in the second loop, the vri and pdp proteins feed back and directly regulate the expression of dclock. results from the mixture bayesian network approach confirm these findings. the identification and interpretation of general and condition-dependent gene co-regulation patterns, using the odds that observations pertain to each mixture component, and the mixture component parameter estimates presented for the gene cyc sub-network can be repeated for all other genes in the network. although the optimal number of mixture components varied across gene sub-networks, similar general and condition-dependent conclusions were obtained for other gene sub-networks in the circadian rhythm pathway and thus are not presented here. this situation may be due to the high degree of interconnectedness among the gene nodes in this network. for networks encompassing less related sub-networks, parameter estimates from each sub-network module should be evaluated. rubin et al. postulated, based on phylogenetic and correlation analyses of the function and domains of two proteins, that the circadian rhythm pathway of the honey bees resembles more the mammalian than the fruit fly counterpart. the similarity between the honey bee and fruit fly circadian rhythm networks found in the present study may be due to absence of reporters associated with genes that differ across insect species, such as the genes considered by rubin et al., on the microarray study under consideration. several indirect relationships among genes reported in the fruit fly kegg pathway were predicted as direct ones due to the absence of microarray reporters corresponding to the intermediate genes. although the identification of relationships among genes constitutes a major step towards a comprehensive description of biological pathways, precise characterization of the relationships is also important. differences between the predicted and reported gene relationships may be due to differences among the experiments considered and insufficient or inaccurate information. additional studies are necessary to ascertain the true nature of these gene relationships. mouse embryo data set and adherens junction pathway yokoyama et al. studied the relationship among various proteins in the adherens junction pathway. using the mixture bayesian network approach, the first depiction of general and condition-specific gene profiles and interactions in the adherens junction pathway of mouse embryos was obtained. figure depicts the adherens junction gene network predicted by the mixture bayesian network approach. this approach predicted directional relationships in which the parent and child gene nodes are identified. however, the adherens-junction mouse pathway in the kegg database included some non-directional gene relationships. because the inferred network was validated against the kegg pathway, non-directional gene relationships are denoted with lines instead of arrows. all direct and indirect relationships between genes predicted by the mixture bayesian network approach were present in the kegg database. consistency between the predicted and reported adherens junction networks was also observed for two additional pathways predicted using the same mouse microarray data set. the parameter estimates and odds that observations pertained to different mixture model components enhanced the understanding of the differences in gene profiles and relationships between mouse embryo treatments. for example, a mixture of three components was best suited to describe the sub-network of gene α-actinin. mixture components one, two and three had weights, and and modeled low, intermediate and high levels of expression of all genes in the α-actinin sub-network, respectively. in addition, gene α-actinin had a negative correlation with parent gene actin in the first mixture component meanwhile this correlation was positive in the other two mixture components. the statistical significance of the difference between correlations from mixture components one and two, and one and three were p-value < × 10- and p-value < × 10- respectively meanwhile the correlations from mixture components two and three were non-significantly different. studies: gene expression omnibus series treatments: control or cnt, methylmercury or mhg, ethanol or eth, 2-chloro-2'deoxyadenosine or cda, oxygen or oxg, pk or pk, and ro5- or ro. mix and mix2: mixture components and respectively. genes in sub-network: actin, rhoa, afadin, α-catenin, vinculin, α-actinin, and zo- consideration of the assignment of treatment observations, regardless of study, to the mixture components for the actin sub-network revealed new and confirmed known biological interpretations. for example, the only three conditions receiving 2-chloro-2'deoxyadenosine or cda that were assigned to the mixture component mix correspond to sampling intervals before and after the critical event of p protein induction expressed in cda-treated embryos between - and h post-exposure and to embryos also treated with ro5-, a presumed agonist that is weakly teratogenic. all the observations during the critical event were assigned to mix these results suggested that sampling time and treatments in addition to cda could have substantial impact on the actin sub-network. the odds that observations pertain to each mixture component also provided insights into the effects of two pharmacologic agents, pk and ro5-. the pk agent can rescue embryos from the effects of toxic agents meanwhile ro is a presumed agonist that can have weak teratogenic effects. the fairly even assignment of pk and ro treatments to both mixture components suggested that the impact of these treatments on the actin sub-network was dependent on the other conditions of the experiments. similarly, the evaluation of the assignment of control samples to mixture components offered information on the behavior of seemingly similar samples. out of ten control samples, seven were better described using mix among the three control samples assigned to mix one corresponded to the mouse sub-strain c57bl/6j and was consistent with the assignment of all other b6j samples that received treatment to mix this finding suggested the potentially high effect of the genotype of the embryo on the relationship among genes in the actin sub-network. the other two control samples assigned to mix corresponded to series gse and and these assignments were consistent with the assignment of all or the vast majority of the samples to mix these assignments are consistent with the strong series effect previously discussed. yeast synchronization data set and cell cycle pathway of the gene relationships predicted by the mixture bayesian approach, nine were confirmed in friedman et al., seven and three were listed in the biogrid and kegg databases respectively, and three and one were reported in the literature and sgd database, respectively. several direct relationsh ips predicted by the mixture bayesian network approach were not reported in the biogrid, kegg, literature, sgd databases or in friedman et al.. likewise, only three out of the directed relationships and four out of the genes in the pathway predicted by friedman et al. are present in the kegg pathway. although experimental studies are necessary to confirm or disprove these findings, additional gene relationships and literature review imply that some are likely to be true relationships. for example, through text mining of biomedical literature, liu et al. found that genes mnn cln and cln were all influenced by the same activators. this finding supports the novel relationships between mnn and either cln or cln detected by the mixture bayesian network approach. had microarray data on the common activators been available, the relationship between these genes could have been further resolved. additional studies are necessary to determine if the unconfirmed relationships are true positives instead of false positives. the mixture bayesian network approach was not able to detect relationships between genes svs and alk and the rest of the genes in the network. these results are consistent with the absence of evidence relating genes svs and alk in the literature or databases. assessment of the mixture model parameter estimates was used to gain insights into the gene sub-networks. for example, a mixture model with two components receiving equal weight was found to describe the sub-network of gene cdc depicted in figure the mean expression of the genes in mixture component one was between a half and a third of the mean expression of the genes in mixture component two. in addition, genes cdc and clb had a positive correlation in mixture component two, meanwhile the correlation was weaker and negative in mixture component one and these correlations were different at p-value < ). the assignment of the samples to the components of the mixture model describing gene sub-networks helped in the understanding of the changes in the network across the different treatments or conditions studied. for example, consideration of the probability that a sample pertained to either one of the two components of the mixture model used to describe the sub-network of gene cdc confirmed some expected assignments and suggested unexpected commonalities between samples. the assignment of elutrition samples collected every minutes to mixture components provided a good example of a slow transition of the network behavior across time. the first six samples collected between and minutes were assigned to one mixture component and six samples obtained from to minutes were assigned to the other mixture component. samples collected at intermediate time points had a similar odds ratio to pertain to either mixture component. a discussion of the assignment of samples in the α-factor, cdc and cdc synchronization is provided in additional file overall, the mixture bayesian network approach detected more confirmed relationships than the approach presented by friedman et al.. however, friedman et al. detected five relationships that were not predicted by the mixture bayesian approach and were not present in the biogrid database. the higher number of relationships predicted by the mixture bayesian network approach relative to the bayesian network implementations of friedman et al. may be due to the flexibility of the mixture model in the mixture bayesian network approach to accommodate fluctuations on the gene patterns and relationships across the multiple conditions evaluated. single, fixed-number, and optimal mixture densities the mixture bayesian network approach to infer gene networks implemented in this study used the em algorithm to obtain parameter estimates and the bic score to identify the optimal number of mixture components supported by the data. this framework offered great flexibility to accommodate a wide range of models and probability density functions across gene sub-networks and within gene sub-network, across the conditions under study. a series of benchmarking tests were conducted to evaluate the benefits of the more flexible and data-driven mixture bayesian network approach for different scenarios. the description of the sub-networks supported by the data based on bic score favored mixture over single distributions. this indicates that for the networks and data sets considered in our study, mixtures offer a better description of the relationship between genes than single distributions. this result was further validated when the assignment of samples from different treatments or conditions to mixture components was investigated. the robustness of the mixture bayesian network learning process granted by the use of mixture distributions and potentially variable numbers of mixture components across the network was investigated. figure presents the bic scores of the networks predicted using the mixture bayesian network approach with the number of mixture components fixed to single or mixture of two or three gaussian densities, or estimated from the data for the circadian rhythm and adherens junction networks previously presented. the bic score of the optimal mixture implementation was higher than the approaches that used a fixed number of mixtures across all subnetworks. the superiority of the mixture bayesian network approach that estimates the optimal number of mixture components for each gene sub-network from the data over approaches with fixed number of components was evident in both networks. in addition, the maximum number of parents in a gene sub-network was set to five in this study, based on a review of pathways reported in the kegg database and associated to the pathways studied. this maximum was supported by the fact that the optimal bic score was always found with one to four parent genes. computational evaluation the computational time required for the prediction of a network across a wide range of number of genes was evaluated. the mixture bayesian network approach was implemented using the c++ programming language. figure shows the computational time required to reconstruct networks of different size, ranging from five to genes sampled at random from the circadian rhythm and adherens junction microarray data sets. five data sets were obtained from the complete microarray data sets for each network size considered, and a network was estimated for each data set and network size. the computing time to predict a network depends in part on the degree of connectivity between the genes studied. because the gene sets were drawn at random from the microarray data sets, the degree of relationship between these genes was expected to vary across data sets within network size and across network sizes studied. thus, both the number of edges predicted in each network and the computational time required to predict the network were recorded. figure and figure present the average computational time and number of edges across the five data sets within network size evaluated for two cases, when the number of mixture components was estimated from the data or fixed at components. these comparisons allowed the assessment of the impact of the mixture component of the algorithm on the computational time and inferred edges. as expected, the higher the number of genes in the network, the higher the computational time required for estimation of the complete network. the lower computational time to infer the adherens junction network compared to the circadian rhythm network may be due to the fewer relationships detected in the first network. in the adherens junction network, the trends of computational time and number of edges across number of genes were similar. in the circadian rhythm network, the computing time increased at a slightly or moderately higher rate than the number of edges across the number of genes. the disparity between the trends of time and number of edges was higher in the adherens junction network with a fixed number of mixture components than when the number of components was estimated from the data. this result may be due to the potential inadequacy of the fixed number of mixture components to describe some or all gene sub-networks. thus, more iterations and more time were required to estimate the model parameters and infer the network structure. furthermore, the nearly constant number of edges inferred in the adherens junction across the number of gene nodes studied suggested that the mixture bayesian network algorithm was not biased towards inferring higher number of edges with higher number of gene nodes network. algorithmic extensions the identification of a bayesian network structure is an np-hard problem because the number of possible candidate structures super-exponentially increases with the number of gene nodes. two structure learning methods that have complementary advantages were evaluated in this study, the sparse candidate algorithm and the simulated annealing method. the sparse candidate algorithm achieves fast structure learning by restricting the search space of potential parent gene nodes. intuitively we would like to restrict the search space as less as possible. in this study, the algorithm was modified to minimize the impact of restricting the search space, while also maintaining efficiency by evaluating additional parent genes for each gene sub-network with less than the maximum number of parent genes. the network structure was also inferred using a simulated annealing algorithm to overcome the potential risk of inferring networks at local maxima instead of global maximum. at each iteration, the simulated annealing algorithm randomly jumped from the current network topology to a nearby topology formed from the addition, deletion, or reversion of an edge in the current topology. simulated annealing can escape from potential local minima by occasionally sampling from less suitable networks and the probability of jumping to another status is empirically adjusted. this procedure is repeated until convergence. the implementation of these two algorithms allowed the evaluation of the mixture bayesian network approach to be independent of the actual algorithm used to learn the topology of the network. the vast majority of the gene relationships was predicted by both the sparse candidate and simulated annealing algorithms. therefore, only results from the sparse candidate algorithm are presented here. the mixture bayesian approach presented here can be extended to include other gene selection criteria, consider other distributional assumptions and prior information on the relationship between genes in a network. these extensions are discussed in additional file  a comprehensive characterization of the circadian rhythm pathway in honey bees and of the adherens junction pathway in mouse embryos was obtained using a mixture bayesian network approach. all the gene relationships detected by the mixture bayesian network based on honey bee gene expression information were consistent with known relationships in the fruit fly circadian rhythm pathway. the mixture bayesian network approach was also able to uncover changes in gene relationships and profiles in the cyc gene sub-network associated with changes on the maturation age, subspecies of honey bee sampled, and subspecies of honey bees in the host hive. likewise, all the gene relationships detected by the mixture bayesian network using mouse-embryo gene expression information were consistent with known relationships in the general mouse adherens junction pathway. all or the vast majority of the samples within experiment were assigned to the same mixture component model in seven out of the nine experiments considered. this assignment of samples to mixture components suggested that experiment had a major impact on the gene relationships and expression profiles in the actin gene sub-network. the application of the mixture bayesian network approach to the well-studied yeast cell cycle confirmed many reported relationships. in addition, the assignment of samples to the components of the mixture model for the cdc gene sub-network provided novel insights into gene associations. for example, the assignment of elutrition samples to the mixture components suggested a slow transition of the network behavior across time. on the other hand, the assignment of cdc15-synchronized samples to mixture components suggested that the networks were instable across most of the time-points considered. this study demonstrated that the mixture bayesian network approach is well-suited to infer gene pathways based on microarray gene expression data. the estimation of the number of mixture components from the data allowed the simultaneous modeling of gene sub-networks with either constant or variable behavior across conditions. this modeling approach has two potential advantages, a) it can increase the accuracy of the inferred relationships and precision of the parameter estimates, and b) it can help identify groups or clusters of observations that support particular models or mixture components. the superior accuracy of the mixture bayesian network approach was demonstrated by the better bic score of the network with the number of mixture components estimated from the data, compared to the bic score of approaches with a fixed number of mixture components. the mixture of gaussian distributions demonstrated in this study naturally accommodated the continuous and typically normally-distributed gene expression measurements. the proposed approach can be extended to use mixtures of other discrete or continuous distributions, or incorporate prior information on the genes in the network. the proposed mixture bayesian network approach allows the detection of general and condition-specific relationships among genes, and the prediction of an overall network structure within the framework of probabilistic inference. the interpretation of parameter estimates and functions can offer a more comprehensive understanding of gene patterns and relationships. insights into the gene patterns, and the relationships between them, can be used to design effective follow-up experiments and investigate the novel relationships predicted by the mixture bayesian network approach. authors' contributions yk implemented the mixture bayesian network algorithm, inferred the three gene networks, contributed to the interpretation of results, and drafted the manuscript. c-xz helped interpret the algorithm and reviewed the manuscript. srz obtained funding for the study, participated in its conception, coordination, interpretation of results, and helped write the manuscript. all authors have read and approved the final version of this manuscript. supplementary material additional file 1 representation of a mixture model for the bayesian sub-network corresponding to child gene genej and parent genes: gene gene..., genep- the figures represent the bayesian sub-network of the jth gene with p- parent genes. click here for file additional file 2 illustration of the network-structure learning process. diagram of the steps of the algorithm to infer the topology of the network. click here for file additional file 3 estimation of mixture parameters. formulae to estimate the parameters using the expectation-maximization algorithm. click here for file additional file 4 assignment of yeast cell cycle conditions to mixture components. distribution of yeast cell cycle treatments across the mixture components of the gene sub-network. click here for file additional file 5 additional algorithmic extensions. potential modifications of the proposed algorithm are discussed. click here for file acknowledgements we would like to thank heather adams, jaebum kim and bruce southey for useful comments. the support of nih/nigms, nsf/itr and, nih nida are greatly appreciated.
first, the identification of the optimal number of mixture components supported by the data offered a robust approach to infer gene relationships and estimate gene expression profiles. second, the classification of conditions and observations into groups that support particular mixture components helped to uncover both gene relationships that are unique or common across conditions. the mixture models offered a comprehensive description of the relationships among genes in complex biological processes or across a wide range of conditions. the mixture bayesian network approach identified all or the vast majority of the gene relationships reported in empirical studies. findings across the three pathways and data sets indicate that the mixture bayesian network approach is well-suited to infer gene pathways based on microarray data. results from the application of mixture bayesian networks substantially augmented the understanding of gene networks and demonstrated the added-value of this methodology to infer gene networks. inference of gene networks typically relies on measurements across a wide range of conditions or treatments. the first reconstructions of a circadian rhythm pathway in honey bees and an adherens junction pathway in mouse embryos were obtained. in addition, general and condition-specific gene relationships, some unexpected, were detected in these two pathways and in a yeast cell-cycle pathway. this approach integrated bayesian network and gaussian mixture models to describe continuous microarray gene expression measurements, and three gene networks were predicted. the mixture parameter estimates and corresponding odds that the gene network inferred for a sample pertained to each mixture component allowed the uncovering of both general and condition-dependent gene relationships and patterns of expression. a comprehensive approach to infer general and condition-dependent gene networks was evaluated. this study demonstrated the two main benefits of learning gene pathways using mixture bayesian networks. although one network structure is predicted, the relationship between genes could vary across conditions.
 inference of gene networks typically relies on measurements across a wide range of conditions or treatments. although one network structure is predicted, the relationship between genes could vary across conditions. a comprehensive approach to infer general and condition-dependent gene networks was evaluated. this approach integrated bayesian network and gaussian mixture models to describe continuous microarray gene expression measurements, and three gene networks were predicted.  the first reconstructions of a circadian rhythm pathway in honey bees and an adherens junction pathway in mouse embryos were obtained. in addition, general and condition-specific gene relationships, some unexpected, were detected in these two pathways and in a yeast cell-cycle pathway. the mixture bayesian network approach identified all or the vast majority of the gene relationships reported in empirical studies. findings across the three pathways and data sets indicate that the mixture bayesian network approach is well-suited to infer gene pathways based on microarray data. furthermore, the interpretation of model estimates provided a broader understanding of the relationships between genes. the mixture models offered a comprehensive description of the relationships among genes in complex biological processes or across a wide range of conditions. the mixture parameter estimates and corresponding odds that the gene network inferred for a sample pertained to each mixture component allowed the uncovering of both general and condition-dependent gene relationships and patterns of expression.  this study demonstrated the two main benefits of learning gene pathways using mixture bayesian networks. first, the identification of the optimal number of mixture components supported by the data offered a robust approach to infer gene relationships and estimate gene expression profiles. second, the classification of conditions and observations into groups that support particular mixture components helped to uncover both gene relationships that are unique or common across conditions. results from the application of mixture bayesian networks substantially augmented the understanding of gene networks and demonstrated the added-value of this methodology to infer gene networks. 