BACKGROUND
the regulation of gene expression is a critical mechanism in the control of biological processes in cellular organisms.
at the transcriptional level, the main regulators contributing to the control are transcription factors , proteins that bind to cis-regulatory elements in the gene promoter regions  <cit> .
by activating or repressing their target genes, tfs can regulate the global gene expression program of a living cell, and form transcriptional regulatory networks  <cit> .
recent studies have identified that micrornas  play an important role in gene regulation at the post-transcriptional level.
the regulation process takes place via mrna cleavage or translational repression, with mirnas binding to the 3’-untranslated regions  of target mrnas through base pairing to complementary sequences  <cit> .
it has also been demonstrated in a body of literature that mirnas regulate a wide range of biological processes in proliferation  <cit> , metabolism  <cit> , differentiation  <cit> , development  <cit> , apoptosis  <cit> , cellular signaling  <cit>  and even cancer development and progression  <cit> .
it is necessary to draw a unified picture for the regulatory relationships between tfs, mirnas and genes.
however, a challenge is that the combined regulations of mirnas and tfs are complicated, since they involve not only the interactions between each regulator and their target genes, but also the interactions between the regulators themselves.
studies of the gene regulatory networks with the presence of both tfs and mirnas will help elucidate the regulatory mechanisms involving both direct and indirect regulatory relationships.
however, it is still highly unlikely for these relationships to be discovered by biological experiments directly, as the process would be extremely costly and time consuming.
on the other hand, well-designed computational approaches may facilitate the understanding of such complex relationships.
previously, researchers studied the co-regulation of tfs and mirnas by finding out their shared downstream targets  <cit> .
the methods used probabilistic models or statistical tests to measure the significance of the shared targets between the regulators, and to remove the insignificant co-regulating interactions that occurred by chance.
gene enrichment analysis was used in  <cit>  to identify significant co-regulation between the transcriptional and post-transcriptional layers.
they found that some biological processes emerged only in co-regulation and that the disruption of co-regulation may be closely related to cancers, suggesting the importance of the co-regulation of mirnas and tfs.
in  <cit>  available predicted targets databases are used to construct the network, and then gene ontology  was used to discover the significant functional co-regulation pairs.
tran et al.
<cit>  proposed a rule based method to discover the gene regulatory modules that consist of mirnas, tfs, and their target genes based on the available predicted target binding information.
le béchec et al.
<cit>  integrated available target prediction databases to construct a regulatory network that involves mirnas, tfs, and mrnas.
this work provides a good resource for exploring the regulatory relationships or identifying the network motifs.
however, target prediction based on sequences have high rate of false discoveries, which affect the quality of the discoveries of the above mentioned methods.
it would be ideal if expression data can be used to refine the discoveries.
roqueiro et al.
<cit>  proposed a method to identify the key regulators  of pathways.
the method used bayesian inference on known pathway structures to infer a set of regulators in the pathway network.
the bayesian network in this method was constructed manually using the known kegg pathways by removing the cycles in the pathways and applying some filtering criteria.
the method drew findings based on existing knowledge and provided a good resource for other methods to validate their results.
however, it may not be good for exploratory study.
recently, huang et al.
<cit>  developed a web tool  for constructing the regulatory networks that include mirna, tfs, and mrnas.
the built networks can be further analysed to identify network motifs.
the method has used both predicted targets and expression data to build the network.
the method integrated the association network based on expression data and the prior network based on sequence data.
however, an edge in this network shows association, which may not indicate a regulation relationship.
a strong association of a and b may be a result of a common regulator which regulates both a and b. zacher et al.
<cit>  proposed a bayesian inference method based on expression data to explain the activity of mirnas and tfs.
however, this approach does not take into account the interactions between mirnas and tfs.
in this paper, we present a framework to construct the complex regulatory network with three components: tfs, mirnas, and target genes.
our approach aims to discover the regulatory relationships of mirnas and tfs on their target genes respectively, as well as the interplay between the two different types of regulators.
the method utilises multiple sources of data, including gene expression data, target information of each regulator based on sequence data, and sample categories .
to test the proposed method, we use the expression data from the nci- <dig> panel of cell lines  <cit> , and investigate the interactions that may involve in the biological process of epithelial-to-mesenchymal transition .
methods
notation and definitions
consider three expression data sets profiling k mirnas, i tfs, and j mrnas across s samples, respectively.
let x={xk},y={yi},z={zj} be the vectors of mirnas, tfs, and mrnas, respectively, where 1≤k≤k,1≤i≤i, and 1≤j≤j.
each sample is labelled by its category, i.e.
the biological condition of the samples, such as cancer or normal.
in this paper, our goal is to discover the interactions between x,y,z  supported by the expression data and under the constraint of target information .
target information for a regulator is the interactions between the regulator and the regulated genes that are predicted based on the sequence data.
we are particularly interested in the interactions between x and y , and network motifs, which are patterns of subgraphs that recur at frequencies much higher than those found in randomised networks  <cit> .
method overview
in the remaining parts of the methods section, we present our framework for constructing the regulatory network with the co-existence of both regulators, tfs and mirnas.
the method aims to produce regulatory networks including mirnas, tfs, and genes that are relevant to diseases.
the overall process is shown in figure  <dig>
there are three main steps in the framework:  data preparation,  network learning and integration, and  network inferences.
in step , we prepare the input for the network structure learning, including collecting target information for tfs and mirnas, normalising expression data, and analysing differentially expressed genes.
at the beginning of step , the target information is transformed into the  <dig> types of network sub-structures , which are used as the initial structure for the bayesian network learning process .
additionally the expression datasets are split according to sample conditions.
the initial structure are evaluated based on the expression profiles in each condition.
the bayesian networks learnt under each condition are integrated using a bootstrapping and averaging procedure.
therefore the result of step  is an integrated global network with three components: mirnas, tfs, and mrnas.
in the network inference step ), we search the global network for the subgraphs that show the interplay between mirnas and tfs, and network motifs that involve at least two regulators.
in the following, we describe each of the three steps in detail.
step : data preparation
epithelial-to-mesenchymal transition  is part of the process of tissue remodeling during embryonic development and wound healing  <cit> , and during carcinogenesis  <cit>  when cancer cells undergo a change transforming into a more invasive tumor  <cit> .
after emt induction, cells lose their epithelial features characterised by the high e-cadherin expression level, and acquire mesenchymal characteristics, including vimentin filaments and a flattened phenotype.
by expressing proteases, cells become more invasive, and they can pass through the underlying basement membrane and migrate.
these are crucial steps in the multistep process of metastasis  <cit> .
data used in this study contain mirna expression profiles for the nci- <dig> panel of  <dig> cancer cell lines obtained from sø et al.
<cit>  .
the mrna expression profiles for nci- <dig> are downloaded from arrayexpress , accession number e-geod- <dig>  cell lines categorised as epithelial  and mesenchymal  are used in this work.
in order to identify all the tf genes in the data sets, we use the list of tf repertoire mined from  <cit> .
this list is then used to query against the mrna expression profiles from nci- <dig> to extract tf expression profiles.
after normalising the expression data of mirnas, tfs, and mrnas, differentially expressed gene analysis is conducted respectively to all the three components, tfs, mirnas, and mrnas.
the differentially expressed genes between epithelial and mesenchymal samples are identified using the limma package of bioconductor  <cit>  with the benjamini and hochberg’s  correction method  <cit> .
<dig> probes of tfs, and  <dig> probes of mrnas are identified as differentially expressed at significant levels .
also  <dig> probes of mirnas are identified with adjusted p-value < <dig> .
the reason for choosing adjusted  <dig>  as the cut-off for mirna differentially expressed analysis is that the b statistic value output from limma changes the value significantly between adjusted p-value < <dig>  and adjusted p-value > <dig> .
this is a good sign for using the value of  <dig>  as a cut-off, and the number of mirnas obtained with this cut-off is also reasonable for analysing the results.
.
the data input to the bayesian network learning in the next step is the expression profiles of three components, mirnas, tfs, and mrnas.
to integrate the data profiled from different platforms, we discretise the expression values of each gene in each sample to binary values  by using the median of each array as the cut-off.
another input to the bayesian network learning is the target information, which is used as the prior knowledge  to reduce the search space in the learning.
mirna targets and tf targets are collected via commonly used databases.
these databases usually predict target genes using sequence data.
in this paper, we are particularly interested in the information of tfs targeting both mrna and mirna genes, and the mirnas targeting mrna and tf genes.
we use transfac  <dig>   <cit>  and the promoter database  <cit>  to generate tf target information.
tf target information for tf-mrna and tf-tf pairs is obtained from crsd  <cit> , the database utilising and integrating six well-known large scale databases, including transfac  <dig>  and promoter databases.
to obtain the tf-mirna target information, we use mir@nt@n downloaded from  <cit> .
meanwhile mirbase v <dig>   <cit>  from the bioconductor package rmir.hs.mirna  <dig>  is used to build the putative target for mirnas.
the detailed results of all target information are shown in the additional file  <dig>
step : bayesian networks structure learning and integration
bayesian network structure learning has been widely used for discovering gene-gene interaction networks  <cit> , but not often for the discoveries of the interactions between regulators and their target genes.
to represent the interactions between regulators  and between the regulators and their target genes in a bayesian network, regulators and target genes are denoted by nodes and regulatory interactions are denoted by directed edges.
when the expression data of regulators and target genes are given, we can use bayesian network structure learning to discover the interactions.
to start the learning process, we use the target information of regulators to initialise the search space.
therefore in this step, we take the expression profiles and target information as the input to produce a network that indicates the interactions between mirna-tf, mirna-mrna, tf-mirna, tf-tf, and tf-mrna.
the following four sub-steps are carried to obtain the network.
step : sample splitting
to explore all possible interactions including up-, down-, and mix- regulations  in different biological conditions, in  <cit>  we developed the “splitting and averaging” strategy for bayesian networks structure learning.
this strategy splits samples in a data set by their categories of biological conditions.
bayesian network structure learning is used to learn the networks from the samples of each condition respectively, and these networks are then integrated or averaged into a single network.
in our current problem, we firstly use this strategy to learn the networks for the epithelial and mesenchymal conditions separately, then obtain the integrated network from the networks learnt under the two conditions.
so in this sub-step, we split each of the three expression datasets according to sample conditions,  <dig> samples in epithelial and  <dig> samples in mesenchymal .
step : creating the initial structure
to learn a bayesian network with n variables or nodes, the search space, if not restricted, will be all the possible networks formed with the variables.
it has been shown that finding the best network from all the networks is np-hard  <cit> .
therefore in this paper, we assume that the regulatory relationships between regulators and their target genes form a bipartite graph.
this assumption reduces the search space significantly.
moreover, we use target information to initialise the network structure and the topology of the network structure is further constrained.
therefore, we are able to discover the optimal solutions using the exhaustive search method in the given search space.
graphically, the target information can be represented using bipartites as illustrated in figure  <dig>  there are  <dig> types of such bipartites or sub-structures, corresponding to our initial knowledge of the interactions of mirna-tf, mirna-mrna, tf-tf, tf-mirna and tf-mrna.
these bipartites are used as the initial structure for the bayesian network learning.
step : bayesian network learning process
each interaction in the initial structure is evaluated based on the expression data, and the high-confidence interactions are retained.
the learning process searches through all possible candidate structures and evaluates the interactions with a bayesian scoring function.
the candidate structures are generated by removing edges from the initial structure one by one.
the scoring function measures the degree of fitness of any candidate structure g to the dataset.
the goal is to select the structure that best fits the data.
in other words, we need to calculate the probability of each candidate structure g given the data d, p. according to bayes’ theorem, we have:
p=ppp
in the above formula, the term p is the same for all candidate structures.
regarding the term p, it is quite common to assume a uniform distribution  <cit> , and thus it is a constant.
therefore, for comparative purposes, it is sufficient to calculate only p for the scoring function.
in this paper, we use the bde  scoring function developed by heckerman et al.
<cit>  in the following.
scoreb=p=∏i=1n∏j=1qiΓ)Γ+mij)×∏k=1riΓ+sijk)Γ)
where:
n is the number of variables , x <dig> x <dig> …,xn,
qi is the number of different instantiations of the parent of a variable xi in g,
ri is the number of possible values of xi in g,
aijk are the hyperparameters for the dirichlet prior distributions of the parameters given the network structure.
sijk are the corresponding observations from data,
nij=∑kaijk, and mijk=∑ksijk.
more details of the bayesian scoring function can be found in  <cit> .
in practice, we use the bayes net toolbox for matlab   <cit> , and the bde scoring function is included in bnt.
in the next step ) we evaluate the confidence levels of the interactions output from the bayesian network structure learning.
for illustration purpose, consider the learning procedure for the interaction between a regulator r <dig> and its target gene g <dig>  assume that in total r <dig> has two targets and the corresponding initial structure is in figure 2a.
the interactions in each of the four possible candidate structures  can be evaluated with the scoring function based on expression data.
the scores will decide if there is no edge between r <dig> and g <dig> or an edge between them.
in this example, there are two structures with an edge between r <dig> and g <dig>  and two structures with no edge between them.
the average score in each of the two cases is calculated and the structure with higher average score /2=- <dig> ) is chosen .
step : integrating and bootstrapping
it is common to have small number of samples in the dataset of a typical microarray experiment, which unfortunately cannot support statistically significant discoveries.
to overcome this problem, bootstrapping  <cit>  is usually used to improve the confidence of discoveries.
especially, in bayesian network structure learning, bootstrapping can be combined with model  average procedure to discover the interactions with high confidence.
in this paper, the averaging procedure is firstly applied to the bayesian network learning process across different candidate structures.
this procedure is then applied to the sample splitting step across different sample conditions to calculate the average score for each interaction.
next, the score of each interaction is averaged over the number of bootstrapping, and the confidence levels are estimated based on a statistical model as illustrated in the next paragraph.
the interactions with high confidence  are selected to form the integrated network
consider again the example about learning the interaction between r <dig> and g <dig>  let n be the number of bootstrapping iterations, qc be the event of learning the interaction on the local data set dc of the cth condition .
as there are only two possible cases of interactions between r <dig> and g <dig>  we approximate the whole learning process of the interaction between r <dig> and g <dig> as a bernoulli experiment.
we denote qc= <dig> when there is an edge between r <dig> and g <dig> , and assume that p=p= <dig> .
qc follows a binomial distribution qc∼b, as the samples drawn with replacement in the bootstrap are independent  <cit> .
at the integration stage by averaging, the interactions from local data sets dc are aggregated, and the interactions of the regulator r <dig> and its target g <dig> learnt through multiple data sets for the c different conditions  also follows a binomial distribution qr1g1∼b.
adopting this statistical model, we are able to extract the learnt interactions at significant levels.
the interaction that has the confidence level greater than the threshold will be included into the integrated global network.
the matlab codes for the whole process is available on request, and the results for the integrated global network is in additional file  <dig>
step : network inference
a challenging problem of bayesian network learning is the complexity of the resulting networks.
bayesian network learning usually produces a large number of interactions, including false discoveries.
in this step, we extract from the global network learnt in the previous step the interplay between mirnas and tfs.
we search the learnt global network for all of the interactions between mirnas and tfs.
the resulting interplay network will help elucidate the complex regulatory mechanism in specific biological conditions.
in addition to discovering the interplay between mirnas and tfs, we use the network motif finding algorithm, cyclus3d  <cit> , to discover the network motifs that involve at least  <dig> regulators.
network motifs are patterns of sub-graphs that recur at frequencies much higher than those found in randomised networks  <cit> .
the randomised networks satisfy the following criteria: 1) they have the same number of nodes as in the real network, 2) each node in a randomised network has the same number of incoming and outgoing edges as the corresponding node has in the real network, 3) the randomised networks used to calculate the significance of n-node subgraphs were generated to preserve the same number of appearances of all -node subgraphs as in the real network.
these criteria ensure the randomised networks have the similar structure to the real network, and ensure that a high-significance pattern is assigned not because it has a highly significant sub-pattern  <cit> .
the resulting motifs can be considered as simple building blocks from which the network is composed  <cit> , and are believed to have specific functions which play critical roles in biological network inference  <cit> .
the results presented in the next section show that this network inferences step can produce a set of interactions and molecules which are highly relevant to the biological condition of the emt datasets.
RESULTS
the output of the method are two types of networks: 1) the interplays between mirnas and tfs, with their details shown in figure  <dig> and figure 4; 2) the results of network motif finding, which are the feed forward loops  that involve at least two regulators .
from figures  <dig>   <dig>  and  <dig>  we can see that the results generated by our method are compact with only a small number of interactions.
these interactions have been shown to be highly relevant to the biological conditions of emt, and also several emt bio-markers which have been confirmed by literature are identified by our method.
in the rest of this section, firstly we present the interactions and bio-markers that have been confirmed from literature, then we describe the enrichment analysis we have conducted to show the relevance of identified genes to emt.
confirmed interactions and bio-markers for emt
previous studies  <cit>  have demonstrated that the mir- <dig> family targets the e-cadherin transcriptional repressors zinc finger e-box binding homeobox  <dig> and zeb <dig> for emt.
these results have confirmed the interactions found using our method , where we see that the hsa-mir- <dig> family  regulates both zeb <dig> and zeb <dig>  these interactions are the important interactions that are involved in the process of inhibition and induction of emt.
figure  <dig> shows the process in detail where genes identified by our method are marked with red bars.
apart from the mir- <dig> family, several important transcription factors that act as the bio-markers for emt are also confirmed by our method.
the two transcription factors, zeb <dig> and zeb <dig>  which are regulated by the mir- <dig> family, are the markers in all three subtypes of emt  <cit> .
another transcription factor that plays a crucial role in emt is snai <dig> .
in fact, all known emt events during development, cancer, and fibrosis appear to be associated with snai activation  <cit> .
our results suggest that snai <dig> indirectly regulates zeb <dig> and zeb <dig> by regulating the mir- <dig> family transcript , and in turn the mir- <dig> family regulates zeb <dig> and zeb <dig> .
this result is consistent with the literature as snai <dig> is confirmed to regulate the mir- <dig> family  <cit> .
the other emt bio-marker identified by our method is ets <dig> .
it has been suggested that ets <dig> is an upstream regulator of zeb <dig> and zeb <dig>  <cit> , and therefore plays a critical role in activating the regulation of emt.
functions of identified genes being highly enriched for emt
the functions of the identified genes  and the pathways which the genes potentially constitute are analysed using genego metacore from genego inc. and the ingenuity pathway analysis .
the genes identified as a result of the network inference step are significantly enriched for several biological functions.
the top functions output from ipa that are known to be critical for emt are gene expression, cellular development, cellular growth and proliferation, cellular movement, and cell death.
moreover, several genes belong to the classes of invasion and migration.
these classes are sub-categories of cellular movement, and they have been confirmed as the functional markers of emt  <cit> .
this suggests that many target genes and their interactive regulators are involved in emt.
table  <dig> shows the genes in the class of invasion and migration together with their p-values.
a significant number of genes identified in the inference step belong to the class of invasion and migration which are emt functional markers.
the results are generated by ipa.
in addition, the pathways which the genes in our results potentially constitute are identified using genego metacore.
the statistically mapped pathways show that they are highly relevant to emt.
there are direct pathways that regulate emt, such as the pathway of development_microrna-dependent inhibition of emt.
this pathway shows the regulation of the mir- <dig> family and other mirnas on the emt bio-markers zeb <dig> and zeb <dig>  and results in the inhibition and induction of emt.
figure  <dig> shows the details of this pathway.
other direct pathways such as the development_tgf-beta-dependent induction of emt via smads, and cell adhesion_tight junctions, are known to play critical roles in the regulatory procedure of emt.
the summary of these pathways and the corresponding p-values are given in table  <dig>
the mapped pathways involve identified genes that are important for emt.
the results are generated by genego metacore.
discussion
during the past few decades, considerable efforts have been made to explore the transcriptional regulatory networks in which transcription factors play the role as a main regulator.
other recent studies have investigated the post-transcriptional regulatory networks with mirnas as the main regulator.
however, with the ultimate goal of achieving a profound understanding of the mechanisms that control gene activities, it is sensible and desirable to find regulatory relationships involving both types of regulators from diverse sources of data.
in this paper, we utilise bayesian network learning in constructing the network, but the integrated global network in general is not a bayesian network.
for instance, one of the requirements for bayesian networks is that the network structure must be a directed acyclic graph .
our integrated global network may include some loops of interactions where two regulators interact with each other, hence it is not a bayesian network.
such cyclic network behaviour is more reasonable in reality, as more and more feedback loops between mirnas and tfs are being reported.
for instance, the zeb/mir- <dig> pair is a feedback loop that regulates emt  <cit> .
therefore, the integrated global network may provide more information than normal bayesian networks which are dags.
in the network inference step, we use network motif finding algorithm to discover the sub-networks that recur at statistically significant level.
interestingly, the results from these small sub-networks still retain several important interactions and molecules relevant to the biological condition of the dataset.
the relationships between the significance in frequency of graphs and biological functions are still open and interesting research topics.
in the dataset used in this paper, the results are supportive for this hypothesis.
an advantage of motif finding is that it produces a manageable number of interactions that can be used for further experimentation.
the results from this paper, therefore, can provide good resources for future validating experiments.
while the network motifs found based on the regulatory network may provide useful patterns to guide biological experiments, these motifs depend on the structure of the regulatory network.
the structure of the regulatory network in this paper is obtained based on the assumption that mirnas and tfs are regulators and mrnas are targets.
however, the knowledge of gene regulatory relationships is still limited and the assumption may not always hold in reality.
when the structure of the regulatory network changes the resulting network motifs may change too.
in the paper, we use the differentially expressed genes as the nodes for the gene regulatory network.
we assume that genes whose expression levels do not change significantly between conditions would not play an important role in the regulatory network.
there may be the case that a gene is the target of two regulators that cancel out each other, resulting in the expression level of the target gene unchanged.
however, to make our method computationally practical we do not consider such cases.
to start the process of bayesian network structure learning, target information is used to initialise the network.
the target information based on sequence data may involve false discoveries.
bayesian network structure learning uses gene expression data to evaluate the confidence level of each interaction  in the initial network, and only the interactions of high confidence are integrated into the global network.
therefore, graphically the set of edges in the global network is a subset of the set of edges in the initial network.
the enrichment analysis shows that the important interactions for emt are retained in the global network, demonstrating the effectiveness of the method.
other high-confidence interactions provide strong hypotheses for experimental validations.
CONCLUSIONS
in this study, we have proposed a framework for inferring complex gene regulatory networks using diverse sources of data, including target information for regulators, expression profiles, and sample categories.
the interplay between regulators and the motifs with which they regulate target genes are revealed in the three-component network, and it is impossible to infer the interplay from any single regulator regulatory networks.
the analysis of the emt datasets has produced several results that have been validated by literature, a number of statistically significant interactions between mirnas and tfs, and novel network motifs.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
tl and jyl conceived this research.
tl designed and performed the experiments.
bl, at and gg provided the source of data and validated the results.
ll verified the learning model.
tl, ll, bl, ks and jyl drafted the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
differentially expressed mirnas, tfs, and mrnas.
limma package from bioconductor is used to identify differentially expressed mirnas, tfs, and mrnas.
click here for file
additional file 2
target information.
the interactions that show the tf and mirna target information.
this information is used to initialise the searching space for bayesian network learning.
click here for file
additional file 3
significant interactions.
all the statistically significant interactions for the complex three-component network.
these interactions represent the regulatory relationships between mir-mrna, mir-tf, tf-mirna, tf-tf, and tf-mrna.
click here for file
acknowledgements
this work has been partially supported by australian research council discovery project dp <dig> and national health and medical research council project grant app <dig>

 the regulation of gene expression is a critical mechanism in the control of biological processes in cellular organisms. at the transcriptional level, the main regulators contributing to the control are transcription factors, proteins that bind to cis-regulatory elements in the gene promoter regions. by activating or repressing their target genes, tfs can regulate the global gene expression program of a living cell, and form transcriptional regulatory networks. recent studies have identified that micrornas play an important role in gene regulation at the post-transcriptional level. the regulation process takes place via mrna cleavage or translational repression, with mirnas binding to the 3’-untranslated regions of target mrnas through base pairing to complementary sequences. it has also been demonstrated in a body of literature that mirnas regulate a wide range of biological processes in proliferation, metabolism, differentiation, development, apoptosis, cellular signaling and even cancer development and progression. it is necessary to draw a unified picture for the regulatory relationships between tfs, mirnas and genes. however, a challenge is that the combined regulations of mirnas and tfs are complicated, since they involve not only the interactions between each regulator and their target genes, but also the interactions between the regulators themselves. studies of the gene regulatory networks with the presence of both tfs and mirnas will help elucidate the regulatory mechanisms involving both direct and indirect regulatory relationships. however, it is still highly unlikely for these relationships to be discovered by biological experiments directly, as the process would be extremely costly and time consuming. on the other hand, well-designed computational approaches may facilitate the understanding of such complex relationships. previously, researchers studied the co-regulation of tfs and mirnas by finding out their shared downstream targets. the methods used probabilistic models or statistical tests to measure the significance of the shared targets between the regulators, and to remove the insignificant co-regulating interactions that occurred by chance. gene enrichment analysis was used in to identify significant co-regulation between the transcriptional and post-transcriptional layers. they found that some biological processes emerged only in co-regulation and that the disruption of co-regulation may be closely related to cancers, suggesting the importance of the co-regulation of mirnas and tfs. in available predicted targets databases are used to construct the network, and then gene ontology was used to discover the significant functional co-regulation pairs. tran et al. proposed a rule based method to discover the gene regulatory modules that consist of mirnas, tfs, and their target genes based on the available predicted target binding information. le béchec et al. integrated available target prediction databases to construct a regulatory network that involves mirnas, tfs, and mrnas. this work provides a good resource for exploring the regulatory relationships or identifying the network motifs. however, target prediction based on sequences have high rate of false discoveries, which affect the quality of the discoveries of the above mentioned methods. it would be ideal if expression data can be used to refine the discoveries. roqueiro et al. proposed a method to identify the key regulators of pathways. the method used bayesian inference on known pathway structures to infer a set of regulators in the pathway network. the bayesian network in this method was constructed manually using the known kegg pathways by removing the cycles in the pathways and applying some filtering criteria. the method drew findings based on existing knowledge and provided a good resource for other methods to validate their results. however, it may not be good for exploratory study. recently, huang et al. developed a web tool for constructing the regulatory networks that include mirna, tfs, and mrnas. the built networks can be further analysed to identify network motifs. the method has used both predicted targets and expression data to build the network. the method integrated the association network based on expression data and the prior network based on sequence data. however, an edge in this network shows association, which may not indicate a regulation relationship. a strong association of a and b may be a result of a common regulator which regulates both a and b. zacher et al. proposed a bayesian inference method based on expression data to explain the activity of mirnas and tfs. however, this approach does not take into account the interactions between mirnas and tfs. in this paper, we present a framework to construct the complex regulatory network with three components: tfs, mirnas, and target genes. our approach aims to discover the regulatory relationships of mirnas and tfs on their target genes respectively, as well as the interplay between the two different types of regulators. the method utilises multiple sources of data, including gene expression data, target information of each regulator based on sequence data, and sample categories. to test the proposed method, we use the expression data from the nci- panel of cell lines, and investigate the interactions that may involve in the biological process of epithelial-to-mesenchymal transition. methods notation and definitions consider three expression data sets profiling k mirnas, i tfs, and j mrnas across s samples, respectively. let x={xk},y={yi},z={zj} be the vectors of mirnas, tfs, and mrnas, respectively, where 1≤k≤k,1≤i≤i, and 1≤j≤j. each sample is labelled by its category, i.e. the biological condition of the samples, such as cancer or normal. in this paper, our goal is to discover the interactions between x,y,z supported by the expression data and under the constraint of target information. target information for a regulator is the interactions between the regulator and the regulated genes that are predicted based on the sequence data. we are particularly interested in the interactions between x and y, and network motifs, which are patterns of subgraphs that recur at frequencies much higher than those found in randomised networks. method overview in the remaining parts of the methods section, we present our framework for constructing the regulatory network with the co-existence of both regulators, tfs and mirnas. the method aims to produce regulatory networks including mirnas, tfs, and genes that are relevant to diseases. the overall process is shown in figure there are three main steps in the framework: data preparation, network learning and integration, and network inferences. in step, we prepare the input for the network structure learning, including collecting target information for tfs and mirnas, normalising expression data, and analysing differentially expressed genes. at the beginning of step, the target information is transformed into the types of network sub-structures, which are used as the initial structure for the bayesian network learning process. additionally the expression datasets are split according to sample conditions. the initial structure are evaluated based on the expression profiles in each condition. the bayesian networks learnt under each condition are integrated using a bootstrapping and averaging procedure. therefore the result of step is an integrated global network with three components: mirnas, tfs, and mrnas. in the network inference step ), we search the global network for the subgraphs that show the interplay between mirnas and tfs, and network motifs that involve at least two regulators. in the following, we describe each of the three steps in detail. step: data preparation epithelial-to-mesenchymal transition is part of the process of tissue remodeling during embryonic development and wound healing, and during carcinogenesis when cancer cells undergo a change transforming into a more invasive tumor. after emt induction, cells lose their epithelial features characterised by the high e-cadherin expression level, and acquire mesenchymal characteristics, including vimentin filaments and a flattened phenotype. by expressing proteases, cells become more invasive, and they can pass through the underlying basement membrane and migrate. these are crucial steps in the multistep process of metastasis. data used in this study contain mirna expression profiles for the nci- panel of cancer cell lines obtained from sø et al.. the mrna expression profiles for nci- are downloaded from arrayexpress, accession number e-geod- cell lines categorised as epithelial and mesenchymal are used in this work. in order to identify all the tf genes in the data sets, we use the list of tf repertoire mined from. this list is then used to query against the mrna expression profiles from nci- to extract tf expression profiles. after normalising the expression data of mirnas, tfs, and mrnas, differentially expressed gene analysis is conducted respectively to all the three components, tfs, mirnas, and mrnas. the differentially expressed genes between epithelial and mesenchymal samples are identified using the limma package of bioconductor with the benjamini and hochberg’s correction method. probes of tfs, and probes of mrnas are identified as differentially expressed at significant levels. also probes of mirnas are identified with adjusted p-value <. the reason for choosing adjusted as the cut-off for mirna differentially expressed analysis is that the b statistic value output from limma changes the value significantly between adjusted p-value < and adjusted p-value >. this is a good sign for using the value of as a cut-off, and the number of mirnas obtained with this cut-off is also reasonable for analysing the results.. the data input to the bayesian network learning in the next step is the expression profiles of three components, mirnas, tfs, and mrnas. to integrate the data profiled from different platforms, we discretise the expression values of each gene in each sample to binary values by using the median of each array as the cut-off. another input to the bayesian network learning is the target information, which is used as the prior knowledge to reduce the search space in the learning. mirna targets and tf targets are collected via commonly used databases. these databases usually predict target genes using sequence data. in this paper, we are particularly interested in the information of tfs targeting both mrna and mirna genes, and the mirnas targeting mrna and tf genes. we use transfac and the promoter database to generate tf target information. tf target information for tf-mrna and tf-tf pairs is obtained from crsd, the database utilising and integrating six well-known large scale databases, including transfac and promoter databases. to obtain the tf-mirna target information, we use mir@nt@n downloaded from. meanwhile mirbase v from the bioconductor package rmir.hs.mirna is used to build the putative target for mirnas. the detailed results of all target information are shown in the additional file step: bayesian networks structure learning and integration bayesian network structure learning has been widely used for discovering gene-gene interaction networks, but not often for the discoveries of the interactions between regulators and their target genes. to represent the interactions between regulators and between the regulators and their target genes in a bayesian network, regulators and target genes are denoted by nodes and regulatory interactions are denoted by directed edges. when the expression data of regulators and target genes are given, we can use bayesian network structure learning to discover the interactions. to start the learning process, we use the target information of regulators to initialise the search space. therefore in this step, we take the expression profiles and target information as the input to produce a network that indicates the interactions between mirna-tf, mirna-mrna, tf-mirna, tf-tf, and tf-mrna. the following four sub-steps are carried to obtain the network. step: sample splitting to explore all possible interactions including up-, down-, and mix- regulations in different biological conditions, in we developed the “splitting and averaging” strategy for bayesian networks structure learning. this strategy splits samples in a data set by their categories of biological conditions. bayesian network structure learning is used to learn the networks from the samples of each condition respectively, and these networks are then integrated or averaged into a single network. in our current problem, we firstly use this strategy to learn the networks for the epithelial and mesenchymal conditions separately, then obtain the integrated network from the networks learnt under the two conditions. so in this sub-step, we split each of the three expression datasets according to sample conditions, samples in epithelial and samples in mesenchymal. step: creating the initial structure to learn a bayesian network with n variables or nodes, the search space, if not restricted, will be all the possible networks formed with the variables. it has been shown that finding the best network from all the networks is np-hard. therefore in this paper, we assume that the regulatory relationships between regulators and their target genes form a bipartite graph. this assumption reduces the search space significantly. moreover, we use target information to initialise the network structure and the topology of the network structure is further constrained. therefore, we are able to discover the optimal solutions using the exhaustive search method in the given search space. graphically, the target information can be represented using bipartites as illustrated in figure there are types of such bipartites or sub-structures, corresponding to our initial knowledge of the interactions of mirna-tf, mirna-mrna, tf-tf, tf-mirna and tf-mrna. these bipartites are used as the initial structure for the bayesian network learning. step: bayesian network learning process each interaction in the initial structure is evaluated based on the expression data, and the high-confidence interactions are retained. the learning process searches through all possible candidate structures and evaluates the interactions with a bayesian scoring function. the candidate structures are generated by removing edges from the initial structure one by one. the scoring function measures the degree of fitness of any candidate structure g to the dataset. the goal is to select the structure that best fits the data. in other words, we need to calculate the probability of each candidate structure g given the data d, p. according to bayes’ theorem, we have: p=ppp in the above formula, the term p is the same for all candidate structures. regarding the term p, it is quite common to assume a uniform distribution, and thus it is a constant. therefore, for comparative purposes, it is sufficient to calculate only p for the scoring function. in this paper, we use the bde scoring function developed by heckerman et al. in the following. scoreb=p=∏i=1n∏j=1qiΓ)Γ+mij)×∏k=1riΓ+sijk)Γ) where: n is the number of variables, x x …,xn, qi is the number of different instantiations of the parent of a variable xi in g, ri is the number of possible values of xi in g, aijk are the hyperparameters for the dirichlet prior distributions of the parameters given the network structure. sijk are the corresponding observations from data, nij=∑kaijk, and mijk=∑ksijk. more details of the bayesian scoring function can be found in. in practice, we use the bayes net toolbox for matlab, and the bde scoring function is included in bnt. in the next step ) we evaluate the confidence levels of the interactions output from the bayesian network structure learning. for illustration purpose, consider the learning procedure for the interaction between a regulator r and its target gene g assume that in total r has two targets and the corresponding initial structure is in figure 2a. the interactions in each of the four possible candidate structures can be evaluated with the scoring function based on expression data. the scores will decide if there is no edge between r and g or an edge between them. in this example, there are two structures with an edge between r and g and two structures with no edge between them. the average score in each of the two cases is calculated and the structure with higher average score /2=- ) is chosen. step: integrating and bootstrapping it is common to have small number of samples in the dataset of a typical microarray experiment, which unfortunately cannot support statistically significant discoveries. to overcome this problem, bootstrapping is usually used to improve the confidence of discoveries. especially, in bayesian network structure learning, bootstrapping can be combined with model average procedure to discover the interactions with high confidence. in this paper, the averaging procedure is firstly applied to the bayesian network learning process across different candidate structures. this procedure is then applied to the sample splitting step across different sample conditions to calculate the average score for each interaction. next, the score of each interaction is averaged over the number of bootstrapping, and the confidence levels are estimated based on a statistical model as illustrated in the next paragraph. the interactions with high confidence are selected to form the integrated network consider again the example about learning the interaction between r and g let n be the number of bootstrapping iterations, qc be the event of learning the interaction on the local data set dc of the cth condition. as there are only two possible cases of interactions between r and g we approximate the whole learning process of the interaction between r and g as a bernoulli experiment. we denote qc= when there is an edge between r and g, and assume that p=p=. qc follows a binomial distribution qc∼b, as the samples drawn with replacement in the bootstrap are independent. at the integration stage by averaging, the interactions from local data sets dc are aggregated, and the interactions of the regulator r and its target g learnt through multiple data sets for the c different conditions also follows a binomial distribution qr1g1∼b. adopting this statistical model, we are able to extract the learnt interactions at significant levels. the interaction that has the confidence level greater than the threshold will be included into the integrated global network. the matlab codes for the whole process is available on request, and the results for the integrated global network is in additional file step: network inference a challenging problem of bayesian network learning is the complexity of the resulting networks. bayesian network learning usually produces a large number of interactions, including false discoveries. in this step, we extract from the global network learnt in the previous step the interplay between mirnas and tfs. we search the learnt global network for all of the interactions between mirnas and tfs. the resulting interplay network will help elucidate the complex regulatory mechanism in specific biological conditions. in addition to discovering the interplay between mirnas and tfs, we use the network motif finding algorithm, cyclus3d, to discover the network motifs that involve at least regulators. network motifs are patterns of sub-graphs that recur at frequencies much higher than those found in randomised networks. the randomised networks satisfy the following criteria: 1) they have the same number of nodes as in the real network, 2) each node in a randomised network has the same number of incoming and outgoing edges as the corresponding node has in the real network, 3) the randomised networks used to calculate the significance of n-node subgraphs were generated to preserve the same number of appearances of all -node subgraphs as in the real network. these criteria ensure the randomised networks have the similar structure to the real network, and ensure that a high-significance pattern is assigned not because it has a highly significant sub-pattern. the resulting motifs can be considered as simple building blocks from which the network is composed, and are believed to have specific functions which play critical roles in biological network inference. the results presented in the next section show that this network inferences step can produce a set of interactions and molecules which are highly relevant to the biological condition of the emt datasets.  the output of the method are two types of networks: 1) the interplays between mirnas and tfs, with their details shown in figure and figure 4; 2) the results of network motif finding, which are the feed forward loops that involve at least two regulators. from figures and we can see that the results generated by our method are compact with only a small number of interactions. these interactions have been shown to be highly relevant to the biological conditions of emt, and also several emt bio-markers which have been confirmed by literature are identified by our method. in the rest of this section, firstly we present the interactions and bio-markers that have been confirmed from literature, then we describe the enrichment analysis we have conducted to show the relevance of identified genes to emt. confirmed interactions and bio-markers for emt previous studies have demonstrated that the mir- family targets the e-cadherin transcriptional repressors zinc finger e-box binding homeobox and zeb for emt. these results have confirmed the interactions found using our method, where we see that the hsa-mir- family regulates both zeb and zeb these interactions are the important interactions that are involved in the process of inhibition and induction of emt. figure shows the process in detail where genes identified by our method are marked with red bars. apart from the mir- family, several important transcription factors that act as the bio-markers for emt are also confirmed by our method. the two transcription factors, zeb and zeb which are regulated by the mir- family, are the markers in all three subtypes of emt. another transcription factor that plays a crucial role in emt is snai. in fact, all known emt events during development, cancer, and fibrosis appear to be associated with snai activation. our results suggest that snai indirectly regulates zeb and zeb by regulating the mir- family transcript, and in turn the mir- family regulates zeb and zeb. this result is consistent with the literature as snai is confirmed to regulate the mir- family. the other emt bio-marker identified by our method is ets. it has been suggested that ets is an upstream regulator of zeb and zeb, and therefore plays a critical role in activating the regulation of emt. functions of identified genes being highly enriched for emt the functions of the identified genes and the pathways which the genes potentially constitute are analysed using genego metacore from genego inc. and the ingenuity pathway analysis. the genes identified as a result of the network inference step are significantly enriched for several biological functions. the top functions output from ipa that are known to be critical for emt are gene expression, cellular development, cellular growth and proliferation, cellular movement, and cell death. moreover, several genes belong to the classes of invasion and migration. these classes are sub-categories of cellular movement, and they have been confirmed as the functional markers of emt. this suggests that many target genes and their interactive regulators are involved in emt. table shows the genes in the class of invasion and migration together with their p-values. a significant number of genes identified in the inference step belong to the class of invasion and migration which are emt functional markers. the results are generated by ipa. in addition, the pathways which the genes in our results potentially constitute are identified using genego metacore. the statistically mapped pathways show that they are highly relevant to emt. there are direct pathways that regulate emt, such as the pathway of development_microrna-dependent inhibition of emt. this pathway shows the regulation of the mir- family and other mirnas on the emt bio-markers zeb and zeb and results in the inhibition and induction of emt. figure shows the details of this pathway. other direct pathways such as the development_tgf-beta-dependent induction of emt via smads, and cell adhesion_tight junctions, are known to play critical roles in the regulatory procedure of emt. the summary of these pathways and the corresponding p-values are given in table the mapped pathways involve identified genes that are important for emt. the results are generated by genego metacore. discussion during the past few decades, considerable efforts have been made to explore the transcriptional regulatory networks in which transcription factors play the role as a main regulator. other recent studies have investigated the post-transcriptional regulatory networks with mirnas as the main regulator. however, with the ultimate goal of achieving a profound understanding of the mechanisms that control gene activities, it is sensible and desirable to find regulatory relationships involving both types of regulators from diverse sources of data. in this paper, we utilise bayesian network learning in constructing the network, but the integrated global network in general is not a bayesian network. for instance, one of the requirements for bayesian networks is that the network structure must be a directed acyclic graph. our integrated global network may include some loops of interactions where two regulators interact with each other, hence it is not a bayesian network. such cyclic network behaviour is more reasonable in reality, as more and more feedback loops between mirnas and tfs are being reported. for instance, the zeb/mir- pair is a feedback loop that regulates emt. therefore, the integrated global network may provide more information than normal bayesian networks which are dags. in the network inference step, we use network motif finding algorithm to discover the sub-networks that recur at statistically significant level. interestingly, the results from these small sub-networks still retain several important interactions and molecules relevant to the biological condition of the dataset. the relationships between the significance in frequency of graphs and biological functions are still open and interesting research topics. in the dataset used in this paper, the results are supportive for this hypothesis. an advantage of motif finding is that it produces a manageable number of interactions that can be used for further experimentation. the results from this paper, therefore, can provide good resources for future validating experiments. while the network motifs found based on the regulatory network may provide useful patterns to guide biological experiments, these motifs depend on the structure of the regulatory network. the structure of the regulatory network in this paper is obtained based on the assumption that mirnas and tfs are regulators and mrnas are targets. however, the knowledge of gene regulatory relationships is still limited and the assumption may not always hold in reality. when the structure of the regulatory network changes the resulting network motifs may change too. in the paper, we use the differentially expressed genes as the nodes for the gene regulatory network. we assume that genes whose expression levels do not change significantly between conditions would not play an important role in the regulatory network. there may be the case that a gene is the target of two regulators that cancel out each other, resulting in the expression level of the target gene unchanged. however, to make our method computationally practical we do not consider such cases. to start the process of bayesian network structure learning, target information is used to initialise the network. the target information based on sequence data may involve false discoveries. bayesian network structure learning uses gene expression data to evaluate the confidence level of each interaction in the initial network, and only the interactions of high confidence are integrated into the global network. therefore, graphically the set of edges in the global network is a subset of the set of edges in the initial network. the enrichment analysis shows that the important interactions for emt are retained in the global network, demonstrating the effectiveness of the method. other high-confidence interactions provide strong hypotheses for experimental validations.  in this study, we have proposed a framework for inferring complex gene regulatory networks using diverse sources of data, including target information for regulators, expression profiles, and sample categories. the interplay between regulators and the motifs with which they regulate target genes are revealed in the three-component network, and it is impossible to infer the interplay from any single regulator regulatory networks. the analysis of the emt datasets has produced several results that have been validated by literature, a number of statistically significant interactions between mirnas and tfs, and novel network motifs. competing interests the authors declare that they have no competing interests. authors’ contributions tl and jyl conceived this research. tl designed and performed the experiments. bl, at and gg provided the source of data and validated the results. ll verified the learning model. tl, ll, bl, ks and jyl drafted the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 differentially expressed mirnas, tfs, and mrnas. limma package from bioconductor is used to identify differentially expressed mirnas, tfs, and mrnas. click here for file additional file 2 target information. the interactions that show the tf and mirna target information. this information is used to initialise the searching space for bayesian network learning. click here for file additional file 3 significant interactions. all the statistically significant interactions for the complex three-component network. these interactions represent the regulatory relationships between mir-mrna, mir-tf, tf-mirna, tf-tf, and tf-mrna. click here for file acknowledgements this work has been partially supported by australian research council discovery project dp and national health and medical research council project grant app 
