Gene networks as a tool to understand transcriptional regulation

INTRODUCTION

Gene regulatory networks, or simply gene networks (GNs), have shown to be a promising approach that the bioinformatics community has been developing for studying regulatory mechanisms in biological systems. These GNs are built from the genome-wide high-throughput gene expression data (Figure 1) that are often available from DNA microarray experiments (Friedman, 2004). Transcriptome data have become increasingly important because they provide an appropriate source for studies of the systemic behavior of complex biological systems, and the GNs are one of the computational tools that have been employed to this end.

A DNA microarray experiment, which measures the levels of genes being expressed, can be both designed to generate time series expression data (sampling at different times under the same conditions) or it can be designed in a static manner, in which one takes snapshots from cells under a variety of physiological conditions (i.e., from different tissues). The microarray generates the mRNA levels of transcripts, and these data can be used to reveal the hidden “program” that controls gene expression in that biological system.

The expression program of a cell determines which group of genes is going to shift from the basal to active transcription, and which genes are going to be inhibited. Several physical and mechanical systems are modeled as stochastic processes susceptible to noisy measures (Roweis and Ghahramani, 1999). We can also borrow this approach to analyze biological systems, as we consider that the system that we want to model is a program that controls gene expression at the molecular level. In this system, we have a number of regulatory mechanisms, represented by the transcriptional control strategies of the cell, such as the synthesis of specialized proteins as transcription factors. We also know that this expression program works differently according to external stimuli and to intracellular conditions, such as the concentration of some metabolites. Furthermore, the output of this system is the mRNA levels of the genes that are involved. These measures are, in turn, highly noisy, due to experimental and computational factors. In this scenario, the question that arises is “What is the structure of the system that has generated this expression data?” or equivalently, “How do genes interact among themselves to generate this output?” The GNs aim to recover these interactions and therefore, the underlying system that works at the transcriptional level, through so-called probabilistic graphical models.

Conceptually, GNs are (un)directed graphs, in which the nodes correspond to the genes and a link between a pair of genes denotes a regulatory interaction between them (Figure 2). Since we are assessing mRNA profiles, a link between two genes represents the relationship that occurs at the transcriptional level, i.e., a transcription factor that is activating its target gene.

For the purpose of inference of GN from data, graphical models provide the mathematical theory to determine the gene interactions. Since the pioneer work of Friedman et al. (2000), several models have been proposed (for a review, see Xia et al. (2004), van Someren et al. (2002), and de Jong (2002)). Bayesian networks (BN) are of special interest for two main reasons: 1) the ability to handle noisy data, such as microarray outputs, and 2) prior knowledge can easily be incorporated into the model.

We developed a framework for GN reconstruction, based on a BN model that captures direct interactions between genes through nonparametric regression and that has multiple search schemes for learning the network. The base functions used for regression are B-splines. This approach of function approximation by nonparametric regression leads to flexibility of the model, since it does not assume any fixed regulation function. The model has a number of search schemes that are based on two search algorithms, known as random hill-climbing and greedy hill-climbing. This framework has many options for initializing the start graph, including empty, random, score-driven, and partial correlation initialization.

METHODS

Before the mathematical enunciation of the model, it is important to discuss some principles behind GN inference. The main objective of an algorithm to infer GNs is to find the correct set of regulators for the genes in the net. Within the BN terminology, this means finding the correct parents for each node. As an example, we assume that expression levels of a regulated gene a are dictated by its parents genes. If gene a is mainly regulated by genes b and c, on a universe of n genes, the goal is to identify, among all possible combinations of n-1 other genes, that the ideal regulators (parents) for a are b and c. We also should find, from expression data, that the mRNA levels of a have their best approximation when we take them as a function of the parents’ expression, that is, E(a) » F(E(b), E(c)). Hence, we need to make use of a mathematical tool to perform this approximation of gene expression given its parents.

The model uses nonparametric regression with B-splines to capture the relationships between genes, as introduced by Imoto et al. (2002) and further developed by Bastos and Guimarães (2005). While some BN models for GNs use a fixed regulation law (Nachman et al., 2004), our approach is flexible on this issue, not restricting the interactions to behavior to a fixed function.

