PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. Nov 2012; 8(11): e1002772.
Published online Nov 8, 2012. doi:  10.1371/journal.pcbi.1002772
PMCID: PMC3493476
Stochastic Modeling of Expression Kinetics Identifies Messenger Half-Lives and Reveals Sequential Waves of Co-ordinated Transcription and Decay
Filippo Cacace,#1 Paola Paci,#1,2* Valerio Cusimano,1 Alfredo Germani,3 and Lorenzo Farina4
1Università Campus Biomedico, Rome, Italy
2CNR-Institute of Systems Analysis and Computer Science (IASI), BioMathLab, Rome, Italy
3Dipartimento di Ingegneria e Scienze dell'Informazione e Matematica (DISIM), Università de L'Aquila, L'Aquila, Italy
4Dipartimento di Informatica e Sistemistica, Sapienza Università di Roma, Rome, Italy
Teresa M. Przytycka, Editor
National Center for Biotechnology Information (NCBI), United States of America
#Contributed equally.
* E-mail: paola.paci/at/iasi.cnr.it
The authors have declared that no competing interests exist.
Conceived and designed the experiments: LF AG. Performed the experiments: FC PP VC. Analyzed the data: PP LF. Wrote the paper: FC PP AG LF. Made the figures: LF.
Received May 5, 2012; Accepted September 23, 2012.
The transcriptome in a cell is finely regulated by a large number of molecular mechanisms able to control the balance between mRNA production and degradation. Recent experimental findings have evidenced that fine and specific regulation of degradation is needed for proper orchestration of a global cell response to environmental conditions. We developed a computational technique based on stochastic modeling, to infer condition-specific individual mRNA half-lives directly from gene expression time-courses. Predictions from our method were validated by experimentally measured mRNA decay rates during the intraerythrocytic developmental cycle of Plasmodium falciparum. We then applied our methodology to publicly available data on the reproductive and metabolic cycle of budding yeast. Strikingly, our analysis revealed, in all cases, the presence of periodic changes in decay rates of sequentially induced genes and co-ordination strategies between transcription and degradation, thus suggesting a general principle for the proper coordination of transcription and degradation machinery in response to internal and/or external stimuli.
The amount of a given transcript in a cell is determined by a fine tuned balance of production and degradation in a complex regulatory network. Regulation of transcription controls when transcription occurs and how much mRNA is created, whereas regulation of degradation controls the rate at which messengers are destroyed. The latter mechanism has recently gained attention due to the increasing evidence of its key role in the overall co-ordination of gene expression. A long lifetime of mRNA enables a cell to produce more proteins from that mRNA. By contrast, a short lifetime rapidly alters protein levels in response to changing needs. Measuring mRNA stability is a complex and expensive experiment and, given the condition-specific response of the degradation pathway, it would be desirable to take advantage of the large variety of expression experiments stored in public databases. To this end, we developed a stochastic model to infer each specific mRNA lifetime from gene expression data. Predictions were validated using malaria data. We then applied our methodology to the reproductive and metabolic cycle of budding yeast and found, in all cases, the presence of a general principle for the proper coordination of transcription and degradation machinery.
Appropriate and timely changes in gene expression are essential for cell life. The transcriptome is finely regulated by a large number of molecular mechanisms able to adjust the balance between mRNA production and degradation. Every aspect of transcript life is subject to elaborate control but, traditionally, the focus of the research has been on transcriptional regulation [1]. However, whereas mRNA abundance results from the dynamic interplay between transcription and degradation, the speed by which cells can adjust their mRNA levels is critically dependent on the rate of mRNA turnover [2]. As a result, small changes in mRNA stability may dramatically drive rapid variations of transcript abundance. Efforts to understand the underlying principles of mRNA decay and transcription co-ordination are very important since the balance between transcription and decay influences most, if not all, the cell responses to endogenous and exogenous signals [3].
The current widespread interest in this topic has been fostered by the finding of specific regulatory mechanisms of mRNA stability such as, for example, RNA binding proteins [4], [5] and small RNAs [6]. Regulation of transcript stability cannot be considered a simple “disposal system” but a sophisticated tool for the proper orchestration of the global cell response to internal and external stimuli [6]. Remarkably, a key role of mRNA stability has been reported in cancer, inflammatory diseases and Alzheimer's [7]. In recent years there has been a surge in empirical studies that measured, on a genome-wide scale in a variety of environmental conditions, messenger half-lives of many organisms, including plants [8] mammals [9] and fungi [2]. The discovery of such new regulatory layer has clarified that, in order to obtain a clear picture of the underlying regulatory machinery, it is necessary to complement the traditional time-course experiment measuring the cell transcriptional response under certain conditions far from steady state with decay rates data under the same condition [10].
Experimental procedures for the evaluation of mRNA decay rates are based on measuring gene expression upon inhibition of transcription [11][13] or on pulse-chase RNA labeling protocols [2], [13][15]. Such protocols are very critical (see Figure 1 for a comparison among different studies summarized in Table 1), since, for instance, transcriptional shut-off blocks growth and has a profound effect on cellular physiology, as well as on mRNA metabolism [2]. In fact, Wang et al. [11] and Grigull et al. [12] datasets show a low value of the Pearson correlation (An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e001.jpg), and no correlation at all can be found (An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e002.jpg) between Munchel et al. [2] and Wang et al. [11] datasets (see Figure 1A and Figure 1B respectively). Despite the same experimental conditions (asynchronous growth), the two half-life independent measurements obtained by Wang et al. and Munchel et al. are uncorrelated, probably due to differences in the shut off protocol (pulse chase for [2] and thermal inactivation for [11]), whereas Grigull et al. and Wang et al. appears significantly correlated, probably due to the same shut off protocol used (thermal inactivation).
Figure 1
Figure 1
Basic comparison statistics among the yeast S. cerevisiae mRNA half-lives during asynchronous growth measured by three independent laboratories.
Table 1
Table 1
Half-lives experimental measurements.
It has been shown that genes having the same biological function [11], [12] are likely to share similar half-life values. Consistently, by averaging using functional groups, we found an increase in correlation between Wang et al. and Grigull et al. datasets (An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e003.jpg, see Supplementary Figure S1A), and still no correlation between Munchel et al. and Wang et al. datasets (An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e004.jpg, see Supplementary Figure S1B).
Here, we developed a stochastic computational model of the expression kinetics to identify condition-specific mRNA stabilities which makes use only of experimental mRNA time profiles. We also assumed that degradation rates are gene-specific but approximately constant over the experiment time course. Predictions of our algorithm, termed DRAGON (Decay RAtes from Gene expressiON), were validated on experimental mRNA abundance [16] and turnover [17] data, both collected during the Intraerythrocytic Developmental Cycle (IDC) of Plasmodium falciparum. The estimations were in line with the experimental measurements. Remarkably, the DRAGON estimated half-lives were consistent with the finding of a peculiar pattern of mean half-life values along the wave of sequentially induced genes in subsequent stages of P. falciparum development. We also applied our methodology to public time-series datasets for which half-lives data, under the same experimental conditions, have not been experimentally measured. In particular we focused on budding yeast reproductive [18], [19] and metabolic cycle data [20]. In fact, for the yeast Saccharomyces cerevisiae, only half-life data under asynchronous growth are publicly available [2], [11], [12]. Our study showed the presence of the same periodic pattern of mean half-life values in all datasets, thus suggesting that such behavior may be a general feature, not limited to the Plasmodium falciparum IDC.
mRNA kinetics and half-life
Experimental evidence suggest that the majority of mRNAs are degraded with a first-order decay rate [2], [21]. This allows to characterize mRNA disappearance time profiles by a first-order rate equation
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e005.jpg
(1)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e006.jpg is the decay rate (or half-life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e007.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e008.jpg), An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e009.jpg is the mRNA concentration and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e010.jpg is the promoter activity (the rate of production of new mRNAs). It is worth noting that, the degradation rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e011.jpg cannot be estimated from the concentration time profile An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e012.jpg for a single gene, since the term An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e013.jpg is not usually available in the typical time-course microarray experiment. The measurement of the promoter activity time profile would require additional experiments (such as those described in [15]) but, in this paper, we will assume that only mRNA abundance time-series data are available. At steady-state An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e014.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e015.jpg are constant so that An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e016.jpg and, consequently, we obtain
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e017.jpg
(2)
From the above equation, it is clear that at steady-state an increase (decrease) in mRNA concentration can be produced either by an increase (decrease) of transcription or by a decreased (increased) value of the decay rate: the two regulatory strategies have therefore an equivalent outcome. As a result, from steady-states measurements, it is hopeless to reveal the relative contribution of transcription and degradation and, most importantly, their co-ordinated activity as well. By contrast, the whole kinetics of induction and relaxation, as measured by time-courses experiments, depends on the degradation and production rate in different ways: increasing (decreasing) the production rate results in a proportionally increased (decreased) mRNA abundance, whereas the rise time (i.e. the time required for the response to rise from An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e018.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e019.jpg of its final value) is not affected. Increasing the decay rate results in a faster rise time both in the induction and relaxation phases, whereas a decrease results in slower rise time [10]. This key point is illustrated in Figure 2A and in Supplementary Figures S2A and S2B.
Figure 2
Figure 2
Kinetics of gene induction.
Another important consequence of half-life specificity is the regulation of the timing of gene induction, as pointed out by Elkon et al. [22]. In fact, an expression wave, i.e. the sequential activation of genes, is usually interpreted as resulting from the corresponding activation of a multi-step transcription factors cascade (as illustrated in Figure 2B). Whereas such mechanism is certainly very important, there is also an alternative way to obtain an expression wave by means of a “stability gradient”. As illustrated in Figure 2C, a single transcription factor may initiate transcription of a set of target genes and their peak of induction can be modulated by a stability “gradient”, i.e. by specifically adjusted decay rates. More precisely, early induced genes would have short half-lives and late responding genes would have long half-lives. Clearly, both mechanism may well act in cells, thus generating a wide spectrum of responses.
Time-courses are a very common design for microarray analysis, which allows researchers to follow the dynamics of the cellular response to perturbations [22]. Such data are available for a very large number of experimental conditions and organisms: only the Stanford Microarray Database includes to date 1545 time course data sets. Among the examples later illustrated in the paper, it is worth mentioning the genome-wide gene expression time-series obtained during the reproductive cell cycle [18], [19], the metabolic cycle [20] and the P. falciparum IDC [16]. The time-series datasets used in this paper are summarized in Table 2.
Table 2
Table 2
Gene expression time-series experimental measurements.
DRAGON–an algorithm for half-life estimation
The goal of the DRAGON methodology is to derive a robust estimate of each mRNA species half-life starting from all available gene expression pairs. The rationale for the algorithm mainly draws on properties of pairs exhibiting a certain degree of common promoter activity (as in [23]). Besides, DRAGON infers common promoter activity using a statistical model that simulates both gene-specific and common effects.
The rate of change of mRNA concentration for a generic pair of genes, say gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e020.jpg and gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e021.jpg, is:
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e022.jpg
(3)
where the symbols An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e023.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e024.jpg represent the mRNA time profiles of the gene pair An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e025.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e026.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e027.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e028.jpg are the promoters activity, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e029.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e030.jpg are the degradation rate of mRNA of gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e031.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e032.jpg, respectively. The terms An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e033.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e034.jpg are not known since we considered the case in which only mRNA abundance is measured.
We modeled promoter activities as the sum of two terms, the first one common to the pair and the other one specific for each gene:
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e035.jpg
(4)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e036.jpg is the common part, scaled by constants An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e037.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e038.jpg, whereas An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e039.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e040.jpg are gene specific independent stochastic processes with zero mean, that is An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e041.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e042.jpg. Equations (3)(4) encompass the case of:
  • An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e043.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e044.jpg fully correlated (correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e045.jpg) for which An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e046.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e047.jpg
  • An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e048.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e049.jpg partially correlated (correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e050.jpg) for which An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e051.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e052.jpg
  • An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e053.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e054.jpg un-correlated (correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e055.jpg) for which either An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e056.jpg or An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e057.jpg.
