Bioinformatics Advance Access originally published online on March 31, 2005
Bioinformatics 2005 21(12):2883-2890; doi:10.1093/bioinformatics/bti415
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
A stochastic differential equation model for quantifying transcriptional regulatory network in Saccharomyces cerevisiae
1Division of Molecular and Genomic Medicine, National Health Research Institutes Zhunan Town, Miaoli County 350, Taiwan
2Bioinformatic Laboratory, Department of Computer Science and Information Engineering, National Taiwan University Taipei 106, Taiwan
*To whom correspondence should be addressed.
 |
Abstract |
|---|
 TOP Abstract 1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Motivation: The explosion of microarray studies has promised to shed light on the temporal expression patterns of thousands of genes simultaneously. However, available methods are far from adequate in efficiently extracting useful information to aid in a greater understanding of transcriptional regulatory network. Biological systems have been modeled as dynamic systems for a long history, such as genetic networks and cell regulatory network. This study evaluated if the stochastic differential equation (SDE), which is prominent for modeling dynamic diffusion process originating from the irregular Brownian motion, can be applied in modeling the transcriptional regulatory network in Saccharomyces cerevisiae.
Results: To model the time-continuous gene-expression datasets, a model of SDE is applied to depict irregular patterns. Our goal is to fit a generalized linear model by combining putative regulators to estimate the transcriptional pattern of a target gene. Goodness-of-fit is evaluated by log-likelihood and Akaike Information Criterion. Moreover, estimations of the contribution of regulators and inference of transcriptional pattern are implemented by statistical approaches. Our SDE model is basic but the test results agree well with the observed dynamic expression patterns. It implies that advanced SDE model might be perfectly suited to portray transcriptional regulatory networks.
Availability: The R code is available on request.
Contact: cykao{at}csie.ntu.edu.tw
Supplementary information: http://www.csie.ntu.edu.tw/~b89x035/yeast/
|
1 INTRODUCTION |
|---|
|
TOPAbstract1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
One of the major challenges in biology is the construction of the complicated mechanism that regulates gene expression. Living cells are controlled by gene-expression programs involving specific promoter sequences and various transcription regulators. These transcriptional programs are modified as cells progressing through the cell-cycle. There are a growing number of examples that the cell-cycle regulated genes are typically those involved in cell-cycle specific processes, which are expressed at peak levels at the time that they are needed for particular functions (Spellman et al., 1998; Lee et al., 2002). However, the molecular mechanisms governing gene-expression regulation remained largely unknown.
Recently, the development of microarray, which permits us to monitor transcriptomes on a genome-wide scale, has dramatically expedited comprehensive understanding of the gene-expression profiles, e.g. how the transcription profiles for genes vary across the cell cycle. However, it is generally believed that analyzing a vast number of microarray datasets does not necessarily provide comprehensive understanding of the role of a given gene in a specific biological process. The unmet needs include efficient translation of existing raw microarray data into a workable model to facilitate greater comprehension of the transcriptional regulatory network. In this study, we present a novel modeling, a stochastic differential equation (SDE) model, to allow quantitative depiction of time-course transcriptional regulation. The model also implements statistical approaches to estimate the contribution of the transcriptional regulators.
In this study, cell-cycle data of budding yeast Saccharomyces cerevisiae (Spellman et al., 1998) was used as a test set. A list of target genes was collected, including 800 cell-cycle regulated genes from Spellman et al. (1998) and 2119 selected genes from Chen et al. (2004). In addition, a candidate pool of potential upstream regulators of target genes consisting of transcription factors, cell cycle control factors and DNA-binding transcriptional regulators as described in literatures (Spellman et al., 1998; Chen et al., 2004; Harbison et al., 2004) was selected. The goal was to fit a generalized linear model using a set of putative regulators to estimate the transcription pattern of a specific target gene. Goodness-of-fit was evaluated by log-likelihood and Akaike Information Criterion (AIC) (Akaike, 1974). Then, the regulatory abilities of selected regulators were estimated by maximum likelihood, and the estimated transcription rate of the target gene was inferred by combining the regulatory functions of selected regulators. Finally, the estimated expression level of the target gene was derived by least squares.
The results from test data fitting curves agreed well with the observed dynamic expression patterns. The first 10 of the best fittings and worst fittings and selected figures of the fitting curves are listed below for comparison. Taken together, test results indicate that our SDE model is good at portraying transcriptional regulatory networks.
|
2 SYSTEMS AND METHODS |
|---|
|
TOPAbstract1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
To construct a systematic study on periodically expressed data, an SDE model was applied to describe the transcriptional process of a target gene. For each target gene, the model associates its transcription rate with the regulatory functions of possible regulators. Then the model iteratively selects regulators from the candidate pool so that the combination of selected regulators may result in the highest log-likelihood; the selection procedure ends when AIC begins to increase. At this stage, the estimation of the contribution and regulatory abilities of selected regulators was carried out.
The calculation and statistical analysis conducted in this paper were implemented by R 2.0.1 (R Development Core Team, http://www.r-project.org/).
2.1 Dynamic stochastic transcription model
In this paper, the modeling of the temporally dynamic gene-transcription process is based on the theory of SDE (Chung and Williams, 1990). The SDE is prominent for modeling dynamic diffusion process that originates from the irregular Brownian motion of pollen particles suspended in water. SDE is widely used for modeling irregular motion, variability or uncertainty due to time series, and have been used to provide powerful solutions in finance, neuroscience, computational cell biology and physics (Tan, 2002). To model the time-continuous gene-expressions datasets, we let Nt denote the exact amount of target gene mRNA at time t. From time t to t + 
t, the dynamic transcription and degradation process is modeled by
 |
is the degradation rate, and 
t,t is the noise or random error captured by normal distribution N(0,
2 t) with zero mean and variance in proportion to t.
Suppose that Wt 
N(0, t), then t, t = Wt and
 |
t 
0, we have the SDE as follows,
 |
 | (1) |
N
t is the quadratic variation of Nt with 
.
2.2 Regulatory functions
At time t, the transcription rate gt of the target gene depends on the n regulators through regulatory functions. Let Xit be the expression level of the i-th regulator at time t. And the regulatory function of a specific regulator i is described as a sigmoid function (von Seggern, 1993):
 | (2) |
i are the mean and deviation of regulator i, respectively. Since the expression level is measured at m time points, there are m samples for time t = t1, t2,...,tm. We estimate µi and i by sample mean 
and sample standard deviation 
for regulator i.
 | (3) |
 |
, and ci is the contribution of regulator i to the transcriptional regulatory network.
2.3 Statistical method
In order to connect the target gene with its n upstream regulators, the model associates the stochastic transcription rate of target gene with the regulatory functions of possible regulators. In statistical sense, the regulators are considered as predictors and the target gene is regarded as response variable. The statistical approach is to fit a linear model using a set of regulatory functions of regulators as the inputs to estimate the dynamic transcription rate of a target gene as the output. It also estimates the contribution and regulatory abilities of selected regulators. For a small t, the previous equation can be approximated by a discrete form as follows:
 |
 |
is distributed as N(0, 1), which represents the source of uncertainty of random error. Thereafter, for a specific target gene regulated by n regulators, m samples are collected at time t = t1, t2,...,tm, we have
 |
 | (4) |
 | (5) |
 |
 |
For simplicity, the above model is written in the matrix notion as
 |
are the parameters to be estimated. Since Ztj are i.i.d. random variables with standard normal distribution, the log-likelihood function of Y is
 |
and 
, since at least one second-order partial derivative is negative, 
, the maximum-likelihood estimators (MLEs) are derived as
 |
 |
in our model, since it is already contained in the c0 term. The selection procedure is used to build a model that might best fit the dataset and is carried out by the AIC (Akaike, 1974) given by
 |
2.4 Computational algorithm
Our method aims to select a set of possible regulators to estimate the dynamic transcription pattern of a target gene. The method searches for possible regulators from the candidate pool and uses an AIC strategy to iteratively add regulators resulting in the highest log-likelihood; the selection ends when an increasing AIC appears. Consequently, the regulatory ability is evaluated by the ML approach and the transcription pattern is derived by least squares (LS) method.
The selection of regulators from the candidate pool is conducted by the AIC forward selection procedure (Weisberg, 1985). In the forward selection (FS), regulators are added at each step and the FS procedure works as follows.
The simple regression model with the single predictor (regulator) that has the biggest log-likelihood with the response (target gene) is used. Then, an additional regulator is added if it meets the criterion:
- It will decrease the AIC more than any other single regulator or has the largest log-likelihood of any of the regulators that are not already in the model.
- Stop when the AIC starts to increase.
- Stop if the number of regulators in model would be larger than some predetermined number, say 10, in our method.
Meanwhile, incomplete data is handled by the fill-in method. Since this data are in time-series manner, it is reasonable to estimate a missing value by the interpolation of its adjacent time-point observations. In addition, if the missing value is at the first or the last time point, an extrapolation estimate is applied.
Thus, the overall algorithm is described as follows:
- Step 1. Handle the incomplete data by fill-in method and calculate the response Y by Equation (4) for all genes in the dataset.
- Step 2. Calculate the regulatory functions of regulators by Equations (2) and (5) for all genes in the candidate pool.
- Step 3. Conduct the AIC–FS approach to select a set of possible regulators to explain the dynamic transcription pattern of a target gene, as well as the parameters estimated by maximum-likelihood method.
- Step 4. Derive LS method to estimate the expression level of the target gene through the selected regulators. The formulae are outlined in the supporting material.
- Step 2. Calculate the regulatory functions of regulators by Equations (2) and (5) for all genes in the candidate pool.
|
3 IMPLEMENTATION |
|---|
|
TOPAbstract1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
3.1 Yeast cell cycle microarray data
To evaluate our model on biological time-course datasets, we applied it to the eukaryote S.cerevisiae (Spellman et al., 1998) data. The cell-cycle process has been studied in great detail in the budding yeast S.cerevisiae, and the dataset contains gene-expression measurements of the mRNA levels of 6178 S.cerevisiae open reading frames (ORFs) at 18 time-points under
factor-based synchronization method. Spellman et al. (1998) identified 800 genes whose expression varied over the different cell-cycle stages by meeting an objective minimum criterion for cell-cycle regulation. The data were transformed into log2 ratio scale, and normalized so that the average log2(ratio) over the course of the experiment was equal to 0. Although there are 800 genes involved in the cell-cycle process of the budding yeast, the number of key regulators that are responsible for the control and regulation of this complex process is much smaller (Li et al., 2004). We performed a genome-wide dynamic analysis to investigate how yeast regulators regulate target transcripts. A list of target genes was constructed by joining 800 regulated genes from Spellman et al. (1998) with 2119 selected genes from Chen et al. (2004) Then, a candidate pool of potential upstream regulators of target genes was constructed by joining transcription factors, cell-cycle control factors and DNA-binding transcriptional regulators described in literature (Spellman et al., 1998; Chen et al., 2004; Harbison et al., 2004). The selection procedure was carried out by the AIC strategy. The lower the AIC score was, the more the pattern displayed by these selected regulators resembled the target gene. Although our test set was limited to a candidate pool of known regulators, the model could be applied to the whole genome pool as well.
3.2 Inference from dynamic stochastic model
Applying estimated coefficients to the selected regulators, the transcriptional regulatory ability was inferred by a combination of regulatory functions, and the transcription pattern of a specific target gene was derived by LS method. The results of our method for 20 target genes selected by Chen et al. (2004) are listed in Table 1. The predicted values, 
, agreed well with the observed expression patterns of six target genes (Table 1 and Fig. 1). Table 2 contains the first 10 best fittings and worst fittings, which were determined on the log-likelihood of Y, the transcription rate minus degradation rate, in our implementation, and selected graphs of fitting curves are presented for comparison in Figure 2. Since the dynamic transcription model here was carried out by the basic theory of SDE, it was not supposed to be powerful enough for identifying all real regulators. However, while we carefully discerned the worst 10 fitting curves, they still captured much of the dynamic patterns of temporal transcription levels (Fig. 2). It was implied that advanced SDE models might be well suited to characterize transcriptional regulatory networks. A summary statistics based on the SDE model is shown in Table 3. In a total of 2341 target genes, the number of regulators ranges mainly from one to six. Very few target genes (5, 0.21%) have 10 regulators. The average number of regulators per gene is 3.0829, and the average percent of activators of the corresponding total regulators are indicated in the fourth column. The percentage composition of the regulators (activators versus repressors) is close to 0.5 for each number of regulators, although the number of activators seems a little less than repressors when total regulators are fewer. In addition, we tested the SDE method on some Spellman's data using fewer experimental points. The sum of square errors (SSE) was still small when using only as little as four time points. The results are available on the website.
|
|
|
|
|
To sum up, our study provides not only a mathematical and statistical method to capture time-course expression information, but also a combination of regulatory functions to characterize regulatory abilities. Overall, the estimated
is pretty close to the observed Xt. In addition, the regulatory parameter cj represents the contribution of regulatory function of regulator i to the target gene, with activation for positive estimation and repression for negative estimation. The test data, complete tables of specific target genes with their selected regulators and fitted curves graphs are available on the website. |
4 DISCUSSION |
|---|
|
TOPAbstract1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Biological systems have been modeled as dynamic systems for a long time, such as immune systems, genetic networks and cell regulatory network (Li et al., 2004). Gene-expression profile data are network interactions and have been analyzed through a variety of clustering methods in an attempt to group genes with similar patterns (Eisen et al., 1998; Tavazoie et al., 1999). Another widely-used application is to organize co-expressed genes by promoter sequence motif (Bussemaker et al., 2001). Their approaches deliberately depict the expression profiles of transcription network. A genome-wide location analysis of gene expression to identify transcription factors binding to promoters (Iyer et al., 2001; Lee et al., 2002) could provide us more information to construct transcriptional regulatory networks. More recently, a few studies in discerning and dissecting periodically expressed cell cycle yeast data are the fuzzy logic approach (Woolf and Wang, 2000), the smooth surface response (SRS) algorithm (Xu et al., 2002), the quantitative characterizing transcriptional regulatory network by Chen et al. (2004) and so on. The fuzzy logic approach and SRS algorithm were proposed to generate a connected network and to construct a model to find triplets of activators, repressors and target genes. Their algorithms well construct parts of gene regulatory machinery in yeast cell cycle. In this paper, we address this question from a dynamic systems point of view. Our approach provides a novel study in characterizing the time series expression level based on the basic Itô rule of SDE theory, as well as identifying possible regulators by common statistical approach. The model explicitly deals with the variable time problem by addressing the term of squared root of time interval in transcription model and sigmoid function for all genes. Estimating by ML and LS, this model can unambiguously portray the regulatory abilities of temporally expressed target genes. Our model not only studies the regulatory relationship between regulators and target genes, but also quantifies the regulatory abilities to the specific target genes. It concerns a method of multi-regulators, but not a triplet model of activator, repressor and target gene.
Our method has several advantages over other implementations. An important feature in the SDE model is that dynamics arise as a natural consequence of time course in the duration of the cell cycle process. The SDE model takes time course into account for addressing continuous-time random error in Brownian motion. The fitted curves adequately depict differently shaped expression patterns while keeping the model parameters as few as possible. Applying more advanced results of SDEs theory, our algorithm can be implemented in more complex dynamic biological systems.
Our approach has, inevitably, certain disadvantages. First, we did not concern interactions between regulators in our model of combinational regulatory functions; undeniably, the real-world regulatory mechanism is much more complicated than as assumed in our model. Second, the AIC forward selection might not be good enough for model selection. A stepwise selection or other criteria for variable selection, such as Mallows Cp (Weisberg, 1985) for selecting a subset, could be applied later. Third, if there are true regulators that have not been identified in literature, they are not included in the candidate pool. As a result, they would never be identified as regulators by our algorithm. All possible genome-wide regression could be considered in the future, though this might be inefficient and impractical. In general, there is no final or best model, only a group of possible models that are all judged as nearly, equally useful. In addition, a basic limitation of this study is the data of mRNA expression level might not exactly reflect the deliberate regulatory mechanism. Modifications, corrections or filtering might be necessary. All of these issues are subjected to further study. Nevertheless, the results indicate that our SDE model can capture the behavior of transcriptional regulatory throughput quite accurately. As noted above, inference under our model for cell-cycle expression is promising and one of the goals in our future studies of cell-cycle gene expression is to construct the sophisticated cell regulatory networks.
|
Acknowledgments |
|---|
We are grateful for the support of the Department of Computer Science and Information Engineering, National Taiwan University, Taipei, Taiwan. This work was supported by grants from the National Health Research Institutes and National Science Council (NRPGM: NSC93–3112-B-400–001) to C.F.H. We thank Dr. Chen-Hsiung Chan of NTU for his critical reading of the manuscript and fruitful discussion.
|
Footnotes |
|---|
Current address: Institute of Statistical Science, Academia Sinica, 128, Academia Rd. Sec. 2, Taipei 115, Taiwan 
Received on February 4, 2005; revised on March 10, 2005; accepted on March 29, 2005
|
REFERENCES |
|---|
|
TOPAbstract1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Akaike, H. (1974) A new look at the statistical model identification. IEEE Trans. Automat. Control, AC-19, 716–723[CrossRef].
Bussemaker, H.J., et al. (2001) Regulatory element detection using correlation with expression. Nat. Genet., 27, 167–171[CrossRef][Web of Science][Medline].
Casella, G. and Berger, R.L. Statistical Inference, (1990) , Belmont Chapter 7 Wadsworth, Inc.
Chen, H.-C., et al. (2004) Quantitative characterization of the transcriptional regulatory network in the yeast cell cycle. Bioinformatics, 20, 1914–1927
Chung, KL and Williams, R.J. Introduction to Stochastic Integration, (1990) 2nd edn. , Boston Birkhäuser.
Eisen, M.B., et al. (1998) Cluster analysis and display of genome-wide expression patterns. Proc. Natl Acad. Sci. USA, 95, 14863–14868
Harbison, C.T., et al. (2004) Transcriptional regulatory code of a eukaryotic genome. Nature, 431, 99–104[CrossRef][Medline].
Iyer, V.R., et al. (2001) Genomic binding sites of the yeast cell-cycle transcription factors SBF and MBF. Nature, 409, 533–538[CrossRef][Medline].
Karatzas, I. and Shreve, S.E. Brownian Motion and Stochastic Calculus, (1991) 2nd edn. , New York Springer-Verlag.
Lee, T.I., et al. (2002) Transcriptional regulatory networks in Saccharomyces cerevisiae. Science, 298, 799–804
Li, F., et al. (2004) The yeast cell-cycle network is robustly designed. Proc. Natl Acad. Sci. USA, 101, 4781–4786
Li, W. and Nyholt, D.R. (2001) Marker selection by Akaike information criterion and Bayesian information criterion. Genet. Epidemiol., 21, Suppl. 1, S272–S277.
Spellman, P.T., et al. (1998) Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell, 9, 3273–3297
Tan, W.-Y. Stochastic Models with Applications to Genetics, Cancers, AIDS and Other Biomedical Systems, (2002) World Scientific Publishing Co. Ltd.
Tavazoie, S., et al. (1999) Systematic determination of genetic network architecture. Nat. Genet., 22, 281–285[CrossRef][Web of Science][Medline].
von Seggern, D. CRC Standard Curves and Surfaces, (1993) , Boca Raton CRC Press, pp. 124.
Weisberg, S. Applied Linear Regression, (1985) chapter 8 John Wiley and Sons Inc.
Woolf, P.J. and Wang, Y. (2000) A fuzzy logic approach to analyzing gene expression data. Physiol. Genomics, 3, 9–15
Xu, H., et al. (2002) A smooth response surface algorithm for constructing a gene regulatory network. Physiol. Genomics, 11, 11–20
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  D. Nam, S. H. Yoon, and J. F. Kim Ensemble learning of genetic networks from time-series expression data Bioinformatics, December 1, 2007; 23(23): 3225 - 3231. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 T. T. Vu and J. Vohradsky Nonlinear differential equation model for quantification of transcriptional regulation applied to microarray data of Saccharomyces cerevisiae Nucleic Acids Res., January 12, 2007; 35(1): 279 - 287. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
X.-w. Chen, G. Anantha, and X. Wang An effective structure learning method for constructing gene networks Bioinformatics, June 1, 2006; 22(11): 1367 - 1374. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||