Bayesian network model for gene networks

Let X = (X₁, X₂, …, X_n)^T be a random n-dimensional vector containing the genes to be analyzed, and assume that G is a directed graph. Under BN theory, genes are random variables, and it is possible to decompose their joint probability into a product of conditional probabilities:

where P_i = (P₁⁽ⁱ⁾ , P₂⁽ⁱ⁾, . . . , P_qi⁽ⁱ⁾ )^T is a qi-dimensional vector of parent variables of X_i in G.

Suppose that s observations x₁, x₂, . . . , x_s of random vector X are performed, and the observations of P_i are denoted p_1i, p_2i, . . . , p_si, where p_ji = (p_j1⁽ⁱ⁾, . . . , p_jqi⁽ⁱ⁾)^T is a qi-dimensional observation vector of parent genes. For example, assume X_ns as an n ´ s matrix, given that each of its components x_i (i = 1, . . . , n) is a vector of length s.

Therefore, X_ns = (x₁, . . . , x_n)^T = (x₍₁₎, . . . , x_(s)) = (x_ij)_{i=1,...,n ; j=1,...,s}, x_i = (x_i1, . . . , x_is), x_j = (x_1j , . . . , x_nj)^T, and x_i^T is the transposed vector of x_i. If X₁, for instance, has a parent vector P₁ = (X₄,X₅)^T, we can obtain p₁₁ = (x₁₄, x₁₅)^T , . . . , p_s1 = (x_s4, x_s5)^T.

Equation 1 still holds if we replace the P probabilities by density functions:

where x_ij is a particular value taken by X_ij. We then need to construct the conditional densities f_i (x_ij ½p_j⁽ⁱ⁾), where i = 1, . . . , n and j = 1, . . . , s. Following Tamada et al. (2003) and Imoto et al. (2002), we used nonparametric regression models for capturing the relations between the expression of gene i on sample j, x_ij and the expression levels of parents on j, p_ji = (p_j1⁽ⁱ⁾ , . . . , p_jqi⁽ⁱ⁾)^T through the equation:

where m_k(k = 1, . . . , qi) are smoothing functions from Â (the set of real numbers) to Â, and e_ij follows a normal distribution with mean 0 and variance s_i². For each function m_k, we assume that:

where {b_1k⁽ⁱ⁾, b_2k⁽ⁱ⁾, . . . , b_Mik⁽ⁱ⁾} is a preconceived set of “basis functions” (such as Fourier series, polynomial bases, B-splines, and others), the coefficients g_1k⁽ⁱ⁾, . . . , g_Mik⁽ⁱ⁾ are unknown parameters, and M_ik is the number of basis functions.

In our case, the basis functions are degree 3 B-splines; the coefficients g_1k⁽ⁱ⁾ , . . . , g_Mik⁽ⁱ⁾ are estimated from data, and M_ik = 20. There is a special algorithm to learn the parameters from data (Eilers and Marx, 1996). We do not detail how B-splines work as function approximation, in nonparametric regression; it is enough to know that they will “construct” a regulation function for each gene with the parent expression values. As we know the parents’ values p_ji = (p_j1⁽ⁱ⁾ , . . . , p_jqi⁽ⁱ⁾)^T, the B-splines generate an estimated expression for gene i, .

It remains to evaluate the quality of the B-spline approximation for a given set of parents P_i. This is done with the probability density function customized for a nonparametric regression model, using B-splines (Imoto and Konishi, 2000), when the ith gene has qi parents:

This result allows us to define the BN:

where q_G = (q₁^T , . . . , qn^T )^T is a parameter vector in graph G, and q_i is a parameter vector in the model of f_i, that is, q_i = (g_i^T, s_i²)^T or q_i = (µ_i, s_i²)^T.

Multiple search schemes and graph selection criterion