Equations (3)(4) can therefore be written for all available gene pairs; thus, for a set of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e058.jpg genes, we have An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e059.jpg pairs to analyze. For each gene pair DRAGON provides an estimate of the time profile of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e060.jpg, of all the parameters, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e061.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e062.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e063.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e064.jpg, and the covariance matrix of the stochastic processes. For each gene we therefore have An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e065.jpg estimates of the decay rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e066.jpg, one for each pair containing that gene.
Notice that equations (3)(4) yield a couple of linear stochastic differential equations. Since measurements of mRNA concentrations are available only at given time points, it is necessary to transform (3)–(4) in a couple of discrete stochastic equations. The exact discretization of (3)–(4) is possible since they are linear [24]. The Kalman filter [25] is used on the resulting discrete equations and a maximum likelihood algorithm is exploited to generate the best possible estimate of the parameters.
A complete description of the mathematical model and of the discretization and parameters estimation procedure is given in the paragraph Stochastic modeling of expression kinetics and Kalman filtering of the Materials and Methods section.
Performance evaluation on malaria IDC experimental data
The IDC is characterized by four morphologic stages: ring, trophozoite, schizont and late schizont. The cycle begins with the red cells invasion by merozoites followed by a remodeling of the host cell in the ring stage [16]. The merozoites then develop into trophozoites. During the schizont stage, after a period of growth, the trophozoite undergoes an asexual dividing process and the parasite is ready for the next round of invasion by new merozoites (late schizont phase).
Bozdech et al. [16], using microarrays, measured genome-wide mRNA abundance profiles across 48 h during one cycle of P. falciparum IDC, collecting one sample per hour. Later on, Shock et al. [17] measured mRNA half-lives of 2774 transcripts of the IDC using chemical inhibitors to reach transcriptional shut-off.
The simultaneous availability of gene expression and decay data during the same biological process (IDC) represents a natural test bed for the validation of the DRAGON algorithm. Therefore, we applied DRAGON on Bozdech et al. dataset to obtain mRNA stability estimations (provided in Supplementary Table S6) to be compared with Shock et al. measurements for performance evaluation. The resulting Pearson correlation between in vitro and in silico measures is An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e067.jpg (P value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e068.jpg), and the first principal component accounts for An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e069.jpg of the variability, thus showing a good performance for DRAGON algorithm (see Figure 3). However, since gene expression and decay data have been measured by different groups, we can speculate that An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e070.jpg of unexplained variability may be partly due to inherent biological variability and to transcriptional inhibition stress. As further analysis, we computed average mRNA half-lives in both studies for functional categories (see Supplementary Figure S3). We found that the two studies are in better accordance when half-lives are averaged for all genes within any given functional category (Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e071.jpg).
Figure 3
Figure 3
DRAGON algorithm validation using P. falciparum IDC data.
Remarkably, Shock et al. in [17], found progressive stage-dependent average increases in mRNA stability and suggested such phenomenon to be a major determinant of mRNA accumulation (see Figure 4A). The same feature is also found using DRAGON estimated half-lives (see Figure 4B). To investigate in further detail the behavior of average half-life of genes sequentially induced during IDC, we computed for each gene the time point corresponding to its peak of expression (see the Data processing paragraph of the Materials and Methods section for details) and selected An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e072.jpg groups of genes having peak of expression at each hourly time points over the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e073.jpg hours monitored by Shock et al. For each gene group we computed half-lives mean and standard deviation and found a high correlation with the corresponding curve obtained using experimental data (Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e074.jpg, P value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e075.jpg; see Figure 4C). Early responding genes are characterized by high instability, whereas late responders are more stable, as also reported by Elkon et al. in [22] when studying mammalian cells. A possible explanation for the presence of stable mRNAs at the schizont stage, suggested by Shock et al., is that it may be important for the merozoite to receive a carefully regulated “starting package”, that would allow rapid activation of the IDC following the next round of invasion [17]. By contrast, the initial low mRNA stability values may be an indication of the fast dynamic remodeling after merozoite invasion [17]. To evaluate the probability of obtaining such behavior by chance, we randomized the gene expression matrix and used DRAGON to estimate half-lives (see Figure 4D). Consistently, the estimation of half-lives using random data does not produce any correlation with experimental data (Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e076.jpg).
Figure 4
Figure 4
Periodic behavior of average half-lives of sequentially induced genes in P. falciparum IDC.
Half-lives estimation during reproductive cycle in S. cerevisiae
Gene expression during yeast cell cycle has been recently measured by Pramila et al. [18] using alpha-factor synchronization and by Orlando et al. [19] using centrifugal elutriation for synchronization. We obtained a high consistency of DRAGON estimations using data for 569 transcripts over replicate datasets (Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e078.jpg for Pramila et al. dataset and Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e079.jpg for the Orlando et al. dataset; see Figure 5A–B). The larger variability in half-lives estimations may be explained by the inconsistencies between replicate time-series in the Pramila et al. dataset with respect to the Orlando et al. dataset (see Figure 5C). All half-lives estimations obtained with the DRAGON algorithm are provided in Supplementary Tables S1 and S2 (Pramila) and in Supplementary Tables S3 and S4 (Orlando). Notwithstanding significant differences in synchronization procedures, we also found a high correlation of DRAGON half-lives estimations over the two datasets (Pearson correlation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e080.jpg, P value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e081.jpg; see Figure 5D) where the first principal component accounts for An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e082.jpg of the overall variability. We can speculate that An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e083.jpg of unexplained variability may be partly due to the different synchronization methods used. In fact, Orlando et al. obtained a cell cycle duration of about 2 hours, 8 samples per cycle [19], whereas Pramila et al. obtained a cell cycle duration of about 1 hour, 12 samples per cycle. Consistently, most of the transcripts during the slower cycle display higher half-lives when compared to the fastest cycle (see Figure 5D).
Figure 5
Figure 5
DRAGON performance over 569 cell cycle regulated genes from the Pramila and Orlando datasets.
GO annotations of genes with extreme half-lives in S. cerevisiae
In this paragraph we briefly discuss functional annotations (done using GOrilla software [26]) of novel predicted half-lives provided by DRAGON algorithm using yeast reproductive and metabolic cycle time series. For the yeast cell cycle we normalized the half-life log-distribution (Z-score), for each dataset, and then computed the geometric mean to obtain a single half-life value for each gene. Notably, the averaging has the effect of reducing the impact of the different synchronization stress response. The list of half-lives normalized values (geometric mean value equal to 1) for common genes to all datasets is provided as Supplementary Table S7 in the Half-life estimation paragraph of the Materials and Methods section.
Unstable genes are enriched with replication fork complex (p-value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e084.jpg) and stable genes (histones HA1-2,HB1-2) are enriched with nucleosome (p-value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e085.jpg). This is consistent with the need of producing a large number of histones during DNA replication process so that stable histone mRNAs contribute to a higher translation efficiency. Moreover, DNA replication timing requires first the formation of the replication fork, then the production of the needed histones for chromatin assembling: such temporal sequence of events is consistent with a rapid turnover of the replication complex genes and a slow turnover of the histone genes (see Supplementary Figure S4). Among unstable genes we also found the G1/S transition cyclins and among stable ones we found G2/M transition cyclins (see Supplementary Figure S5). In this case, the temporal sequence of events is the progression of the cell cycle from DNA replication to mitosis.
For the yeast metabolic cycle (half-lives estimations using DRAGON algorithm are provided in Supplementary Table S5) we found many stable mRNA species involved in the organic acid and arginine metabolism and protein catabolic processes. Among unstable messengers we found genes involved in DNA repair (p-value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e086.jpg), DNA metabolism (p-value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e087.jpg) and chromatin silencing (p-value An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e088.jpg).
Periodic behavior of average half-lives of sequentially induced genes
The increasing pattern of average half-life found during P. falciparum IDC (shown in Figure 4A) motivated us to investigate whether a periodicity could be found also in other cyclical biological processes. We focused on the reproductive cell cycle and the metabolic cycle in Saccharomyces cerevisiae, for which high resolution time series measurements are available on public repositories (see Table 2).
To study if a periodic pattern of average half-life of sequentially induced genes exists along the cell cycle progression, for each gene we computed the time points at which maximal expression is attained (see the Data processing paragraph of the Materials and Methods section for details). Thus, we obtained, for each time point, the list of genes having expression peak value at that time and computed the corresponding mean and variance of DRAGON estimated half-lives values. Indeed, we found a cyclic behavior along sequentially induced genes in both datasets (see Figure 6A for the Pramila et al. dataset and Figure 6B for the Orlando et al. dataset). Synchronization methods, cell cycle duration and number of samples are different between the two cited studies, but, reassuringly, the phases of the cell cycle at which mean half-life is minimal or maximal is consistent. In fact, for both datasets we observed a cyclical increase of mean half-life from G1 phase to M phase and a subsequent decrease back to G1. The figure clarifies that the minimal mean half-life is reached at the G1/S transition, whereas the maximal value correspond to the M/G2 phase for both cycles and datasets. The latter is consistent with the observation that, in higher eukaryotes, mitosis is accompanied by global repression of nuclear RNA synthesis [27], indicating that mRNAs must be stable to be inherited from daughter cells.
Figure 6
Figure 6
Periodic behavior of average half-life during the reproductive and metabolic cycle.
The yeast metabolic cycle has been recently studied by Tu et al. [20] using a continuous culture system, after a brief starvation period, the culture spontaneously began persistent respiratory cycles of about 5 hours. In the same study, a genome-wide microarray gene expression measurement was performed. Samples were taken every 25 minutes for 3 consecutive cycles. Using DRAGON algorithm we estimated half-lives using data of 1043 transcripts. Surprisingly, also in this case we found a cyclical pattern for mean half-life of sequentially induced genes. The maximum peak is located at the RC phase and the minimum peak located at RB phase (see Figure 6C).
Integrated analysis–sequential waves of co-ordinated transcription and decay
Recently, the appearance of a number of studies has revealed the fundamental role of stability regulation in shaping appropriate cell response [1]. A key point has been recently addressed by Shalem et al. [10], who have shown the dynamic co-ordinated interplay between transcription and degradation. They have found in yeast two basic regulatory strategies in response to stress. More precisely, they measured changes of mRNA abundance and decay rates in a yeast population subjected to oxidative and DNA damage stress. By grouping genes according to the time point at which the maximal (minimal) fold change is attained and combining normalized (mean and variance) mRNA abundance and decay rate data, they constructed a “stability versus folding” (SF) diagram where change in mRNA stability relative to a reference state (mean value in our case) is plotted against the maximal fold change. Using yeast expression time-course data obtained in response to an oxidative stress and a DNA damage, they were able to reveal two different strategies: a) a “counteracting regulation” strategy (see Figure 7A), characterized by genes in which an increase (decrease) in degradation rates counteracts a increase (decrease) in mRNA abundance, i.e. repressed genes are stabilized and induced genes are destabilized; b) a “synergistic regulation” strategy (see Figure 7B), characterized by genes in which an increase (decrease) in degradation rates is associated with an decrease (increase) in mRNA abundance, i.e. induced genes are stabilized and repressed genes are destabilized.
Figure 7
Figure 7
Illustration of the counteracting and synergistic regulatory strategy.
Shalem et al. also found that, progressing from early time points forward, the negative correlation (counteracting) was replaced with a positive correlation (synergistic). Such co-ordination strategy may permit crosstalk between different steps of mRNA biogenesis, providing a mechanism to control the order and timing of events [28]. The work of Shalem et al. has shown the importance of combining expression data with decay rates under the same experimental condition to reveal the underlying strategy of co-ordination of the two “regulatory arms”, namely transcription and degradation. Uncovering such relationships is certainly a fundamental task, since the underlying reciprocal influences between mRNA production and degradation are largely unexplored [10]. The DRAGON algorithm, by estimating half-lives directly from gene expression data under specific conditions, allows the computational integration of mRNA abundance and decay rates data, making this powerful combined analysis possible when experimentally measured half-lives are not available.
We computed SF diagrams for P. falciparum IDC, yeast cell cycle (Pramila et al. dataset) and metabolic cycle (shown in Figure 8). In panels A,C and E each blue dot corresponds to a Pearson correlation of the SF diagram at the peak time point indicated on the x-axis, for the three datasets. In panels B,D and F the SF diagrams corresponding to the correlation values indicated by the arrows in panels A,C and E are displayed. The arrows point to maximal negative (red dots in panels B,D,F and red arrows in panels A,C,E) and maximal positive correlation values (green dots in panels B,D,F and green arrows in panels A,C,E).
Figure 8
Figure 8
Temporal progress on the stability/folding (SF) diagrams.
Strikingly, in all cases we reached the same conclusions of Shalem et al., namely we found that early induced genes show counteracting regulation, whereas late induced genes show a synergistic regulation.
Advantages and disadvantages of the method
The main advantage of the DRAGON algorithm consists in the estimation of the mRNA half-lives directly from gene expression time-course during condition-specific experiments. Moreover it estimates the correlation among promoter activities between pairs of genes. Another advantage of the algorithm lies in its robustness. Specifically, we observed that even if the accuracy of the absolute values of the estimated half-lives can be influenced by many factors (such as the number of points in the time series, the accuracy of the measurements, the time interval between samples, the choice of the thresholds for the outliers, etc.), the ranking of half-lives is insensitive to the factors mentioned above.
The main disadvantages are the following: DRAGON can work only with time-series under the same experimental condition and cannot handle steady-state values under different conditions. As a general rule a reliable estimate requires at least 10–12 time samples, i.e. a number significantly larger than the number of parameters to be estimated (this rule is not obviously always applicable as the required number of points depends strongly on the signal to noise ratio) and a sampling time not larger than the expected average half-life. If no information is available about the correlation of promoter activities, as a rule of thumb, a set of at least 50–100 time series must be processed together in order to have reliable half-lives estimates. One basic hypothesis is that the half-life of a transcript is approximately constant during the time course of the experiment, thus a substantial change of its value would yield an unreliable estimate. These problems can be handled by performing more measurements using a shorter sample time, or by considering moving time windows. The computational overhead can be significant: for a sample of 1000 time series there are An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e089.jpg pairs to analyze, requiring a computation of about 150 hours on a medium-speed single-processor machine capable of analyzing 2 pairs per second.
Conclusion
Our analysis supports and strengthens Shalem et al. conclusions about the coordination of transcriptional and mRNA degradation in the cell in response to stress. We have demonstrated that during periodic processes, such as the P. falciparum IDC, the reproductive cell cycle and the metabolic cycle, the alternative interplays between changes in mRNA stability and changes in mRNA abundance are activated by periodically switching from a counteracting to a synergistic regulation. In light of these results, the classical vision of periodic processes as the result of serial transcription factor sequential activation, should be re-considered from a broader point of view by including post-transcriptional regulation and coordination.
Stochastic modeling of expression kinetics and Kalman filtering
We defined as An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e090.jpg the time profile of the expression of gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e091.jpg at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e092.jpg. The underlying conservation equation simply stems from the observation that the rate of change of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e093.jpg with time, i.e. its time derivative An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e094.jpg, must equal the difference between the production and degradation term. Based on experimental evidence [2], the degradation is well described by a first order term. The dynamics of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e095.jpg-th transcript is therefore described by
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e096.jpg
(5)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e097.jpg is the mRNA decay rate of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e098.jpg-th messenger. This value is linked to the half-life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e099.jpg of the transcript by the relation An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e100.jpg. An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e101.jpg is the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e102.jpg-th gene promoter activity regulated by transcription factors. Such regulation occurs by triggering or suppressing the transcription of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e103.jpg-th gene, thus we have An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e104.jpg. Moreover, the observed measure An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e105.jpg is also a noisy time-series, thus we have
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e106.jpg
(6)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e107.jpg is the standard deviation of measurements white noise An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e108.jpg (see supplementary material Text S1 for an example of the identification procedure). We considered a generic pair of expression time profiles characterized by the presence of two terms: a stochastically correlated promoter activity An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e109.jpg and a gene-specific term An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e110.jpg. We then considered the case:
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e111.jpg
(7)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e112.jpg is a scaling factor accounting for the relative contribution to the overall promoter activity regarding gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e113.jpg. The term An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e114.jpg models the part of the promoter activity which is not common to the pair. We model this part by means of a noise term, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e115.jpg which is assumed to be a white noise. The common part An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e116.jpg is modeled as a Wiener process:
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e117.jpg
(8)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e118.jpg is white noise. Thus An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e119.jpg. The complete mathematical dynamic model for two transcripts An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e120.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e121.jpg, together with their respective measurement equations, is
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e122.jpg
(9)
We can rewrite the linear dynamic system (9) using a compact matrix notation
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e123.jpg
(10)
where
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e124.jpg
and
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e125.jpg
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e126.jpg
Since the dynamic system (10) is linear, it can be exactly discretized (see [24]) for a given time interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e127.jpg, corresponding to the time interval between two consecutive measurements. The An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e128.jpg-th measurements corresponds to An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e129.jpg, thus in the discretized system we can use An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e130.jpg in place of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e131.jpg, to keep the notation simple.
The solution of the linear dynamic system (10) is
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e132.jpg
(11)
and its discretized form is
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e133.jpg
(12)
where
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e134.jpg
and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e135.jpg is the covariance matrix defined by
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e136.jpg
(13)
The unknown parameters of the model to be estimated are An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e137.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e138.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e139.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e140.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e141.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e142.jpg. The state variables of the system are An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e143.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e144.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e145.jpg. For each given choice of the parameters we used the Kalman filter [25] to estimate of the state variables.
The Kalman filter equation uses a feedback control strategy. It contains a prediction term for projecting forward (in time) the current state to obtain the a priori estimate, and a correcting term for incorporating a new measurement into the a priori estimate to obtain an improved a posteriori estimate
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e146.jpg
(14)
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e147.jpg is the prediction Kalman gain that depends on the parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e148.jpg of the stochastic equation.
For each choice of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e149.jpg we run the Kalman filter. A probability value is associated to the resulting estimation. These values measures the probability that the current parametrization of the model generates the measured time series. Denoting by An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e150.jpg the innovation of the stochastic process, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e151.jpg is a sequence of independent gaussian random variables with covariance An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e152.jpg. The optimal set An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e153.jpg of parameters if therefore chosen according to a maximum likelihood criterion as the choice corresponding to the maximum of the a priori probability density of the innovation sequence. This corresponds to the minimum of the likelihood function
A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e154.jpg
where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e155.jpg is the number of samples. We are interested in the half life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e156.jpg of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e157.jpg-th messenger. To use all the available information and make the method robust with respect to measurement and estimation errors, we have designed the following algorithm (see Supplementary Figure S6):
  • Given a set of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e158.jpg mRNA time profiles, perform the maximum likelihood estimation for every pair An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e159.jpg and compute the corresponding An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e160.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e161.jpg.
  • For each pair An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e162.jpg compute the ratio matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e163.jpg whose elements are the ratios between the half-lives of gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e164.jpg and gene An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e165.jpg. The matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e166.jpg is generally not symmetric due to the presence of outliers and numerical sensitivity. Thus we defined the ratio estimation by row as An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e167.jpg and by column as An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e168.jpg. The matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e169.jpg contains all the ratios An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e170.jpg on the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e171.jpg-th row, and all the ratios An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e172.jpg on the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e173.jpg-th column. Let us denote An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e174.jpg the sum of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e175.jpgth row, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e176.jpg the sum of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e177.jpg column, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e178.jpg the mean operator, that is,
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e179.jpg
    (15)
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e180.jpg
    (16)
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e181.jpg
    (17)
  • Given An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e182.jpg, delete outliers to obtain a final matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e183.jpg. First, compute the probability density (using a smoothing kernel approach) of all the entries An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e184.jpg and delete those values below a probability of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e185.jpg of occurring in the distribution. Second, since ideally An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e186.jpg, we considered as outliers those pairs such that An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e187.jpg.
  • On the resulting An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e188.jpg matrix compute for each transcript An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e189.jpg two estimates of its half-life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e190.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e191.jpg, using equations (15), (16)and (17). We obtained
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e192.jpg
    (18)
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e193.jpg
    (19)
    This computation requires the value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e194.jpg. When this value is known for the group of transcripts under analysis the measured value can be used. Otherwise, letting An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e195.jpg one can obtain half-life values that are relative to the average half-life of the group. However, we have followed a third approach. All the results reported in this paper have been obtained by replacing An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e196.jpg with the geometric mean of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e197.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e198.jpg, that is
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e199.jpg
    (20)
    The final estimate of the half-life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e200.jpg for the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e201.jpg-th gene is computed as the weighted average of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e202.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e203.jpg using as weights the respective variances An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e204.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e205.jpg as follows
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e206.jpg
    (21)
    where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e207.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e208.jpg are the standard deviation of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e209.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e210.jpg, respectively.
  • We considered as a quality index for each estimated half-life An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e211.jpg the following:
    A mathematical equation, expression, or formula.
 Object name is pcbi.1002772.e212.jpg
    where An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e213.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e214.jpg are noise variances of the discrete system (12) and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e215.jpg (see equation (13)) is the mutual covariance of the state noise between time series An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e216.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e217.jpg. Thus, high values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e218.jpg imply the presence of a correlation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e219.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e220.jpg in equation (9). We removed the half-lives having a An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e221.jpg value smaller than the An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e222.jpg percentile of its distribution.