this method firstly utilises bayesian network structure learning to construct a regulatory network from multiple sources of data: gene expression profiles of mirnas, tfs and mrnas, target information based on sequence data, and sample categories. we have designed a framework to infer gene regulatory networks involving both tfs and mirnas from multiple sources of data, including gene expression data, target information, and sample categories. we apply the proposed framework to the data sets of epithelial-to-mesenchymal transition. then, in order to produce more meaningful results for further biological experimentation and research, the method searches the learnt network to identify the interplay between mirnas and tfs and applies a network motif finding algorithm to further infer the network. regulatory mechanisms of the two main regulators are of great interest to biologists and may provide insights into the causes of diseases. currently, it is very difficult to study the regulatory mechanisms that involve both mirnas and tfs in a biological lab. however, networks involving a single type of regulator may not fully reveal the complex gene regulatory mechanisms, for instance, the way in which a tf indirectly regulates a gene via a mirna. the results elucidate the complex gene regulatory mechanism for emt which involves both tfs and mirnas. moreover, the results generated by our method are compact, involving a small number of interactions which have been proved highly relevant to emt. results on the emt data sets have shown that the proposed approach is able to produce compact and meaningful gene regulatory networks that are highly relevant to the biological conditions of the data sets. however, the interplay between mirnas and tfs in a regulatory network still remains unearthed. previous research has been mostly directed at inferring either mirna or tf regulatory networks from data. we propose a framework to learn from heterogeneous data the three-component regulatory networks, with the presence of mirnas, tfs, and mrnas. transcription factors and micrornas are primary metazoan gene regulators.
 transcription factors and micrornas are primary metazoan gene regulators. regulatory mechanisms of the two main regulators are of great interest to biologists and may provide insights into the causes of diseases. however, the interplay between mirnas and tfs in a regulatory network still remains unearthed. currently, it is very difficult to study the regulatory mechanisms that involve both mirnas and tfs in a biological lab. even at data level, a network involving mirnas, tfs and genes will be too complicated to achieve. previous research has been mostly directed at inferring either mirna or tf regulatory networks from data. however, networks involving a single type of regulator may not fully reveal the complex gene regulatory mechanisms, for instance, the way in which a tf indirectly regulates a gene via a mirna.  we propose a framework to learn from heterogeneous data the three-component regulatory networks, with the presence of mirnas, tfs, and mrnas. this method firstly utilises bayesian network structure learning to construct a regulatory network from multiple sources of data: gene expression profiles of mirnas, tfs and mrnas, target information based on sequence data, and sample categories. then, in order to produce more meaningful results for further biological experimentation and research, the method searches the learnt network to identify the interplay between mirnas and tfs and applies a network motif finding algorithm to further infer the network. we apply the proposed framework to the data sets of epithelial-to-mesenchymal transition. the results elucidate the complex gene regulatory mechanism for emt which involves both tfs and mirnas. several discovered interactions and molecular functions have been confirmed by literature. in addition, many other discovered interactions and bio-markers are of high statistical significance and thus can be good candidates for validation by experiments. moreover, the results generated by our method are compact, involving a small number of interactions which have been proved highly relevant to emt.  we have designed a framework to infer gene regulatory networks involving both tfs and mirnas from multiple sources of data, including gene expression data, target information, and sample categories. results on the emt data sets have shown that the proposed approach is able to produce compact and meaningful gene regulatory networks that are highly relevant to the biological conditions of the data sets. this framework has the potential for application to other heterogeneous datasets to reveal the complex gene regulatory relationships. 