To learn the GN from data, we need a search algorithm that starts with an initial solution and traverses, using a criterion, across the space of possible networks to find the optimal network. To guide the search, we use scoring functions, also called evaluation criteria, which assign scores to each network according to its adherence to the observed data. Most of the evaluation criteria used to analyze a network use a common idea, the maximization (or minimization) of the posterior probability,

where p(G) is the prior probability of the graph G, and p(X | G) is the probability of the observations given the graph. The term p(X) (prior probability of the data) is constant and not related to p(G), and therefore will not be taken into account in the model’s evaluation. We use the Bayesian Information Criterion score, which has a term to measure the accuracy of the model in predicting data and another term to penalize the complexity of the model. For that reason, the more complex the graph structure, the worse its evaluation.

The space of solutions determined by all the acyclic directed graphs (a BN) grows exponentially with the number of nodes (Chickering, 1996). Under this scenario, we have implemented two heuristic search algorithms to seek for the optimal net, random hill-climbing and greedy hill-climbing (Ott et al., 2004). Briefly, both random hill-climbing and greedy hill-climbing make local changes in the network, through the following operations: 1) add a parent, 2) remove a parent, 3) reverse link direction, and 4) nothing. The difference between them is that the random approach chooses the modification at random, and only accepts it if it improves the overall score, while the greedy one evaluates the best possible move, applying the operators for the current network. Our implementation uses the greedy hill-climbing with many random restarts in order to better explore the search space.

The framework also has a variety of initialization options that define the starting graph for the learning procedure. We can begin the search with an empty graph, a random graph, a score-driven generated graph, and a graph that starts with interactions found by the partial correlation analysis (de La Fuente et al., 2004).

Next we present the results on applying the model to artificial and real biological expression data.

RESULTS AND DISCUSSION

Comparing search schemes

To compare the multiple search schemes available in the framework, we have used an artificial GN with 10 genes (Imoto et al., 2002). The artificial GN was sampled 100 times, each one being perturbed by a Gaussian noise (to mimic microarray experiments). This network, as well as the regulation functions that associate the genes, can be seen in Figure 3.

For random hill-climbing, we tested five different configurations for the initial network:

A. Score-driven: the parents of each gene are defined by the Bayesian information criterion score;

B. Random: the number of parents is sampled from a uniform distribution and they are chosen randomly among the genes;

C. Partial correlation: we input the initial graph with the links identified by the partial correlation analysis;

D. Partial correlation + score-driven: the links from the correlation analysis plus a score-driven initialization;

E. Partial correlation + random: the links from the correlation analysis, followed by a random initialization.

We performed 100 executions with each scheme and evaluated the distribution of the score for the learned networks (Figure 4). The lower the score, the better the learned GN. We see from this simulation that configurations A, B, and D produced solutions with scores in the interval [3.5 ´ 10³, 3.55 ´ 10³], which suggests that they had analogous performances. Configuration E generated solutions within two intervals. The best scores were achieved by scheme C, in which we only used the links found by the correlation analysis to start the graph, producing scores ranging from 3.31 ´ 10³ to 3.36 ´ 10³.

The same procedure was used with the greedy hill-climbing. After all simulations, we applied a t-test to decide which configuration was the best. The t-test points to scheme C, that is random hill-climbing and partial correlation initialization, as the best configuration to begin the search with.

Pheromone response pathway

As an application using a real biological system, we have applied the BN model to the Saccharomyces cerevisiae pheromone response pathway. This signaling pathway is responsible for a series of complex physiological changes in preparation for mating in haploid yeast cells, including changes in the expression of about 200 genes, arrest in the G1 phase of the cell cycle and membrane fusion of mating partners (Bardwell, 2005).

The transmission of the signal starts when the pheromone peptide is recognized by cell surface receptors, leading to a signaling cascade that ultimately activates genes needed for mating. STE2 acts as the specific MATa cell pheromone receptor, while STE3 is the receptor for the MATa cells (Lewin, 2000). After pheromone recognition, the G-protein coupled to the pheromone receptor (a trimeric complex formed by GPA1, STE4 and STE18) recruits its effectors in order to propagate the signal to an MAPK signaling cascade that finally activates the STE12 transcription factor, the major player that will command the DNA transcription. The G-protein effectors are the STE5/STE11 and FAR1/CDC24 complexes as well as the STE20 protein kinase (Bardwell, 2005). A schematic view of the mating pathway is shown in Figure 5.