Datasets
Public experimental data used throughout the paper are described in Table 1 (experimental half-lives measurements) and in Table 2 (gene expression time series). Pramila et al. in [18] and Orlando et al. in [19] experimentally measured genome-wide gene expression data during the reproductive cell cycle. We considered the ranking provided by the combined test developed by de Lichtenberg et al. [29] for each replicate for the two datasets and, among the list of 1000 genes with highest ranking, we selected those common to all datasets. We ended up with a list of 569 genes that we used for half-life estimation. Tu et al. in [20] experimentally measured genome-wide gene expression data during the metabolic cell cycle. We selected 1000 genes with the best periodicity score according to [20]. Of the 1000 genes, DRAGON estimated half-lives are 939.
Data processing
Half-lives determination of genes induced during each stage of P. falciparum IDC
Shock et al. in [17] experimentally measured genome-wide values of decay rates for each gene in each of the four stages of the IDC. To obtain a single half-life value for each messenger, we performed a k-means clustering of microarray gene expression data [16] by considering 5 stages (according to [16]: early ring, ring, trophozoite, schizont and late schizont). Then, we merged the early ring cluster with the ring cluster to obtain the same stages as in Shock et al.. Among the 4488 genes in [16] we chose 1000 genes with the best periodicity score (power signal/power total ratio) according to [16]. Of the 1000 genes, DRAGON estimated half-lives are 967, available experimental half-lives are 675. Both data are available over a set of 616 genes.
Expression peak timing estimation
To estimate peak timing, for a given noisy gene expression time profile, we preliminary performed the smoothing algorithm presented by Bar-Joseph et al. in [30]. The algorithm employs two parameters: grid An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e223.jpg (number of spline curves) and classes An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e224.jpg (number of classes to use for clustering). In particular, for Pramila datasets we used An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e225.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e226.jpg, for Orlando datasets we used An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e227.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e228.jpg, for Tu dataset we used An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e229.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e230.jpg, for Malaria dataset we used An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e231.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1002772.e232.jpg.
Half-life estimations
DRAGON estimated half-lives are provided as supplementary materials Tables S1, S2, S3, S4, S5, S6, S7, described Supporting Information section.
Matlab code will be provided upon request.
Additional data and information can be found at web site http://www.dis.uniroma1.it/~farina/dragon.
Figure S1
Functional categories analysis in yeast S. cerevisiae during asynchronous growth measured by three laboratories. Three genome-wide studies are considered: Grigull et al., Wang et al. and Munchel et al. (A) Average mRNA half-lives in both studies Wang et al. and the Grigull et al. datasets for 111 functional categories from the yeast GO Biological Process database (http://www.geneontology.org) that are represented in the set of 2863 transcripts by 5 or more members. (B) compare, in the same way, the Munchel et al. and the Wang et al. datasets.
(TIF)
Figure S2
Kinetics of gene induction. Panels A–B show in silico experiments to illustrate some basic features of gene induction kinetics. The reference time profile with unity steady-state is plotted in black. The “ON” and “OFF” regions correspond to the turning “ON” or “OFF” of the promoter activity. (A) Induction kinetic of transcripts having the same half-life value and, as a consequence, the same speed of response. The higher (or lower) steady-state value of the red and blue time profiles is due only to an increased (or decreased) transcription rate. (B) Induction kinetic of transcripts having different half-lives. The time profile plotted in red corresponds to an unstable transcript. It has a faster induction and relaxation profile but a lower steady-state value. By contrast, the blue one has an higher half-life value, resulting in a higher steady state value but a slower response. The example illustrates that, to obtain both a fast response and an high steady-state value, the regulatory strategy must destabilize transcriptionally up-regulated genes.
(TIF)
Figure S3
Functional categories analysis for DRAGON estimations using P. falciparum IDC data. Average mRNA half-lives in both studies, DRAGON iestimations versus and experimentally measured by Shock et al. half-lives, for 12 functional categories from the P. falciparum GO annotation database (http://www.geneontology.org) that are represented in the set of 616 transcripts by 5 or more members.
(TIF)
Figure S4
GO annotations of genes with extreme half-lives in S. cerevisiae DNA replication timing requires first the formation of the replication fork, then the production of the needed histones for chromatin assembling: such temporal sequence of events is consistent with a rapid turnover of the replication complex genes and a slow turnover of the histone genes.
(TIF)
Figure S5
GO annotations of genes with extreme half-lives in S. cerevisiae Among unstable genes we also found the G1/S transition cyclins and among stable ones we found G2/M transition cyclins. In this case, the temporal sequence of events is the progression of the cell cycle from DNA replication to mitosis.
(TIF)
Figure S6
DRAGON algorithm pipeline.
(TIF)
Table S1
DRAGON estimated half-lives using alpha30 dataset.
(XLSX)
Table S2
DRAGON estimated half-lives using alpha38 dataset.
(XLSX)
Table S3
DRAGON estimated half-lives using Orlando replicate 1 dataset.
(XLSX)
Table S4
DRAGON estimated half-lives using Orlando replicate 2 dataset.
(XLSX)
Table S5
DRAGON estimated half-lives using metabolic cycle dataset.
(XLSX)
Table S6
DRAGON estimated half-lives using P. falciparum dataset.
(XLSX)
Table S7
DRAGON estimated normalized half-lives using all yeast datasets.
(XLSX)
Text S1
Example of parameter estimation through the Kalman filter.
(PDF)
Acknowledgments
We would like to thank Prof. Pino Macino and Dr. Teresa Colombo for helpful suggestions on the manuscript and the “Consorzio interuniversitario per le Applicazioni di Supercalcolo Per Università e Ricerca” (CASPUR) for providing computing resources and Dr. Alessandro Federico for support.
The authors are partially supported by The Epigenomics Flagship Project (Progetto Bandiera Epigenomica) EPIGEN funded by Italian Ministry of Education, University and Research (MIUR) and the National Research Council of Italy (CNR).
Funding Statement
The authors are partially supported by The Epigenomics Flagship Project (Progetto Bandiera Epigenomica) EPIGEN funded by Italian Ministry of Education, University and Research (MIUR) and the National Research Council of Italy (CNR). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
1. Garneau N, Wilusz J, Wilusz C (2007) The highways and byways of mrna decay. Nat Rev Mol Cell Bio 8: 113–126. [PubMed]
2. Munchel S, Shultzaberger R, Takizawa N, Weis K (2011) Dynamic profiling of mrna turnover reveals gene-specific and system-wide regulation of mrna decay. Mol Biol Cell 22: 2787–2795. [PMC free article] [PubMed]
3. Keene J (2010) The global dynamics of rna stability orchestrates responses to cellular activation. BMC Biol 8: 95. [PMC free article] [PubMed]
4. Gerber A, Herschlag D, Brown P (2004) Extensive association of functionally and cytotopically related mrnas with puf family rna-binding proteins in yeast. PLoS Biol 2: 342–354. [PMC free article] [PubMed]
5. Hogan D, Riordan D, Gerber A, Herschlag D, Brown P (2008) Diverse rna-binding proteins interact with functionally related sets of rnas, suggesting extensive regulatory system. PLoS Biol 6: e255. [PMC free article] [PubMed]
6. Houseley J, Tollervay D (2009) The many pathways of rna degradation. Cell 136: 763–776. [PubMed]
7. Cheneval D, Kastelic T, Fuerst P, Parker C (2010) A review of methods to monitor the modulation of mrna stability: a novel approach to drug discovery and therapeutic intervention. J Biomol Screen 15: 609–622. [PubMed]
8. Narsai R, Howell K, Millar A, O'Toole N, Small I, et al. (2007) Genome-wide analysis of mrna decay rates and their determinants in arabidopsis thaliana. Plant Cell 19: 3418–3436. [PubMed]
9. Sharova L, Sharov A, Nedorezov T, Piao Y, Shaik N, et al. (2009) Database for mrna half-life of 19977 genes obtained by dna microarray analysis of pluripotent and differentiating mouse embryonic cells. DNA Res 16: 45–58. [PMC free article] [PubMed]
10. Shalem O, Dahan O, Martinez M, Furman I, Segal E, et al. (2008) Transient transcriptional re- sponses to stress are generated by opposing effects of mrna production and degradation. Mol Syst Biol 4: 223. [PMC free article] [PubMed]
11. Wang Y, Liu C, Storey J, Tibshirani R, Hershlag D, et al. (2002) Precision and functional specificity in mrna decay. Proc Natl Acad Sci U S A 99: 5860–5865. [PubMed]
12. Grigull J, Mnaimneh S, Pootoolal J, Robinson M, Hughes T (2004) Genome-wide analysis of mrna stability using transcription inhibitors and microarrays reveals posttranscriptional control of ribosome biogenesis factors. Mol Cell Biol 25: 5534–5547. [PMC free article] [PubMed]
13. Coller J (2008) Methods to determine mrna half-life in saccharomyces cerevisiae. Methods Enzymol 448: 267–284. [PubMed]
14. Rabani M, Levin J, Fan L, Adiconis X, Raychowdhury R, et al. (2011) Metabolic labeling of rna uncovers principles of rna production and degradation dynamics in mammalian cells. Nature Biotechnol 29: 436–442. [PMC free article] [PubMed]
15. Garcia-Martinez J, Aranda A, Perez-Ortin J (2004) Genomic run-on evaluates transcription rates for all yeasts genes and identifies gene regulatory mechanisms. Mol Cell 15: 303–313. [PubMed]
16. Bozdech Z, Llinas M, Pulliam B, Wong E, Zhu J, et al. (2003) The transcriptome of the intraery- throcytic developmental cycle of Plasmodium falciparum. PLoS Biol I: E5. [PMC free article] [PubMed]
17. Shock J, Fischer K, DeRisi J (2007) Whole-genome analysis of mrna decay in Plasmodium falci- parum reveals a global lengthening of mrna half-life during the intra-erythrocytic developmental cycle. Genome Biol 8: R134. [PMC free article] [PubMed]
18. Pramila T, Wu W, Miles S, Noble W, Breeden L (2006) The forkhead transcription factor hcm1 regulates chromosome segregation genes and fills the s-phase gap in the transcriptional circuitry of the cell cycle. Genes Dev 20: 2266–2278. [PubMed]
19. Orlando D, Lin C, Bernard A, Wang J, Socolar J, et al. (2008) Global control of cell-cycle tran- scription by coupled cdk and network oscillators. Nature 453: 944–947. [PMC free article] [PubMed]
20. Tu B, Kudlicki A, Rowicka M, McNight S (2005) Logic of the yeast metabolic cycle: temporal compartmentalization of cellular processes. Science 310: 1152–1158. [PubMed]
21. Ross J (1995) mrna stability in mammalian cells. Microbiol Rev 59: 423–450. [PMC free article] [PubMed]
22. Elkon R, Zlotorynski E, Zeller K, Agami R (2010) Major role for mrna stability in shaping the kinetics of gene induction. BMC Genomics 11: 259. [PMC free article] [PubMed]
23. Farina L, Santis AD, Salvucci S, Morelli G, Ruberti I (2008) Embedding mrna stability in corre-lation analysis of time-series gene expression data. PLoS Comput Biol 4: e1000141. [PMC free article] [PubMed]
24. Kailath T (1980) Linear systems. Prentice-Hall. 682 p.
25. Kalman R (1960) A new approach to linear filtering and prediction problems. Transactions of the ASME–Journal of Basic Engineering 82: 35–45.
26. Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z (2009) Gorilla: A tool for discovery and visualization of enriched go terms in ranked gene lists. BMC Bioinformatics 2009: 10–48. [PMC free article] [PubMed]
27. Shermoen A, O'Farrell P (1991) Progression of the cell cycle through mitosis leads to abortion of nascent transcripts. Cell 67: 303–310. [PMC free article] [PubMed]
28. Dahan O, Gingold H, Pilpel Y (2011) Regulatory mechanisms and networks couple the different phases of gene expression. Trends Genet 27: 316–322. [PubMed]
29. de Lichtenberg U, Jensen L, Fausboll A, Jensen T, Bork P, et al. (2004) Comparison of computa-tional methods for the identification of cell cycle-regulated genes. Bioinformatics 21: 1164–1171. [PubMed]
30. Bar-Joseph Z, Gerber G, Gifford D, Jaakkola T, Simon I (2003) Continuous representations of time-series gene expression data. J Comput Biol 10: 341–356. [PubMed]
Articles from PLoS Computational Biology are provided here courtesy of
Public Library of Science