We generated GNs for two microarray datasets. As the first dataset, we used a subset of 12 genes from the cell cycle expression data generated by Spellman et al. (1998). These data have 18 samples (time series experiment). Ten of the genes, STE2 or STE3, STE18, STE20, STE5, STE11, STE7, FUS3, GPA1, FAR1, and STE12, act directly on the mating pathway, while AGA1 and FUS1 are required for cell fusion (Bardwell, 2005). Figure 6 shows the GNs (GN-1 and GN-2) inferred from this dataset, using STE2 and STE3, respectively. The nodes were colored depending on the known function of the gene (see Table 1 for the mapping). For this reason, it is desirable that genes with related functions occur near each other in the graph.

For each experiment, we ran 600 executions (i.e., 600 learned networks) and the output network was defined by a “voting criterion”. In conformity with this criterion, we chose a threshold, and only the links that appeared above this threshold across all executions composed the resultant network. In our experiments, we used a threshold value of 55% of the trials. The program run took about 20 h in a Pentium 4 with 2.4 GHz and 1 GB RAM.

From Figure 6, we can observe that in both networks (GN-1 and GN-2), STE12 has an important position, because it has the highest number of children. Thus, the recovered networks suggest that this gene has an important role in transcriptional regulation. Indeed, as indicated by the biological literature, STE12 is the most important transcription factor that is activated by the pheromone response pathway (Zeitlinger et al., 2003), and it influences the transcription of a number of other genes.

For GN-1, the network successfully predicts known regulatory interactions between STE12 and genes STE2, FUS3 and FAR1, which are positively acting components of the pathway (Bardwell, 2005). It did not capture known interactions with AGA1, a gene involved in cell fusion (White and Rose, 2001), nor were interactions found with GPA1, a negative regulator of the pathway (Bardwell, 2005). However, GN-1 reports three interesting relations of STE12 with STE20, STE5 and STE11. STE5 is an adaptor protein that binds to and activates STE11, the upstream kinase on the MAPK cascade. STE20 is another protein in the neighborhood that is another activator for STE11. GN-1 indicates some kind of regulation (activation or inhibition) of these genes by STE12, which we did not find evidenced in a specialized database (SCPD, 2005). We hypothesize that STE12 is in someway upregulating these genes, because MAPK targets will act towards maintaining the transcription of STE12 target genes (e.g., inhibition of Dig1 and Dig2, which are inhibitors of STE12, thus benefiting the maintainability of mating genes transcription). The diagram in Figure 7 provides an illustration of this scenario. The links around STE12, inferred by GN-1, are highlighted (red) and the links in black depict interactions at protein level. This diagram also allows one to note that GNs encode protein-DNA interactions besides only protein-protein relationships, as the case of the signaling pathway in Figure 5.

GN-1 also indicates that FAR1 expression is controlled by STE12 and GPA1. One of the functions of FAR1 is to promote G1 cell-cycle arrest. Given the fact that GPA1 is a negative regulator of the pheromone response (see above), this link (GPA1, FAR1) is a feasible hypothesis, and it remains to be determined by which mechanism GPA1 down-regulates FAR1.

STE18 is the gene with the second-most regulatory connections in GN-1 (Figure 6). STE18 expresses one of the subunits of the G-protein complex, bound to the pheromone receptor. GN-1 shows that mRNA expression of STE18 affects the regulation of STE7 and STE11 (both protein kinases), as well as STE12, the main transcription factor.

GN-2, which contains STE3 instead of STE2, has also revealed STE12 as the most important gene in the transcriptional regulation of the pathway. GN-2 has found the same interactions around STE12 in the first experiment, with STE3 (in place of STE2), FUS3, FAR1, STE20, STE5 (inverted), and STE11, previous justified. It has also found two interactions evidenced in the literature (Bardwell, 2005) not found by GN-1: a direct interaction (STE12, GPA1) and a second-level interaction (STE12, AGA1).

The second dataset (Hartemink et al., 2002) had 12 genes that are directly or closely involved with the mating pathway, and it was used to infer a third network GN-3. The expression data consisted of 20 observations of the two types of haploid yeast cells (a and a), under several conditions, including exposure to different nutritive media, and exposure to various types of stress, such as heat, oxidative species, excessive acidity, and excessive alkalinity.

GN-3 is presented in Figure 8. STE12 again is the most connected gene, but here some important regulatory connections, such as (STE12, FUS1) and (STE12, STE2) were not discovered, even though several others were identified, including (STE12, FAR1) and (STE12, GPA1). This example shows that even for the same subset of genes, the design of the experiment is a key issue for gene network inference. The first dataset was generated in a time series fashion, whereas in the second one, there were changes in the conditions, but not in time. To properly analyze the dynamics of expression, it is better to use time series microarray data (Bar-Joseph, 2004).

A last important comment can be made about the interaction (STE12, FAR1), which appears in all three networks. According to information from the SCPD yeast promoter database (SCPD, 2005), FAR1 has only one transcription factor, called MCM1, which is found by a motif tool. Nonetheless, our networks hint that STE12 also regulates FAR1 by some mechanism, which could be through the binding of STE12 to an upstream region of that gene, as depicted in Figure 9. Indeed, a chromatin immuno-precipitation assay revealed that STE12 also binds to the upstream region of FAR1 (Ren et al., 2000), although it does not exhibit a motif for STE12. This fact best exemplifies what kind of knowledge about transcriptional control one could extract from GNs, and in that manner rationally design wet experiments.

CONCLUDING REMARKS

We have introduced a BN model for GNs, and we have tested it with both artificial and real biological networks. We analyzed the yeast pheromone response pathway, and we demonstrated the usefulness of GNs as a computational approach for the analysis of transcriptional regulation. In summary, a GN can be used, among other things, to i) define transcriptional factors (activators and inhibitors) for a target gene and ii) find co-regulated genes. The intention of the efforts for developing both theories and software for network analysis is that these networks could provide useful clues about biological systems, thus helping with the design and refinement of wet experiments.

The BN model is suitable for small networks. A learning scheme that scales-up to a large number of variables should be investigated, and is a future goal. Nowadays, the finding of an efficient reconstruction method with no constraints in the number of nodes using BN is a cutting-edge problem (Bar-Joseph, 2004).

We are aware of the limitations of gene networks as a way to understand the behavior of a biological system, in terms of phenotype. These occur because 1) there is a low correlation between expression level and protein level (Ideker et al., 2001), and 2) much of the cellular regulation occurs post-translationally, and genomic-scale technologies to measure protein levels are still at beginning stages of development (Rice and Stolovitzky, 2004). For this reason, and to attenuate this weakness, a trend in this field of network analysis is information fusion. According to this concept, the algorithms to learn the networks should incorporate other sources of biological data, such as location data, sequence data and protein-protein interactions. Some researchers have already proposed models in this direction (e.g., Yamanashi et al., 2004; Bernard and Hartemink, 2005) and this will become an active research topic in the next years.

ACKNOWLEDGMENTS

We thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for financial support. The authors also thank A.J. Hartemink for providing the yeast microarray data.

REFERENCES

Bar-Joseph Z (2004). Analyzing time series gene expression data. Bioinformatics 20: 2493-2503.

Bardwell L (2005). A walk-through of the yeast mating pheromone response pathway. Peptides 26: 339-350.

Bastos G and Guimarães KS (2005). A simpler Bayesian network model for genetic regulatory network inference. Proceeding of International Joint on Neural Networks, Montreal, Canada.

Bernard A and Hartemink AJ (2005). Informative structure priors: Joint learning of dynamic regulatory networks from multiple types of data. Pac. Symp. Biocomput., Hawaii, USA, pp. 459-470.

Chickering D (1996). Learning Bayesian networks is NP-complete. In: Learning from data: AI and statistics V (Fisher D and Lenz HJ, eds.). Springer-Verlag, New York, NY, USA, pp. 121-130.

de Jong H (2002). Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol. 9: 67-103.

de la Fuente A, Bing N, Hoeschele I and Mendes P (2004). Discovery of meaningful associations in genomic data using partial correlation coefficients. Bioinformatics 20: 3565-3574.

Eilers PHC and Marx BD (1996). Flexible smoothing with B-splines and penalties. Stat. Sci. 11: 89-121.

Friedman N (2004). Inferring cellular networks using probabilistic graphical models. Science 303: 799-805.

Friedman N, Linial M, Nachman I and Pe’er D (2000). Using Bayesian networks to analyze expression data. J. Comput. Biol. 7: 601-620.

Hartemink AJ, Gifford DK, Jaakkola TS and Young RA (2002). Combining location and expression data for principled discovery of genetic regulatory network models. Pac. Symp. Biocomput. pp. 437-449.

Ideker T, Thorsson V, Ranish JA, Christmas R, et al. (2001). Integrated genomic and proteomic analyses of a systematically perturbed metabolic network. Science 292: 929-934.

Imoto S and Konishi S (2000). B-spline nonparametric regression models and information criteria. Proceedings of 2nd Int. Symp. on Frontiers of Time Series Model, pp. 240-241.

Imoto S, Goto T and Miyano S (2002). Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression. Pac. Symp. Biocomput., pp. 175-186.

Lewin B (2000). Genes VII. Oxford University Press, New York, NY, USA.

Nachman I, Regev A and Friedman N (2004). Inferring quantitative models of regulatory networks from expression data. Bioinformatics 20 (Suppl 1): I-248-I-256.

Ott S, Imoto S and Miyano S (2004). Finding optimal models for small gene networks. Pac. Symp. Biocomput., Hawaii, USA, pp. 557-567.

Ren B, Robert F, Wyrick JJ, Aparicio O, et al. (2000). Genome-wide location and function of DNA binding proteins. Science 290: 2306-2309.

Rice JJ and Stolovitzky G (2004). Making the most of it: pathway reconstruction and integrative simulation using the data at hand. Biosilico 2: 70-77.

Roweis S and Ghahramani Z (1999). A unifying review of linear Gaussian models. Neural Comput. 11: 305-345.

Schafer J and Strimmer K (2005). An empirical Bayes approach to inferring large-scale gene association networks. Bioinformatics 21: 754-764.

SCPD (2005). The promoter database of Saccharomyces cerevisiae. Available at: http://rulai.cshl.edu/SCPD/. Accessed August, 2005.

Spellman PT, Sherlock G, Zhang MQ, Iyer VR, et al. (1998). Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell 9: 3273-3297.

STKE (2005). Science Signal Transduction Knowledge Environment (STKE). Available at: http://stke.sciencemag.org/cgi/content/full/sci;306/5701/1508/FIG1. Accessed August, 2005.

Tamada Y, Kim S, Bannai H, Imoto S, et al. (2003). Estimating gene networks from gene expression data by combining Bayesian network model with promoter element detection. Bioinformatics 19 (Suppl 2): II-227-II-236.

van Someren EP, Wessels LF, Backer E and Reinders MJ (2002). Genetic network modeling. Pharmacogenomics 3: 507-525.

White JM and Rose MD (2001). Yeast mating: getting close to membrane merger. Curr. Biol. 11: R16-R20.

Xia Y, Yu H, Jansen R, Seringhaus M, et al. (2004). Analyzing cellular biochemistry in terms of molecular networks. Annu. Rev. Biochem. 73: 1051-1087.

Yamanishi Y, Vert JP and Kanehisa M (2004). Protein network inference from multiple genomic data: a supervised approach. Bioinformatics 20 (Suppl 1): I-363-I-370.

Zeitlinger J, Simon I, Harbison CT, Hannett NM, et al. (2003). Program-specific distribution of a transcription factor dependent on partner transcription factor and MAPK signaling. Cell 113: 395-404.