<title>Introduction</title> <p>Gene expression is a highly regulated process composed of two fundamental biological events: synthesis and degradation. Transcription regulation is achieved by modulating the frequency of transcription initiation and, although the most studied, this event represents just the first of the many complex stages leading to a mature mRNA. Recent experimental work is beginning to shed light on the complex architecture underlying mRNA degradation pathways by identifying the factors and enzymes involved. Therefore, it is now widely accepted that mRNA decay contribution to the control of gene expression is not simply a biological waste-disposal system, but a key player for the temporal coordination of cellular functions. Moreover, a number of highly complex and sophisticated specific mechanisms have been identified ##REF##17245413##[1]##. Such mechanisms include the interaction with mRNA binding proteins ##REF##15024427##[2]## and the nonsense-mediated mRNA decay pathway ##UREF##0##[3]##, both able to affect the accumulation of hundreds of transcripts.</p> <p>Recent technologies, such as microarrays, are able to provide measurements of mRNA abundance over time under different experimental conditions. In order to decipher the intricate regulatory network underlying the highly coordinate cell behavior, effective computational methods have been developed to take advantage of gene expression data. The basic idea underlying such methods stems from the experimental observation that genes are organized in groups showing similar time profiles ##REF##9843569##[4]## (called “clusters”). These groups often share some common biological features, such as the same cellular function or the presence of a common motif at their promoter regions ##UREF##1##[5]## where transcription factors (TFs) can bind and possibly turn them on or off in a coordinated manner, when needed. For this reason, it is now widely accepted that <italic>co-expression</italic> is a good indication for <italic>co-regulation</italic>
##REF##9847082##[6]##–##REF##11340206##[8]##, meaning that whenever two genes display similar time profiles it is likely that they are both targets of the same transcription factor(s). The search for co-regulated genes depends on association metrics used by clustering algorithms ##UREF##1##[5]##,##REF##9843981##[9]##,##REF##10077610##[10]## and gene network inference algorithms ##UREF##2##[11]##–##UREF##3##[13]##. Therefore, measuring the degree of co-expression of genes is a fundamental step for data analysis, and in fact, many similarity measures have been proposed in the literature ##UREF##4##[14]##. Among those available to quantitatively measure simultaneous expression, we will refer to the usual R<sup>2</sup> value obtained from a linear regression model between two given gene expression time profiles denoted by <italic>m<sub>A</sub></italic>(<italic>t</italic>) and <italic>m<sub>B</sub></italic>(<italic>t</italic>). Their co-varying degree is therefore measured as the fraction of the total variance explained by the regression <italic>m<sub>A</sub></italic>(<italic>t</italic>) = <italic>c</italic>
<sub>1</sub>
<italic>m<sub>B</sub></italic>(<italic>t</italic>)+<italic>c</italic>
<sub>2</sub>. Such coefficient, indicated in this paper as the <italic>simultaneous</italic> R<sup>2</sup> of the corresponding gene pair, is the square of the Pearson correlation and takes values between 0 and 1.</p> <p>In order to infer the gene regulatory network, several laboratories have combined microarray data with protein-DNA interaction data, taking advantage of ChIP-on-chip experiments ##REF##16522208##[15]##. Such studies have shown that the same transcription factor (or combinations of) may target genes with very different expression time profiles, even in the same experimental condition. For example, the targets of the yeast cell cycle transcriptional regulators MBF/SBF display expression peak times that span from early G1 to late S. Moreover, delays have been recently observed between putatively co-regulated genes ##REF##11743722##[16]##,##REF##12054769##[17]##. One fundamental biological mechanism underlying such temporal spread is certainly combinatorial regulation of transcription factors. In fact, various TFs can modulate target response by cooperating or competing for DNA binding. Consequently, new computational techniques have recently appeared in the literature to tackle this problem ##UREF##5##[18]##–##REF##17551510##[25]##. However, combinatorial regulation is not the only mechanism responsible for peak time delay, as other regulation layers are active throughout transcript life and impact its abundance over time. One such additional regulation layers is certainly the post-transcriptional one, that is the stability properties of transcripts that may specifically contribute to the determination of their timing and amount during cell response to various internal and/or external signals. Strikingly, recent genome-wide measurement of the yeast transcripts half-lives ##REF##11972065##[26]##,##REF##15169913##[27]## has shown functional specificity in mRNA decay. Together, these results pointed to a general relationship between physiological function and mRNA decay rate thus providing strong evidence that precise control of mRNA turnover is a fundamental feature of gene expression programs in yeast ##REF##11972065##[26]## and in many other organisms.</p> <p>Here we focus on the development of a novel computational tool aiming to uncover co-regulated genes through transcriptional and post-transcriptional regulatory mechanisms. To this purpose, starting from the computational approach developed by Farina <italic>et al.</italic>
##REF##17370425##[28]##, we introduce a new relationship between gene pairs, called <italic>lead-lag relationship</italic>. The term “lead-lag” has been taken from the field of control systems engineering where the same relationship holds between the input and the output of the so called “lead-lag compensator”, which is the fundamental building block for the design of automatic control systems ##UREF##7##[29]##. In a biological perspective, the lead-lag relationship should be referred to genes under a common regulatory signal (“input”) involved in the same biological function (“output”) as, for example, in the dynamic multi sub-units complex formation ##UREF##8##[30]##,##REF##11779829##[31]##. Using yeast cell-cycle time-series gene expression data, we demonstrate that this new similarity metric is able to capture the dynamics of gene expression, including varying mRNA stability and delays. Thus, the predictive power of lead-lag R<sup>2</sup> for the identification of co-regulated genes is significantly higher than that of standard similarity measures, allowing the selection of a large number of entirely new putatively co-regulated genes. Furthermore, the lead-lag metric can also be used to uncover the relationship between gene espression time-series and the formation of protein complexes.</p>
<title>Materials and Methods</title> <title>Computation of the Simultaneous and Lead-Lag R<sup>2</sup> Between Gene Expression Time Profiles</title> <p>The mRNA relative abundance time course data obtained from cell populations experiments for gene A and B is denoted by <italic>m<sub>A</sub></italic> and <italic>m<sub>B</sub></italic>, respectively. The <italic>simultaneous R</italic>
<sup>2</sup>, is the usual squared Pearson correlation coefficient which measures the fraction of the total variance explained by a linear fit between the two variables <italic>m<sub>A</sub></italic> and <italic>m<sub>B</sub></italic>, that iswhere <italic>η</italic> accounts for intrinsic and extrinsic noise.</p> <p>The rationale behind the <italic>lead-lag R</italic>
<sup>2</sup> is the following. We considered two genes, <italic>A</italic> and <italic>B</italic>, subject to the same regulatory signal (promoter activity) – possibly of different strength – due to the presence at their promoters of the same TF complex in its active state. Moreoever, we assumed that the change in mRNA levels due to the degradation rate could be reasonably well captured by a first order rate kinetics ##REF##7565413##[53]##, and consequently the dynamic equation that includes both synthesis and degradation is the followingwhere the two variables <italic>m<sub>A</sub></italic> and <italic>m<sub>B</sub></italic> measure gene expression on a linear scale (fold induction), <italic>P<sub>X</sub></italic> is the promoter activity time profile of the TF complex relative to gene <italic>X</italic>, <italic>α<sub>X</sub></italic> is its maximal strength, <italic>k<sub>X</sub></italic> is the degradation rate (<italic>k<sub>X</sub></italic> = log(2)/<italic>t<sub>1/2</sub></italic>) and <italic>η<sub>X</sub></italic> accounts for intrinsic and extrinsic noise. In order to remove size effects, the common signal between the promoter activities of the two genes is indicated as <italic>p</italic>(<italic>t</italic>) and is such thatso that we getFrom the second equation of (2) we haveand substituting it into the first equation of (2) we obtainthat can be rewritten asBy evaluating the time integral of both sides we finally get:whereand coefficient <italic>c<sub>5</sub></italic> accounts for the integration constant. The lead-lag <italic>R</italic>
<sup>2</sup> is the fraction of the total variance explained by model (3). Note that the lead-lag <italic>R</italic>
<sup>2</sup> depends on the time order of the data, whereas the simultaneous <italic>R</italic>
<sup>2</sup> remains the same after a time shuffling of the data. Moreover, it is worth emphasizing that model (3) may well describe other biologically relevant mechanisms, such as time-shifted profiles as shown in the supporting information file ##SUPPL##0##Text S1##. In this case, obviously, the coefficients <italic>c<sub>i</sub></italic>, which depend on the underlying model, will change accordingly. Any pair of time profiles, satisfying model (3) will be said to have a lead-lag relationship and a good fit to (3) can be obtained also in situations different from those assumed to derive it. This property is very useful since it provides flexibility in modeling different biological phenomena resulting from the presence of a common regulatory signal.</p> <p>The reason for the term “lead-lag” is due to the fact that two signals satisfying model (3) also define the transfer function of a “lead-lag compensator” widely used in control systems engineering. Assuming, for the sake of simplicity, the signals devoid of linear trends and noise (<italic>c</italic>
<sub>4</sub> = <italic>c</italic>
<sub>5</sub> = <italic>δ</italic> = 0), model (3) in the Laplace domain is as follows:which can be rewritten as:so that the resulting transfer function between <italic>m<sub>A</sub></italic>(<italic>s</italic>) and <italic>m<sub>B</sub></italic>(<italic>s</italic>) is that of a lead-lag compensator:
</p> <title>A Direct Formula for Computation of the Lead-Lag R<sup>2</sup> from Gene Expression Data</title> <p>Let the available experimental time series of two genes <italic>A</italic> and <italic>B</italic> be composed of <italic>N</italic>&gt;5 samples taken at times <italic>t<sub>1</sub></italic>,…,<italic>t<sub>N</sub></italic>. Model (3)can be rewritten, using matrix notation, as followswhere Δ collects all of the noise terms andso that the least square estimation of the parameter vector isAccordingly, the goodness of fit to model (3) is measured bywhere the norm used is the usual Euclidean norm. It is important to note that the lead-lag R<sup>2</sup> can be computed directly from gene expression data and values near unity indicates that the model well fits the available time series.</p> <title>Numerical computation of time integral</title> <p>Given a gene expression time profile [<italic>mRNA</italic>]<italic><sup>t</sup></italic> measured at times <italic>t<sub>1</sub></italic>,…,<italic>t<sub>N</sub></italic>, we computed its time integral in two steps. First, we used a piecewise cubic Hermite interpolation formula to obtain, for each time interval, 4 more samples. Over the interpolated time series we computed the integral by using a 2-points closed Newton-Cotes formula (trapezoidal rule).</p> <title>Datasets</title> <title>Cell cycle regulated genes</title> <p>We considered the extended list of 1159 cell cycle regulated genes reported in reference ##REF##15513999##[38]##. Each gene in this list has been considered as cell-cycle regulated in at least one of the six methods reported in reference ##REF##15513999##[38]##. We used such an extended list in order to have a sufficiently large dataset for our statistical analysis.</p> <title>Gene expression datasets</title> <p>We considered yeast cell cycle data measured by three independent groups ##REF##9843569##[4]##,##REF##16912276##[43]##,##REF##9702192##[54]##. The data from the Spellman <italic>et al.</italic> group consist of genome-wide gene expression data during the yeast cell cycle using three different synchronization methods. We denoted as ELU, the elutriation based dataset composed of one cell cycle, as ALPHA, the pheromone α arrest factor based dataset composed of two cell cycles and as CDC15 the temperature sensitive CDC15 mutant based dataset composed of three cell cycles. Only two cell cycles of the CDC15 dataset could be used due to the large number of missing data. The dataset in Cho <italic>et al.</italic>
##REF##9702192##[54]##, denoted by CDC28, is composed of two cell cycle and synchronized using a temperature sensistive CDC28 mutant. The last dataset has been downloaded from the authors website ##REF##16912276##[43]## and is composed of three genome-wide gene expression measurement during the yeast cell cycle using alpha factor synchronization. We denoted such dataset, composed of two cell cycles each, as ALPHA_28, ALPHA_30 and ALPHA_38. Two data sets, ALPHA_30 and ALPHA_38, are dye swap technical replicates.</p> <title>Transcription factors dataset</title> <p>We considered the main cell cycle TFs (SWI4, SWI6, MBP1, NDD1, FKH1, FKH2, MCM1, ACE2, SWI5, YOX1) according to Bahler ##REF##16285853##[39]##, and as targets, those genes included in the McIsaac <italic>et al.</italic> dataset ##REF##16522208##[15]## with a stringent threshold for DNA binding (<italic>p</italic>-value&lt;0.001). The MacIsaac <italic>et al.</italic> dataset contained 660 of the 1159 cell cycle regulated genes. Therefore, we ended up with a list of 660 genes available for the subsequent computational analysis.</p> <title>Half-lives dataset</title> <p>We used half-life genome-wide measurements of the yeast transcript measured by Wang <italic>et al.</italic>
##REF##11972065##[26]## and by Kuai <italic>et al.</italic>
##REF##16166263##[45]##.</p> <title>Integration of gene expression datasets</title> <p>For each dataset, we computed the simultaneous and lead-lag–<italic>R</italic>
<sup>2</sup> for all possible pairs using <italic>N</italic> = 660 genes, that is we computed such parameters for <italic>N</italic>(<italic>N</italic>−1)/2 = 217470 pairs. More precisely, the <italic>R</italic>
<sup>2</sup> values were computed for each cell cycle in each dataset, thus obtaining 13 values for each gene pair (ELU: 1 cell cycle, ALPHA: 2 cell cycles, CDC15: 2 cell cycles, CDC28: 2 cell cycles, ALPHA_28: 2 cell cycles, ALPHA_30: 2 cell cycles and ALPHA_38: 2 cell cycles). The average dataset has been constructed by computing the R<sup>2</sup> values for each cycle and for each dataset, for a total amount of 13 cycles. The mean R<sup>2</sup> value for each genes pair was obtained by computing the mean of the 13 available values. In case of missing data in the original dataset, computation of the mean R<sup>2</sup> value was performed only when at least 8 out of 13 cycles were available. Such data were used to compute the diagram showed in ##FIG##2##Figure 3B##. The values obtained by averaging all 13 cycles provided us with a single value for each gene pair and they were used to compute the ROC curve shown in ##FIG##2##Figure 3A##. Cell cycle data with missing values were removed from the dataset.</p>
<title>Results/Discussion</title> <title>Specific Features of Transcript Degradation Regulation Versus Transcription Regulation</title> <p>To clarify the specific features of gene regulation at the mRNA stability level, it is worth thinking of the case when two genes are turned on at the same time by the same transcriptional signal, and the newly synthesized transcripts of both genes are degraded at the same rate. Consequently, differences in their gene expression profile will be determined only by the response of the two genes to the transcriptional signal (<italic>i.e.</italic> different affinities of the transcription factor to promoter regions). A computer simulation of this situation is depicted in ##FIG##0##Figure 1A## where two genes are expressed following a first-order kinetics (see ##SUPPL##0##Text S1##. for details on the equations used for the simulation). The transcription is turned on at the same time for both genes but with a different rate: the first gene is transcribed more rapidly than the second one (##FIG##0##Figure 1Aa##). Their degradation rate is the same (##FIG##0##Figure 1Ab##) and therefore the two gene expression profiles differ only for the magnitude of the response, whereas preserving the shape of the curve (##FIG##0##Figure 1Ac##). In this case, the normalized time profiles are identical (##FIG##0##Figure 1Ad##) and therefore the simultaneous R<sup>2</sup> is maximal (R<sup>2</sup> = 1). Indeed, the “converse” situation is very different. ##FIG##0##Figure 1Ca–b## illustrates the case in which the two genes are transcribed at the same rate while their degradation decreases at the same time but with a different rate. The two profiles do not have the same shape (##FIG##0##Figure 1Cc##). As a consequence, the corresponding simultaneous R<sup>2</sup> will not be maximal (R<sup>2</sup>&lt;1) as can be seen from differences in the normalized profiles (##FIG##0##Figure 1Cd##).</p> <p>Such considerations illustrate that the impact of stability regulation on time profiles is quantitatively and – most importantly – <italic>qualitatively</italic> different from that of transcription regulation. It is therefore not surprising that specific systems biology computational tools have begun to appear in the literature ##REF##17370425##[28]##,##REF##16317069##[32]##. The different impact of mRNA stability regulation versus transcription regulation results from the fact that the rate of mRNA degradation is proportional to the substrate concentration but the rate of production is <italic>not</italic>
##REF##8375615##[33]##. Such behaviour is reasonably well captured by a first order rate equation. In fact, messengers half-lives are experimentally measured usually by fitting a single exponential decay function to the time profiles observed after transcriptional shut-off ##REF##11972065##[26]##.</p> <p>Another important issue is that the differences of transcription rate regulation with respect to degradation rate regulation cannot be clearly seen by simply looking at the long term behavior of the response, <italic>i.e.</italic> at steady state values. In fact, the final amount of mRNA upon a prolonged regulatory signal equals the ratio transcription rate/degradation rate so that, from this perspective, a <italic>N</italic>-fold increase of transcription rate is equivalent to a <italic>N</italic>-fold decrease in degradation rate (and viceversa). An example of such behavior can be seen by comparing ##FIG##0##Figure 1Ac## with ##FIG##0##Figure 1Cc##: the steady state values are the same in both cases but the overall shape of the response (its “dynamics”) is very different.</p> <p>Such “loss of correlation” phenomenon due to differential stability regulation can be further understood by considering a time varying rates, resulting in a transient mRNA time profile, as shown in ##FIG##0##Figure 1Ba–d and 1Da–d##. Again, an increase in the rate of transcription results in an increased response displaying a highly correlated temporal profile (##FIG##0##Figure 1Bd##), whereas a decrease in the degradation rate results in a low correlated temporal profile with a shift in peak time, as shown by ##FIG##0##Figure 1Da–d##. It is plain that, by combining time varying transcription and degradation rates, a large variety of dynamic time patterns may be generated. It is important to note that peak timing regulation may also stem from time delays, as shown in ##FIG##1##Figure 2B##, which can be generated by different biological mechanisms such as transcriptional combinatorial regulation, cascade regulations ##REF##11572776##[34]##, feedforward motifs and single input motifs ##REF##12399590##[35]##. Time delays in gene expression data have been studied using delay correlation analysis ##REF##11743722##[16]##,##UREF##9##[36]##,##REF##15289483##[37]##.</p> <p>The scenario depicted above naturally leads to the possibility that co-regulation may involve both the transcriptional and post-transcriptional machinery. Therefore, a large variety of temporal profiles can be obtained by combining any of those shown in ##FIG##1##Figure 2##.</p> <title>The Lead-Lag Relationship</title> <p>In this paper we consider a novel relationship between gene expression time profiles which includes also the possible presence of mRNA stability variations as a further mechanism to modulate transcript abundance over time. Such new coordinated relationship will be called <italic>lead-lag relationship</italic>. Such terminology is borrowed from the field of system and control engineering where it refers to the basic building block for the realization of a regulatory device able to provide optimal properties to a given process and called “lead-lag compensator” ##UREF##7##[29]##. In order to identify lead-lag relationships, we propose a quantitative measure between gene expression time profiles, called <italic>lead-lag R</italic>
<sup>2</sup>, able to incorporate in a single parameter such relationship and consequently potentially enhancing the predictive power of gene expression analysis for the identification of putatively co-regulated genes. In fact, we aim to study here the possibility that an high lead-lag R<sup>2</sup> between expression time profiles of two given genes is a good indication for the presence of a common regulation mechanism.</p> <p>The lead-lag R<sup>2</sup> is quantitatively defined by a linear multiple regression model among the two given gene expression time profiles <italic>m<sub>A</sub></italic>(<italic>t</italic>) and <italic>m<sub>B</sub></italic>(<italic>t</italic>) and the area under curve until time <italic>t</italic> (<italic>i.e.</italic> their time integral over time):and measured by the <italic>lead-lag</italic> R<sup>2</sup>, that is the fraction of the total variance explained by the above multiple regression model. Such coefficient is computed directly from at least 6 time points of gene expression data and takes values between 0 and 1. The rationale behind such new relationship stems from a simple mathematical model conceived to capture, from gene expression time series data, those genes which are co-regulated at the transcriptional level having an equal or different mRNA stability.</p> <p>It is worth noting that the simultaneous relationship is also a particular lead-lag relationship (just set <italic>c<sub>2</sub></italic> = <italic>c<sub>3</sub></italic> = <italic>c<sub>4</sub></italic> = 0) so that the magnitude of the lead-lag R<sup>2</sup> is always larger or equal than that of the simultaneous R<sup>2</sup>. In the following we will show that the magnitude of the increase from simultaneous R<sup>2</sup> to lead-lag R<sup>2</sup> is specific for each gene pair and that it is statistically correlated both to the presence of a common transcriptional signal and to differences between the half-lives. More details of the lead-lag R<sup>2</sup> and its numerical computation are given in the <xref ref-type="sec" rid="s3">Materials and Methods</xref> section.</p> <title>Predicting Co-Regulation from Lead-Lag Relationships</title> <p>The mathematical model used to define the lead-lag R<sup>2</sup> is based on the assumption that co-regulated genes have the same transcriptional signal (promoter activity) and equal or different transcript stabilities. Consequently, we postulated that two given genes showing a lead-lag relationship (namely, with high lead-lag R<sup>2</sup> values) are likely to be regulated by common transcription factors. To test this hypothesis, we selected a list of 1159 genes indicated as cell-cycle regulated in at least one out of six yeast genome-wide studies ##REF##15513999##[38]##. We then used a large integrated dataset of yeast cell-cycle data generated by three independent groups using different synchronization methods and composed of 7 datasets (13 cell cycles for each gene, see <xref ref-type="sec" rid="s3">Materials and Methods</xref> for details). We considered as “gold standard” the transcriptional regulatory network recently published by MacIsaac and collegues ##REF##16522208##[15]##. Such reconstructed network is very reliable since the authors combined complementary strategies to improve the ability to identify the specificity of transcriptional regulators from genome-wide chromatin immunoprecipitation data. The Mc Isaac <italic>et al.</italic> dataset consists of a list of targets for 203 TFs using different conservative criteria. Among those available 203 TFs, we selected a <italic>p</italic>-value for binding of 0.001 obtaining a list of 3107 genes, containing 660 of the genes in the list of the cell cycle regulated ones. We then choose the 10 TFs widely recognized as having a fundamental role during the cell cycle ##REF##16285853##[39]##: SWI4, SWI6, MBP1, NDD1, FKH1, FKH2, MCM1, ACE2, SWI5 and YOX1. Using this data, we could assess the effectiveness of our approach by computing true and false positive rates and ROC curves. To this end, we evaluated the lead-lag R<sup>2</sup> for each gene pair in the dataset (<italic>N</italic>(<italic>N</italic>−1)/2 pairs, <italic>N</italic> = 660) and considered as putatively co-regulated those pairs whose R<sup>2</sup> values were over a threshold <italic>t<sub>high</sub></italic> and, as putatively non co-regulated, those pairs whose R<sup>2</sup> values were below a threshold <italic>t<sub>low</sub></italic>. Gene pairs with scores between thresholds were not considered. In order to construct a ROC curve we used varying thresholds: as an upper threshold <italic>t<sub>high</sub></italic> for co-regulation we selected the value corresponding to percentiles <italic>p</italic> ranging from 50<sup>th</sup> to 90<sup>th</sup> with a step of 10 and, as a lower threshold <italic>t<sub>low</sub></italic> for non-coregulation, we selected the value corresponding to the “symmetric” percentile 100−<italic>p</italic>. For each threshold we could compute true positives, true negatives, false positives, false negatives and therefore construct a ROC curve (##FIG##2##Figure 3A##, green plots) where all the R<sup>2</sup> values have been averaged over the 7 datasets. The average dataset has been constructed by computing the R<sup>2</sup> values for each cycle and for each dataset, for a total amount of 13 cycles. The mean R<sup>2</sup> value for each genes pair was obtained by computing the mean of the 13 available values. In case of missing data in the original dataset, computation of the mean R<sup>2</sup> value was performed only when at least 8 out of 13 cycles were available. Each class of putatively co-regulated gene pairs was obtained by selecting those pairs exceeding the upper thresholds corresponding to the percentiles from 50<sup>th</sup> to 90<sup>th</sup> with a step of 10 of the R<sup>2</sup> distribution. Therefore, true positives are those pairs of the class having at least one common transcription factor according to the Mc Isaac <italic>et al.</italic> dataset (<italic>p</italic>-value for binding&lt;0.001), whereas false positives are those pairs in the class without a common transcription factor (<italic>p</italic>-value for binding&gt;0.001). Analogously, true negatives and false negatives were computed within the class of gene pairs having the lower thresholds corresponding to the percentiles from 50<sup>th</sup> to 10<sup>th</sup> with a step of −10 of the R<sup>2</sup> distribution.</p> <p>To evaluate the performance of predictions obtained with the lead-lag R<sup>2</sup> we repeated the same analysis using the simultaneous R<sup>2</sup> as a similarity measure between two given genes (##FIG##2##Figure 3A## blue plots). The results clearly show that the lead-lag R<sup>2</sup> certainly outperforms standard analysis based only on simultaneous relationships, increasing the true positive rate from 80% to 95%. The fraction of false positives also slightly decrease but remains relatively high possibly due to the fact that we have considered highly conservative criteria and selected the targets of only those 10 transcription factors having a major regulatory role during the cell cycle. It is worth noting that the performances on the average dataset are much better than the average of all performances (see ##FIG##2##Figure 3B##) thus showing that an integrative approach using multiple independent datasets is always the best choice, whenever applicable. Moreover, given the large number of datasets considered, we can also conclude that the results obtained are largely independent of the noise and the stress response induced by the synchronization methods. Finally, we note that the above results remain valid even if we consider the transcriptional network presented by McIsaac <italic>et al.</italic>
##REF##16522208##[15]## using different selection criteria for DNA binding (see ##SUPPL##0##Text S1##).</p> <title>mRNA Half-Lives and Lead-Lag R<sup>2</sup>
</title> <p>The peculiarity of the lead-lag relationship between two given genes relies on the presence of a common regulatory signal driving the expression of transcripts with equal or different mRNA half-lives. For this reason, we investigated whether co-regulated gene pairs having an high lead-lag R<sup>2</sup> values are significantly enriched with differential transcript's stabilities. Half-life values are not available during the cell cycle and in the same experimental conditions used for establishing cell synchronization. Nevertheless, genome-wide half-lives data for un-synchronized cells were published recently by Wang <italic>et al.</italic>
##REF##11972065##[26]##. Using DNA microarrays, the authors precisely measured the decay of each yeast mRNA in YPD medium, after thermal inactivation of a temperature-sensitive RNA polymerase II. Such half-life measurements were not obtained during the cell cycle, so that we do not expect an exact agreement with the actual ones. Nevertheless, by considering a large number of gene pairs (16740) it appears reasonable that, on average, the half-life ratios between gene pairs may not vary significantly. Therefore, we used such available data for a statistical evaluation of the presence of gene pairs with high lead-lag R<sup>2</sup> values with respect to the simultaneous R<sup>2</sup> among those co-regulated pairs having large half-life ratios.</p> <p>To this end, we considered all possible gene pairs having, at least, one common transcription factor according to the MacIsaac <italic>et al.</italic> dataset ##REF##16522208##[15]## using a <italic>p</italic>-value for binding less that 0.001 and considered five half-life ratio bins: less than 2-fold, from 2-fold to 3-fold, from 3-fold to 4-fold, from 4-fold to 5-fold and more than 5-fold. We computed the simultaneous R<sup>2</sup> and also the difference between the lead-lag R<sup>2</sup> and simultaneous R<sup>2</sup> for all the gene pairs in each of the half-life bins. Such difference is used in order to select that part of the lead-lag R<sup>2</sup> value which is not due to the simultaneous espression of the gene pair. Therefore, we got a distribution of values for each half-life ratio bin and computed the corresponding mean value and standard deviation.</p> <p>
##FIG##3##Figure 4## shows the results of the above described computation. ##FIG##3##Figure 4A## makes clear that the highest and the lowest half-life bin display very different lead-lag minus simultaneous R<sup>2</sup> mean values. To further support this feature, we performed a <italic>t</italic> test and found that the increase of the mean value of the distribution in the first and the last bin of the lead-lag R<sup>2</sup> minus the simultaneous R<sup>2</sup> is indeed significant (95% confidence level, <italic>p</italic>-value 10<sup>−9</sup>). The simultaneous R<sup>2</sup> also shows a mildly significant decrease (95% confidence level, <italic>p</italic>-value 0.03) of the mean values between the first and the last bin. To further evaluate the statistical significance of this analysis we computed the <italic>Z</italic>-score corresponding to 100000 randomizations of the half-life measurements. The results are shown in the scatterplot of ##FIG##3##Figure 4B## and they provide computational evidence that the lead-lag R<sup>2</sup> of gene pairs is statistically correlated to their half-life ratios. In fact, a high positive <italic>Z</italic>-score (about 5) corresponds to the highest half-life ratio bin and a negative <italic>Z</italic>-score (about −5) corresponds to the first half-life ratio bin. On the other hand, <italic>Z</italic> scores for the simultaneous R<sup>2</sup> are all within the values −3 and 3 and therefore the observed difference of the mean values between the first and the last bin is not significantly affected by the randomizations. This scenario is consistent with the biological process underlying the mathematical model used to define the lead-lag R<sup>2</sup> thus showing that our analysis well captures the effects of post-transcriptional control on gene expression time profiles during the cell-cycle.</p> <title>Comparison to Other Similarity Measures</title> <p>The results presented so far have clearly shown that lead-lag correlation analysis outperforms the usual simultaneous correlation analysis (squared Pearson coefficient) for the prediction of co-regulation, <italic>i.e.</italic> the presence of a common transcription factor, from gene expression time profiles. As previously discussed, truly co-regulated genes do often display large differences of gene expression time profiles, <italic>e.g.</italic> peak shifts, delays or other kinds of nonlinear relationships. In this paragraph, we consider other similarity measures relevant to the analysis of gene expression data and compare their performances with those obtained using the lead-lag R<sup>2</sup>. In particular, we used 5 similarity measures other than the lead-lag: Spearman's rank, Kendall's tau, cosine, dynamic time-warped and time-delayed correlation, all squared to capture inverted relationships also. Spearman's rank, Kendall's tau and cosine correlation are the most common choices for the analysis of gene expression data in the presence of nonlinear relationships between time series, but they do not take into account the time ordering of data. By contrast, time-warped and time-delayed correlation have been specifically developed to analyze gene expression time profiles. The time-delayed correlation analysis has been proposed by Schmitt <italic>et al.</italic>
##REF##15289483##[37]## where, for any genes pair, a R<sup>2</sup> value is obtained by selecting the highest simultaneous R<sup>2</sup> over all admissible time delays between profiles. The dynamic time-warped correlation has been recently used by Aach and Church ##REF##11395426##[40]## and Hermans and Tsiporkova ##REF##17060360##[41]## for the alignment of gene expression time series obtained in experiments using different cell synchronization methods. These two works are both based, for gene-to-gene comparisons, on the Dynamic Time Warping (DTW) algorithm developed by Sankoff and Kruskal ##UREF##10##[42]##. Accordingly, we defined a time-warped R<sup>2</sup> by selecting the highest simultaneous R<sup>2</sup> over all the possible time warped paths. For any similarity measure, we performed the same analysis reported in a previous section using the same data, and the results are shown in ##FIG##4##Figure 5## where sensitivity (panel A) and specificity (panel B) are reported for each threshold.</p> <p>First of all, the cosine correlation analysis produces the poorest performances, very close to a random choice, and therefore such similarity measure is not reported in ##FIG##4##Figure 5##. On both panels we note that simultaneous, Spearman's rank and Kendall's tau provide comparable results which are clearly below the performances of the group of methods which take into account the time ordering of data. In this group the lead-lag correlation analysis shows the best performances, both in terms of sensitivity and specificity.</p> <title>Examples of Lead-Lag Analysis Using Yeast Cell Cycle Gene Expression Data</title> <p>In this section we present some examples of “typical” lead-lag relationships using the most recent yeast cell cycle data ##REF##16912276##[43]## and discuss their biological relevance. The complete list of gene pairs exceeding the 95<sup>th</sup> percentile of the distribution for each of the R<sup>2</sup> values considered in this paper is provided in the supporting information file ##SUPPL##0##Text S1##.</p> <title>Key cell cycle regulators under common transcription factors</title> <p>The budding yeast cell cycle is characterized by consecutive waves of expression of key regulators such as cyclins and transcription factors ##UREF##11##[44]##. CLB6, a G1/S-phase cyclin, has a lead-lag relationship with GIN4 as shown in ##FIG##5##Figure 6A##, a gene encoding a key component involved in transitioning to the next stage of the cycle ##REF##11572776##[34]##. The lead-lag relationship suggests the presence of a common transcription factor and, consistently, the two genes are both targets of the transcription factor complex MBF/SBF according to ChIP experiments ##REF##16522208##[15]##. Moreover, the time profiles shown in ##FIG##5##Figure 6A## indicate also the possibility that the transcriptional signal is turned on and then quickly turned off, so that the subsequent behaviour of the two genes is mainly determined by the degradation process alone. Accordingly, transcripts stabilities – as measured after transcriptional shut-off ##REF##16166263##[45]## – significantly differ in value.</p> <p>Cell Division Cycle 6 (CDC6) is a component of the pre-replicative complex essential for the initiation of DNA replication, normally expressed at the end of mitosis. It has a lead-lag relationships with ASH1 (##FIG##5##Figure 6B##) which encodes a GATA-like transcription factor localized at daughter cells where it serves to repress the late G1-specific transcription of HO and preventing mating-type switching ##UREF##12##[46]##. Our analysis suggests the presence of a common regulatory signal. In fact, both genes are key regulators of separate biological processes that are simultaneously activated by the SWI5 transcription factor ##REF##16522208##[15]##,##REF##7641697##[47]##,##REF##10409653##[48]##. Moreover, our analysis also suggests that the CDC6 transcript is fairly unstable. Consistently, the CDC6 protein is unstable ##REF##7641697##[47]##.</p> <p>SWI5 encodes a key transcription factor that activates transcription of genes expressed at the M/G1 boundary and in G1 phase of the cell cycle. NCE102 is a non-classical export protein involved in alternative clearance/detoxification pathway to eliminate damaged material ##UREF##13##[49]##. They display a lead-lag relationship (##FIG##5##Figure 6C##) and, in fact, they are both targets of the M-phase activating complex FKH2/NDD1 according to ChIP experiments ##REF##16522208##[15]## and large differences of their half-life values are observed after transcription inhibition ##REF##11972065##[26]##. The gene expression profiles shown in ##FIG##5##Figure 6C## reflect the prototypical situation of peak delay depicted in ##FIG##1##Figure 2C##.</p> <p>YOX1 is a transcription factor involved in the repression of ECB acitivity ##UREF##12##[46]## thus contributing to move the cycle forward. YOX1 shows a lead-lag relationship with MNN1 (##FIG##5##Figure 6D##), a gene encoding a cell wall glycoprotein ##REF##16672383##[50]##. Consistently, they are both targets of SBF according to ChIP experiments ##REF##16522208##[15]## and the transcripts have different half-life values ##REF##11972065##[26]##. Moreover, looking at the time profiles depicted in ##FIG##5##Figure 6D##, one may argue that during the second half of the cycle another transcription factor is active at the MNN1 promoter.</p> <p>All the above examples consist of pairs of genes that are under the control of the same transcription factor and that show differential mRNA stability values consistent with their lead-lag relationship (except for CDC6 transcript whose experimental half-life is not available). Moreover, it is worth noting that large differences in half-lives value (as in the cases shown in ##FIG##5##figure 6C and 6D##) significantly affect the overall time profiles producing also an evident peak shift.</p> <p>Finally, it is worth noting that the lead-lag relationship is symmetrical and, therefore, it does not provide information about which gene is “lead” and which is “lag”. However, such information can be easily obtained by visual inspection. In fact, from ##FIG##5##Figure 6A##, the lead gene is the one with the steepest decaying profile having, consistently, a smaller half-life. Moreover, from ##FIG##5##Figure 6C##, one can see that the gene with the larger half-life displays a delayed peak and, therefore, it corresponds to the lag gene.</p> <title>Dynamic formation of the replication complex</title> <p>Many studies have focused on the relationship between gene expression time courses and the formation of protein complexes. Interestingly, Jansen <italic>et al.</italic>
##REF##11779829##[31]## suggested to classify protein complexes as either <italic>permanent</italic> or <italic>transient</italic>, with permament ones being maintained through most cellular conditions. They also found that, generally, permanent complexes tend to have simultaneously correlated gene expression while transient ones do not. Moreover, they also noted that subunits of the same protein complex may show significant simultaneous expression. In particular, they studied gene expression of the replication complex in yeast and found a very low simultaneous correlation among subunits, not significantly different from a random control ##REF##11779829##[31]##. However, they also found two sub-complexes – the MCM complex and the DNA polymerases δ and ε complex – showing much greater simultaneous correlation.</p> <p>Using gene expression time profiles during one cell cycle (##REF##16912276##[43]##, dataset, alpha_38 time series) for the genes encoding MCM proteins (MCM cluster) and DNA polymerases and ε (POL cluster), we computed simultaneous and lead-lag R<sup>2</sup> and the scatterplots of the resulting values for gene pairs belonging to the two different sub-complexes are shown in ##FIG##6##Figure 7##, panel B. As a negative control we used a group of 5 simultaneously expressed genes (R<sup>2</sup>&gt;0.7) coding for proteins of the cytoplasmic ribosomal large subunit (RPL4A, RPL4B, RPL1A, RPL1B, RPP0 denoted by RIB cluster). Ribosomal proteins are under the transcriptional control of IFH1/FHL1 ##REF##15128447##[51]##,##REF##15616568##[52]## whereas the replication complex is regulated by the transcription factors MBF/SBF ##REF##15128447##[51]##. The scatterplot reported in ##FIG##6##Figure 7## of the simultaneous <italic>vs.</italic> lead-lag R<sup>2</sup> values shows that, whereas the POL/MCM pairs display high values of lead-lag R<sup>2</sup> and low values of simultaneous R<sup>2</sup>, the control pairs POL/RIB and MCM/RIB display a very different pattern spread over a larger range thus denoting the absence of any meaningful relationship.</p> <p>
##FIG##6##Figure 7## makes very clear that the gene expression of the MCM and the DNA polymerases δ and ε subcomplexes is significantly simultaneously correlated within the same group whereas such correlation dramatically drops if we consider pairs of genes belonging to different subcomplexes. In fact, the sample distribution of the simultaneous R<sup>2</sup> between the two clusters is spread over the range [0,0.5] thus showing the absence of any significant level of simultaneous correlation. By contrast, the between clusters lead-lag R<sup>2</sup> histogram is concentrated in the highest part of the range close to 1. The high values of the lead-lag R<sup>2</sup> strongly suggest that mRNA stability may play a fundamental role in the dynamic formation of multiple protein complexes. Accordingly, the average half-life measured after transcriptional shut-off of the MCM group is 14±6 min (the “lead” genes). and that of the POL group is 19±6 min ##REF##11972065##[26]## (the “lag” genes). The presence of lead-lag relationships between transient sub-complexes is briefly discussed in the supporting information file ##SUPPL##0##Text S1##.</p> <title>Conclusions</title> <p>The expression of genes in the cell is to a large extent controlled at the level of mRNA accumulation. One key point in the analysis of gene expression dynamics is that mRNA abundance is determined by two regulated processes: transcription and degradation both specifically affecting transcript levels. Computational analysis of genome-wide expression time series has shown that clusters of co-expressed (<italic>i.e.</italic> simultaneously correlated) profiles often provide clues for the presence of common transcription factors regulating both genes. Such computational analysis (known as “clustering”) is very useful since it allows the prediction of the underlying regulatory actions based exclusively on the available gene expression data obtained from a given experiment. The rationale behind such belief is a sort of a “guilty by association” approach: genes' products appearing and disappearing at the same time are likely to have some common transcriptional regulation. Nevertheless, it may well be the case that the same transcriptional signal regulating two (or more) genes may yield quite different outcomes on each transcript. In fact, a number of biological events following transcription may selectively affect cytoplasmic mRNA abundance, such as, for example, the activity of the enzymatic machinery involved in mRNA processing and degradation. In order to address this issue, we provided a novel computational methodology that, based exclusively on the available gene expression data, is able to effectively predict co-regulation even with variation in the dynamic response due to mRNA stability differences. Moreover, our approach also captures the relation of simultaneous or time shifted co-expression so that it provides a single integrative general index – the <italic>lead-lag R<sup>2</sup></italic>−able to uncover the presence of a common regulatory signal underlying gene expression time dynamics also at the post-transcriptional level.</p> <p>In order to test the validity of our approach on real data, we used yeast genome-wide cell-cycle expression time series obtained by several independent groups using different synchronization methods. In fact, by doing so, we could integrate the available cell cycle data and obtain a much more reliable aggregated dataset. We considered those gene pairs with the highest lead-lag R<sup>2</sup> values and found the prediction for the presence of a common transcription factor to be highly consistent with protein-DNA binding data (ChIP experiments). Our results clearly indicate that co-regulation is not generally equivalent to simultaneous expression.</p> <p>We believe that the same analysis can be successfully used to predict post-transcriptional regulation, <italic>i.e.</italic> the presence of a common mechanism able to stabilize or de-stabilize specific transcripts, as for the members of the PUF proteins family ##REF##15024427##[2]##. Moreover, we envisage the possibility that our methodology could be used on different data and organisms and thus providing a computational support to the understanding of transcriptional and post-transcriptional networks, given the recent growing interest in the post-transcriptional regulation layer ##REF##17245413##[1]## of gene expression (miRNA) and its role in many diseases, such as cancer. Finally, the characterization of the replication complex in terms of lead-lag relationships among gene expression time profiles of its sub-complexes suggests the possibility that our analysis could be effectively used as a tool for predicting the formation of transient multiple protein complexes.</p>
<title>Results/Discussion</title> <title>Specific Features of Transcript Degradation Regulation Versus Transcription Regulation</title> <p>To clarify the specific features of gene regulation at the mRNA stability level, it is worth thinking of the case when two genes are turned on at the same time by the same transcriptional signal, and the newly synthesized transcripts of both genes are degraded at the same rate. Consequently, differences in their gene expression profile will be determined only by the response of the two genes to the transcriptional signal (<italic>i.e.</italic> different affinities of the transcription factor to promoter regions). A computer simulation of this situation is depicted in ##FIG##0##Figure 1A## where two genes are expressed following a first-order kinetics (see ##SUPPL##0##Text S1##. for details on the equations used for the simulation). The transcription is turned on at the same time for both genes but with a different rate: the first gene is transcribed more rapidly than the second one (##FIG##0##Figure 1Aa##). Their degradation rate is the same (##FIG##0##Figure 1Ab##) and therefore the two gene expression profiles differ only for the magnitude of the response, whereas preserving the shape of the curve (##FIG##0##Figure 1Ac##). In this case, the normalized time profiles are identical (##FIG##0##Figure 1Ad##) and therefore the simultaneous R<sup>2</sup> is maximal (R<sup>2</sup> = 1). Indeed, the “converse” situation is very different. ##FIG##0##Figure 1Ca–b## illustrates the case in which the two genes are transcribed at the same rate while their degradation decreases at the same time but with a different rate. The two profiles do not have the same shape (##FIG##0##Figure 1Cc##). As a consequence, the corresponding simultaneous R<sup>2</sup> will not be maximal (R<sup>2</sup>&lt;1) as can be seen from differences in the normalized profiles (##FIG##0##Figure 1Cd##).</p> <p>Such considerations illustrate that the impact of stability regulation on time profiles is quantitatively and – most importantly – <italic>qualitatively</italic> different from that of transcription regulation. It is therefore not surprising that specific systems biology computational tools have begun to appear in the literature ##REF##17370425##[28]##,##REF##16317069##[32]##. The different impact of mRNA stability regulation versus transcription regulation results from the fact that the rate of mRNA degradation is proportional to the substrate concentration but the rate of production is <italic>not</italic>
##REF##8375615##[33]##. Such behaviour is reasonably well captured by a first order rate equation. In fact, messengers half-lives are experimentally measured usually by fitting a single exponential decay function to the time profiles observed after transcriptional shut-off ##REF##11972065##[26]##.</p> <p>Another important issue is that the differences of transcription rate regulation with respect to degradation rate regulation cannot be clearly seen by simply looking at the long term behavior of the response, <italic>i.e.</italic> at steady state values. In fact, the final amount of mRNA upon a prolonged regulatory signal equals the ratio transcription rate/degradation rate so that, from this perspective, a <italic>N</italic>-fold increase of transcription rate is equivalent to a <italic>N</italic>-fold decrease in degradation rate (and viceversa). An example of such behavior can be seen by comparing ##FIG##0##Figure 1Ac## with ##FIG##0##Figure 1Cc##: the steady state values are the same in both cases but the overall shape of the response (its “dynamics”) is very different.</p> <p>Such “loss of correlation” phenomenon due to differential stability regulation can be further understood by considering a time varying rates, resulting in a transient mRNA time profile, as shown in ##FIG##0##Figure 1Ba–d and 1Da–d##. Again, an increase in the rate of transcription results in an increased response displaying a highly correlated temporal profile (##FIG##0##Figure 1Bd##), whereas a decrease in the degradation rate results in a low correlated temporal profile with a shift in peak time, as shown by ##FIG##0##Figure 1Da–d##. It is plain that, by combining time varying transcription and degradation rates, a large variety of dynamic time patterns may be generated. It is important to note that peak timing regulation may also stem from time delays, as shown in ##FIG##1##Figure 2B##, which can be generated by different biological mechanisms such as transcriptional combinatorial regulation, cascade regulations ##REF##11572776##[34]##, feedforward motifs and single input motifs ##REF##12399590##[35]##. Time delays in gene expression data have been studied using delay correlation analysis ##REF##11743722##[16]##,##UREF##9##[36]##,##REF##15289483##[37]##.</p> <p>The scenario depicted above naturally leads to the possibility that co-regulation may involve both the transcriptional and post-transcriptional machinery. Therefore, a large variety of temporal profiles can be obtained by combining any of those shown in ##FIG##1##Figure 2##.</p> <title>The Lead-Lag Relationship</title> <p>In this paper we consider a novel relationship between gene expression time profiles which includes also the possible presence of mRNA stability variations as a further mechanism to modulate transcript abundance over time. Such new coordinated relationship will be called <italic>lead-lag relationship</italic>. Such terminology is borrowed from the field of system and control engineering where it refers to the basic building block for the realization of a regulatory device able to provide optimal properties to a given process and called “lead-lag compensator” ##UREF##7##[29]##. In order to identify lead-lag relationships, we propose a quantitative measure between gene expression time profiles, called <italic>lead-lag R</italic>
<sup>2</sup>, able to incorporate in a single parameter such relationship and consequently potentially enhancing the predictive power of gene expression analysis for the identification of putatively co-regulated genes. In fact, we aim to study here the possibility that an high lead-lag R<sup>2</sup> between expression time profiles of two given genes is a good indication for the presence of a common regulation mechanism.</p> <p>The lead-lag R<sup>2</sup> is quantitatively defined by a linear multiple regression model among the two given gene expression time profiles <italic>m<sub>A</sub></italic>(<italic>t</italic>) and <italic>m<sub>B</sub></italic>(<italic>t</italic>) and the area under curve until time <italic>t</italic> (<italic>i.e.</italic> their time integral over time):and measured by the <italic>lead-lag</italic> R<sup>2</sup>, that is the fraction of the total variance explained by the above multiple regression model. Such coefficient is computed directly from at least 6 time points of gene expression data and takes values between 0 and 1. The rationale behind such new relationship stems from a simple mathematical model conceived to capture, from gene expression time series data, those genes which are co-regulated at the transcriptional level having an equal or different mRNA stability.</p> <p>It is worth noting that the simultaneous relationship is also a particular lead-lag relationship (just set <italic>c<sub>2</sub></italic> = <italic>c<sub>3</sub></italic> = <italic>c<sub>4</sub></italic> = 0) so that the magnitude of the lead-lag R<sup>2</sup> is always larger or equal than that of the simultaneous R<sup>2</sup>. In the following we will show that the magnitude of the increase from simultaneous R<sup>2</sup> to lead-lag R<sup>2</sup> is specific for each gene pair and that it is statistically correlated both to the presence of a common transcriptional signal and to differences between the half-lives. More details of the lead-lag R<sup>2</sup> and its numerical computation are given in the <xref ref-type="sec" rid="s3">Materials and Methods</xref> section.</p> <title>Predicting Co-Regulation from Lead-Lag Relationships</title> <p>The mathematical model used to define the lead-lag R<sup>2</sup> is based on the assumption that co-regulated genes have the same transcriptional signal (promoter activity) and equal or different transcript stabilities. Consequently, we postulated that two given genes showing a lead-lag relationship (namely, with high lead-lag R<sup>2</sup> values) are likely to be regulated by common transcription factors. To test this hypothesis, we selected a list of 1159 genes indicated as cell-cycle regulated in at least one out of six yeast genome-wide studies ##REF##15513999##[38]##. We then used a large integrated dataset of yeast cell-cycle data generated by three independent groups using different synchronization methods and composed of 7 datasets (13 cell cycles for each gene, see <xref ref-type="sec" rid="s3">Materials and Methods</xref> for details). We considered as “gold standard” the transcriptional regulatory network recently published by MacIsaac and collegues ##REF##16522208##[15]##. Such reconstructed network is very reliable since the authors combined complementary strategies to improve the ability to identify the specificity of transcriptional regulators from genome-wide chromatin immunoprecipitation data. The Mc Isaac <italic>et al.</italic> dataset consists of a list of targets for 203 TFs using different conservative criteria. Among those available 203 TFs, we selected a <italic>p</italic>-value for binding of 0.001 obtaining a list of 3107 genes, containing 660 of the genes in the list of the cell cycle regulated ones. We then choose the 10 TFs widely recognized as having a fundamental role during the cell cycle ##REF##16285853##[39]##: SWI4, SWI6, MBP1, NDD1, FKH1, FKH2, MCM1, ACE2, SWI5 and YOX1. Using this data, we could assess the effectiveness of our approach by computing true and false positive rates and ROC curves. To this end, we evaluated the lead-lag R<sup>2</sup> for each gene pair in the dataset (<italic>N</italic>(<italic>N</italic>−1)/2 pairs, <italic>N</italic> = 660) and considered as putatively co-regulated those pairs whose R<sup>2</sup> values were over a threshold <italic>t<sub>high</sub></italic> and, as putatively non co-regulated, those pairs whose R<sup>2</sup> values were below a threshold <italic>t<sub>low</sub></italic>. Gene pairs with scores between thresholds were not considered. In order to construct a ROC curve we used varying thresholds: as an upper threshold <italic>t<sub>high</sub></italic> for co-regulation we selected the value corresponding to percentiles <italic>p</italic> ranging from 50<sup>th</sup> to 90<sup>th</sup> with a step of 10 and, as a lower threshold <italic>t<sub>low</sub></italic> for non-coregulation, we selected the value corresponding to the “symmetric” percentile 100−<italic>p</italic>. For each threshold we could compute true positives, true negatives, false positives, false negatives and therefore construct a ROC curve (##FIG##2##Figure 3A##, green plots) where all the R<sup>2</sup> values have been averaged over the 7 datasets. The average dataset has been constructed by computing the R<sup>2</sup> values for each cycle and for each dataset, for a total amount of 13 cycles. The mean R<sup>2</sup> value for each genes pair was obtained by computing the mean of the 13 available values. In case of missing data in the original dataset, computation of the mean R<sup>2</sup> value was performed only when at least 8 out of 13 cycles were available. Each class of putatively co-regulated gene pairs was obtained by selecting those pairs exceeding the upper thresholds corresponding to the percentiles from 50<sup>th</sup> to 90<sup>th</sup> with a step of 10 of the R<sup>2</sup> distribution. Therefore, true positives are those pairs of the class having at least one common transcription factor according to the Mc Isaac <italic>et al.</italic> dataset (<italic>p</italic>-value for binding&lt;0.001), whereas false positives are those pairs in the class without a common transcription factor (<italic>p</italic>-value for binding&gt;0.001). Analogously, true negatives and false negatives were computed within the class of gene pairs having the lower thresholds corresponding to the percentiles from 50<sup>th</sup> to 10<sup>th</sup> with a step of −10 of the R<sup>2</sup> distribution.</p> <p>To evaluate the performance of predictions obtained with the lead-lag R<sup>2</sup> we repeated the same analysis using the simultaneous R<sup>2</sup> as a similarity measure between two given genes (##FIG##2##Figure 3A## blue plots). The results clearly show that the lead-lag R<sup>2</sup> certainly outperforms standard analysis based only on simultaneous relationships, increasing the true positive rate from 80% to 95%. The fraction of false positives also slightly decrease but remains relatively high possibly due to the fact that we have considered highly conservative criteria and selected the targets of only those 10 transcription factors having a major regulatory role during the cell cycle. It is worth noting that the performances on the average dataset are much better than the average of all performances (see ##FIG##2##Figure 3B##) thus showing that an integrative approach using multiple independent datasets is always the best choice, whenever applicable. Moreover, given the large number of datasets considered, we can also conclude that the results obtained are largely independent of the noise and the stress response induced by the synchronization methods. Finally, we note that the above results remain valid even if we consider the transcriptional network presented by McIsaac <italic>et al.</italic>
##REF##16522208##[15]## using different selection criteria for DNA binding (see ##SUPPL##0##Text S1##).</p> <title>mRNA Half-Lives and Lead-Lag R<sup>2</sup>
</title> <p>The peculiarity of the lead-lag relationship between two given genes relies on the presence of a common regulatory signal driving the expression of transcripts with equal or different mRNA half-lives. For this reason, we investigated whether co-regulated gene pairs having an high lead-lag R<sup>2</sup> values are significantly enriched with differential transcript's stabilities. Half-life values are not available during the cell cycle and in the same experimental conditions used for establishing cell synchronization. Nevertheless, genome-wide half-lives data for un-synchronized cells were published recently by Wang <italic>et al.</italic>
##REF##11972065##[26]##. Using DNA microarrays, the authors precisely measured the decay of each yeast mRNA in YPD medium, after thermal inactivation of a temperature-sensitive RNA polymerase II. Such half-life measurements were not obtained during the cell cycle, so that we do not expect an exact agreement with the actual ones. Nevertheless, by considering a large number of gene pairs (16740) it appears reasonable that, on average, the half-life ratios between gene pairs may not vary significantly. Therefore, we used such available data for a statistical evaluation of the presence of gene pairs with high lead-lag R<sup>2</sup> values with respect to the simultaneous R<sup>2</sup> among those co-regulated pairs having large half-life ratios.</p> <p>To this end, we considered all possible gene pairs having, at least, one common transcription factor according to the MacIsaac <italic>et al.</italic> dataset ##REF##16522208##[15]## using a <italic>p</italic>-value for binding less that 0.001 and considered five half-life ratio bins: less than 2-fold, from 2-fold to 3-fold, from 3-fold to 4-fold, from 4-fold to 5-fold and more than 5-fold. We computed the simultaneous R<sup>2</sup> and also the difference between the lead-lag R<sup>2</sup> and simultaneous R<sup>2</sup> for all the gene pairs in each of the half-life bins. Such difference is used in order to select that part of the lead-lag R<sup>2</sup> value which is not due to the simultaneous espression of the gene pair. Therefore, we got a distribution of values for each half-life ratio bin and computed the corresponding mean value and standard deviation.</p> <p>
##FIG##3##Figure 4## shows the results of the above described computation. ##FIG##3##Figure 4A## makes clear that the highest and the lowest half-life bin display very different lead-lag minus simultaneous R<sup>2</sup> mean values. To further support this feature, we performed a <italic>t</italic> test and found that the increase of the mean value of the distribution in the first and the last bin of the lead-lag R<sup>2</sup> minus the simultaneous R<sup>2</sup> is indeed significant (95% confidence level, <italic>p</italic>-value 10<sup>−9</sup>). The simultaneous R<sup>2</sup> also shows a mildly significant decrease (95% confidence level, <italic>p</italic>-value 0.03) of the mean values between the first and the last bin. To further evaluate the statistical significance of this analysis we computed the <italic>Z</italic>-score corresponding to 100000 randomizations of the half-life measurements. The results are shown in the scatterplot of ##FIG##3##Figure 4B## and they provide computational evidence that the lead-lag R<sup>2</sup> of gene pairs is statistically correlated to their half-life ratios. In fact, a high positive <italic>Z</italic>-score (about 5) corresponds to the highest half-life ratio bin and a negative <italic>Z</italic>-score (about −5) corresponds to the first half-life ratio bin. On the other hand, <italic>Z</italic> scores for the simultaneous R<sup>2</sup> are all within the values −3 and 3 and therefore the observed difference of the mean values between the first and the last bin is not significantly affected by the randomizations. This scenario is consistent with the biological process underlying the mathematical model used to define the lead-lag R<sup>2</sup> thus showing that our analysis well captures the effects of post-transcriptional control on gene expression time profiles during the cell-cycle.</p> <title>Comparison to Other Similarity Measures</title> <p>The results presented so far have clearly shown that lead-lag correlation analysis outperforms the usual simultaneous correlation analysis (squared Pearson coefficient) for the prediction of co-regulation, <italic>i.e.</italic> the presence of a common transcription factor, from gene expression time profiles. As previously discussed, truly co-regulated genes do often display large differences of gene expression time profiles, <italic>e.g.</italic> peak shifts, delays or other kinds of nonlinear relationships. In this paragraph, we consider other similarity measures relevant to the analysis of gene expression data and compare their performances with those obtained using the lead-lag R<sup>2</sup>. In particular, we used 5 similarity measures other than the lead-lag: Spearman's rank, Kendall's tau, cosine, dynamic time-warped and time-delayed correlation, all squared to capture inverted relationships also. Spearman's rank, Kendall's tau and cosine correlation are the most common choices for the analysis of gene expression data in the presence of nonlinear relationships between time series, but they do not take into account the time ordering of data. By contrast, time-warped and time-delayed correlation have been specifically developed to analyze gene expression time profiles. The time-delayed correlation analysis has been proposed by Schmitt <italic>et al.</italic>
##REF##15289483##[37]## where, for any genes pair, a R<sup>2</sup> value is obtained by selecting the highest simultaneous R<sup>2</sup> over all admissible time delays between profiles. The dynamic time-warped correlation has been recently used by Aach and Church ##REF##11395426##[40]## and Hermans and Tsiporkova ##REF##17060360##[41]## for the alignment of gene expression time series obtained in experiments using different cell synchronization methods. These two works are both based, for gene-to-gene comparisons, on the Dynamic Time Warping (DTW) algorithm developed by Sankoff and Kruskal ##UREF##10##[42]##. Accordingly, we defined a time-warped R<sup>2</sup> by selecting the highest simultaneous R<sup>2</sup> over all the possible time warped paths. For any similarity measure, we performed the same analysis reported in a previous section using the same data, and the results are shown in ##FIG##4##Figure 5## where sensitivity (panel A) and specificity (panel B) are reported for each threshold.</p> <p>First of all, the cosine correlation analysis produces the poorest performances, very close to a random choice, and therefore such similarity measure is not reported in ##FIG##4##Figure 5##. On both panels we note that simultaneous, Spearman's rank and Kendall's tau provide comparable results which are clearly below the performances of the group of methods which take into account the time ordering of data. In this group the lead-lag correlation analysis shows the best performances, both in terms of sensitivity and specificity.</p> <title>Examples of Lead-Lag Analysis Using Yeast Cell Cycle Gene Expression Data</title> <p>In this section we present some examples of “typical” lead-lag relationships using the most recent yeast cell cycle data ##REF##16912276##[43]## and discuss their biological relevance. The complete list of gene pairs exceeding the 95<sup>th</sup> percentile of the distribution for each of the R<sup>2</sup> values considered in this paper is provided in the supporting information file ##SUPPL##0##Text S1##.</p> <title>Key cell cycle regulators under common transcription factors</title> <p>The budding yeast cell cycle is characterized by consecutive waves of expression of key regulators such as cyclins and transcription factors ##UREF##11##[44]##. CLB6, a G1/S-phase cyclin, has a lead-lag relationship with GIN4 as shown in ##FIG##5##Figure 6A##, a gene encoding a key component involved in transitioning to the next stage of the cycle ##REF##11572776##[34]##. The lead-lag relationship suggests the presence of a common transcription factor and, consistently, the two genes are both targets of the transcription factor complex MBF/SBF according to ChIP experiments ##REF##16522208##[15]##. Moreover, the time profiles shown in ##FIG##5##Figure 6A## indicate also the possibility that the transcriptional signal is turned on and then quickly turned off, so that the subsequent behaviour of the two genes is mainly determined by the degradation process alone. Accordingly, transcripts stabilities – as measured after transcriptional shut-off ##REF##16166263##[45]## – significantly differ in value.</p> <p>Cell Division Cycle 6 (CDC6) is a component of the pre-replicative complex essential for the initiation of DNA replication, normally expressed at the end of mitosis. It has a lead-lag relationships with ASH1 (##FIG##5##Figure 6B##) which encodes a GATA-like transcription factor localized at daughter cells where it serves to repress the late G1-specific transcription of HO and preventing mating-type switching ##UREF##12##[46]##. Our analysis suggests the presence of a common regulatory signal. In fact, both genes are key regulators of separate biological processes that are simultaneously activated by the SWI5 transcription factor ##REF##16522208##[15]##,##REF##7641697##[47]##,##REF##10409653##[48]##. Moreover, our analysis also suggests that the CDC6 transcript is fairly unstable. Consistently, the CDC6 protein is unstable ##REF##7641697##[47]##.</p> <p>SWI5 encodes a key transcription factor that activates transcription of genes expressed at the M/G1 boundary and in G1 phase of the cell cycle. NCE102 is a non-classical export protein involved in alternative clearance/detoxification pathway to eliminate damaged material ##UREF##13##[49]##. They display a lead-lag relationship (##FIG##5##Figure 6C##) and, in fact, they are both targets of the M-phase activating complex FKH2/NDD1 according to ChIP experiments ##REF##16522208##[15]## and large differences of their half-life values are observed after transcription inhibition ##REF##11972065##[26]##. The gene expression profiles shown in ##FIG##5##Figure 6C## reflect the prototypical situation of peak delay depicted in ##FIG##1##Figure 2C##.</p> <p>YOX1 is a transcription factor involved in the repression of ECB acitivity ##UREF##12##[46]## thus contributing to move the cycle forward. YOX1 shows a lead-lag relationship with MNN1 (##FIG##5##Figure 6D##), a gene encoding a cell wall glycoprotein ##REF##16672383##[50]##. Consistently, they are both targets of SBF according to ChIP experiments ##REF##16522208##[15]## and the transcripts have different half-life values ##REF##11972065##[26]##. Moreover, looking at the time profiles depicted in ##FIG##5##Figure 6D##, one may argue that during the second half of the cycle another transcription factor is active at the MNN1 promoter.</p> <p>All the above examples consist of pairs of genes that are under the control of the same transcription factor and that show differential mRNA stability values consistent with their lead-lag relationship (except for CDC6 transcript whose experimental half-life is not available). Moreover, it is worth noting that large differences in half-lives value (as in the cases shown in ##FIG##5##figure 6C and 6D##) significantly affect the overall time profiles producing also an evident peak shift.</p> <p>Finally, it is worth noting that the lead-lag relationship is symmetrical and, therefore, it does not provide information about which gene is “lead” and which is “lag”. However, such information can be easily obtained by visual inspection. In fact, from ##FIG##5##Figure 6A##, the lead gene is the one with the steepest decaying profile having, consistently, a smaller half-life. Moreover, from ##FIG##5##Figure 6C##, one can see that the gene with the larger half-life displays a delayed peak and, therefore, it corresponds to the lag gene.</p> <title>Dynamic formation of the replication complex</title> <p>Many studies have focused on the relationship between gene expression time courses and the formation of protein complexes. Interestingly, Jansen <italic>et al.</italic>
##REF##11779829##[31]## suggested to classify protein complexes as either <italic>permanent</italic> or <italic>transient</italic>, with permament ones being maintained through most cellular conditions. They also found that, generally, permanent complexes tend to have simultaneously correlated gene expression while transient ones do not. Moreover, they also noted that subunits of the same protein complex may show significant simultaneous expression. In particular, they studied gene expression of the replication complex in yeast and found a very low simultaneous correlation among subunits, not significantly different from a random control ##REF##11779829##[31]##. However, they also found two sub-complexes – the MCM complex and the DNA polymerases δ and ε complex – showing much greater simultaneous correlation.</p> <p>Using gene expression time profiles during one cell cycle (##REF##16912276##[43]##, dataset, alpha_38 time series) for the genes encoding MCM proteins (MCM cluster) and DNA polymerases and ε (POL cluster), we computed simultaneous and lead-lag R<sup>2</sup> and the scatterplots of the resulting values for gene pairs belonging to the two different sub-complexes are shown in ##FIG##6##Figure 7##, panel B. As a negative control we used a group of 5 simultaneously expressed genes (R<sup>2</sup>&gt;0.7) coding for proteins of the cytoplasmic ribosomal large subunit (RPL4A, RPL4B, RPL1A, RPL1B, RPP0 denoted by RIB cluster). Ribosomal proteins are under the transcriptional control of IFH1/FHL1 ##REF##15128447##[51]##,##REF##15616568##[52]## whereas the replication complex is regulated by the transcription factors MBF/SBF ##REF##15128447##[51]##. The scatterplot reported in ##FIG##6##Figure 7## of the simultaneous <italic>vs.</italic> lead-lag R<sup>2</sup> values shows that, whereas the POL/MCM pairs display high values of lead-lag R<sup>2</sup> and low values of simultaneous R<sup>2</sup>, the control pairs POL/RIB and MCM/RIB display a very different pattern spread over a larger range thus denoting the absence of any meaningful relationship.</p> <p>
##FIG##6##Figure 7## makes very clear that the gene expression of the MCM and the DNA polymerases δ and ε subcomplexes is significantly simultaneously correlated within the same group whereas such correlation dramatically drops if we consider pairs of genes belonging to different subcomplexes. In fact, the sample distribution of the simultaneous R<sup>2</sup> between the two clusters is spread over the range [0,0.5] thus showing the absence of any significant level of simultaneous correlation. By contrast, the between clusters lead-lag R<sup>2</sup> histogram is concentrated in the highest part of the range close to 1. The high values of the lead-lag R<sup>2</sup> strongly suggest that mRNA stability may play a fundamental role in the dynamic formation of multiple protein complexes. Accordingly, the average half-life measured after transcriptional shut-off of the MCM group is 14±6 min (the “lead” genes). and that of the POL group is 19±6 min ##REF##11972065##[26]## (the “lag” genes). The presence of lead-lag relationships between transient sub-complexes is briefly discussed in the supporting information file ##SUPPL##0##Text S1##.</p> <title>Conclusions</title> <p>The expression of genes in the cell is to a large extent controlled at the level of mRNA accumulation. One key point in the analysis of gene expression dynamics is that mRNA abundance is determined by two regulated processes: transcription and degradation both specifically affecting transcript levels. Computational analysis of genome-wide expression time series has shown that clusters of co-expressed (<italic>i.e.</italic> simultaneously correlated) profiles often provide clues for the presence of common transcription factors regulating both genes. Such computational analysis (known as “clustering”) is very useful since it allows the prediction of the underlying regulatory actions based exclusively on the available gene expression data obtained from a given experiment. The rationale behind such belief is a sort of a “guilty by association” approach: genes' products appearing and disappearing at the same time are likely to have some common transcriptional regulation. Nevertheless, it may well be the case that the same transcriptional signal regulating two (or more) genes may yield quite different outcomes on each transcript. In fact, a number of biological events following transcription may selectively affect cytoplasmic mRNA abundance, such as, for example, the activity of the enzymatic machinery involved in mRNA processing and degradation. In order to address this issue, we provided a novel computational methodology that, based exclusively on the available gene expression data, is able to effectively predict co-regulation even with variation in the dynamic response due to mRNA stability differences. Moreover, our approach also captures the relation of simultaneous or time shifted co-expression so that it provides a single integrative general index – the <italic>lead-lag R<sup>2</sup></italic>−able to uncover the presence of a common regulatory signal underlying gene expression time dynamics also at the post-transcriptional level.</p> <p>In order to test the validity of our approach on real data, we used yeast genome-wide cell-cycle expression time series obtained by several independent groups using different synchronization methods. In fact, by doing so, we could integrate the available cell cycle data and obtain a much more reliable aggregated dataset. We considered those gene pairs with the highest lead-lag R<sup>2</sup> values and found the prediction for the presence of a common transcription factor to be highly consistent with protein-DNA binding data (ChIP experiments). Our results clearly indicate that co-regulation is not generally equivalent to simultaneous expression.</p> <p>We believe that the same analysis can be successfully used to predict post-transcriptional regulation, <italic>i.e.</italic> the presence of a common mechanism able to stabilize or de-stabilize specific transcripts, as for the members of the PUF proteins family ##REF##15024427##[2]##. Moreover, we envisage the possibility that our methodology could be used on different data and organisms and thus providing a computational support to the understanding of transcriptional and post-transcriptional networks, given the recent growing interest in the post-transcriptional regulation layer ##REF##17245413##[1]## of gene expression (miRNA) and its role in many diseases, such as cancer. Finally, the characterization of the replication complex in terms of lead-lag relationships among gene expression time profiles of its sub-complexes suggests the possibility that our analysis could be effectively used as a tool for predicting the formation of transient multiple protein complexes.</p>
<title>Conclusions</title> <p>The expression of genes in the cell is to a large extent controlled at the level of mRNA accumulation. One key point in the analysis of gene expression dynamics is that mRNA abundance is determined by two regulated processes: transcription and degradation both specifically affecting transcript levels. Computational analysis of genome-wide expression time series has shown that clusters of co-expressed (<italic>i.e.</italic> simultaneously correlated) profiles often provide clues for the presence of common transcription factors regulating both genes. Such computational analysis (known as “clustering”) is very useful since it allows the prediction of the underlying regulatory actions based exclusively on the available gene expression data obtained from a given experiment. The rationale behind such belief is a sort of a “guilty by association” approach: genes' products appearing and disappearing at the same time are likely to have some common transcriptional regulation. Nevertheless, it may well be the case that the same transcriptional signal regulating two (or more) genes may yield quite different outcomes on each transcript. In fact, a number of biological events following transcription may selectively affect cytoplasmic mRNA abundance, such as, for example, the activity of the enzymatic machinery involved in mRNA processing and degradation. In order to address this issue, we provided a novel computational methodology that, based exclusively on the available gene expression data, is able to effectively predict co-regulation even with variation in the dynamic response due to mRNA stability differences. Moreover, our approach also captures the relation of simultaneous or time shifted co-expression so that it provides a single integrative general index – the <italic>lead-lag R<sup>2</sup></italic>−able to uncover the presence of a common regulatory signal underlying gene expression time dynamics also at the post-transcriptional level.</p> <p>In order to test the validity of our approach on real data, we used yeast genome-wide cell-cycle expression time series obtained by several independent groups using different synchronization methods. In fact, by doing so, we could integrate the available cell cycle data and obtain a much more reliable aggregated dataset. We considered those gene pairs with the highest lead-lag R<sup>2</sup> values and found the prediction for the presence of a common transcription factor to be highly consistent with protein-DNA binding data (ChIP experiments). Our results clearly indicate that co-regulation is not generally equivalent to simultaneous expression.</p> <p>We believe that the same analysis can be successfully used to predict post-transcriptional regulation, <italic>i.e.</italic> the presence of a common mechanism able to stabilize or de-stabilize specific transcripts, as for the members of the PUF proteins family ##REF##15024427##[2]##. Moreover, we envisage the possibility that our methodology could be used on different data and organisms and thus providing a computational support to the understanding of transcriptional and post-transcriptional networks, given the recent growing interest in the post-transcriptional regulation layer ##REF##17245413##[1]## of gene expression (miRNA) and its role in many diseases, such as cancer. Finally, the characterization of the replication complex in terms of lead-lag relationships among gene expression time profiles of its sub-complexes suggests the possibility that our analysis could be effectively used as a tool for predicting the formation of transient multiple protein complexes.</p>
<p>Analyzed the data: LF ADS SS GM IR. Wrote the paper: LF ADS GM IR. Conceived methodology: LF ADS. Contributed to the development and validation of the method: GM IR.</p> <p>Current methods for the identification of putatively co-regulated genes directly from gene expression time profiles are based on the similarity of the time profile. Such association metrics, despite their central role in gene network inference and machine learning, have largely ignored the impact of dynamics or variation in mRNA stability. Here we introduce a simple, but powerful, new similarity metric called <italic>lead-lag R<sup>2</sup></italic> that successfully accounts for the properties of gene dynamics, including varying mRNA degradation and delays. Using yeast cell-cycle time-series gene expression data, we demonstrate that the predictive power of lead-lag R<sup>2</sup> for the identification of co-regulated genes is significantly higher than that of standard similarity measures, thus allowing the selection of a large number of entirely new putatively co-regulated genes. Furthermore, the lead-lag metric can also be used to uncover the relationship between gene expression time-series and the dynamics of formation of multiple protein complexes. Remarkably, we found a high lead-lag R<sup>2</sup> value among genes coding for a transient complex.</p> <title>Author Summary</title> <p>Microarrays provide snapshots of the transcriptional state of the cell at some point in time. Multiple snapshots can be taken sequentially in time, thus providing insight into the dynamics of change. Since genome-wide expression data report on the abundance of mRNA, not on the underlying activity of genes, we developed a novel method to relate the expression pattern of genes, detected in a time-series experiment, using a similarity measure that incorporates mRNA decay and called <italic>lead-lag R<sup>2</sup></italic>. We used the lead-lag R<sup>2</sup> similarity measure to predict the presence of common transcription factors between gene pairs using an integrated dataset consisting of 13 yeast cell-cycles. The method was benchmarked against six well-established similarity measures and obtained the best true positive rate result, around 95%. We believe that the lead-lag analysis can be successfully used also to predict the presence of a common mechanism able to modulate the degradation rate of specific transcripts. Finally, we envisage the possibility to extend our analysis to different experimental conditions and organisms, thus providing a simple off-the-shelf computational tool to support the understanding of the transcriptional and post-transcriptional regulation layer and its role in many diseases, such as cancer.</p>
<title>Supporting Information</title>
<p>The authors thank Tim Gardner, Arun Krishnan, Feng He, and Alessandro Giuliani for their critical reading of the manuscript and constructive suggestions.</p>
<fig id="pcbi-1000141-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g001</object-id><label>Figure 1</label><caption><title>Combination of transcriptional and post-transcriptional regulation.</title><p>Gene expression time profiles obtained by simulations using various types of regulation of the transcription rate and degradation rate. Panels A and B refer to the situation in which two genes have the same constant degradation rate but different transcription rate signals (persistent or transient). Panels C and D refer to the situation in which two genes have the same constant transcription rate but different degradation rate signals (persistent or transient). Normalized time profiles are linearly scaled such that their values remain bounded between 1 and 2, <italic>i.e.</italic> by setting to 1 the lowest value and to 2 the highest value so that the peak-to-peak amplitude is set to 1.</p></caption></fig> <fig id="pcbi-1000141-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g002</object-id><label>Figure 2</label><caption><title>Relationships between gene expression time profiles.</title><p>Typical behaviors of genes related by a simultaneous (Panel A), time-delayed (Panel B) and peak-delayed (Panel C) expression pattern. The time delayed profiles may be the result of different switching times in transcriptional activation whereas the peak delayed profile may be the result of the same transcriptional signal and different mRNA stabilities. Normalization is performed by setting to 1 the lowest value and to 2 the highest value so that the peak-to-peak amplitude is set to 1.</p></caption></fig> <fig id="pcbi-1000141-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g003</object-id><label>Figure 3</label><caption><title>Predicting co-regulation: ROC curves.</title><p>The ROC curves for co-regulation prediction were computed using 10 transcription factors involved in the cell cycle and assuming as true targets the DNA binding data provided by MacIsaac <italic>et al.</italic>
##REF##16522208##[15]##, with <italic>p</italic>-value for binding 0.001. Two varying thresholds have been used for constructing the ROC curves; the higher percentile <italic>p</italic> (ranging from 50% to 90%) for the prediction of co-regulation, and the symmetric (100−<italic>p</italic>) percentile (ranging from 50% to 10%) for no co-regulation. Panel A shows ROC curves corresponding to R<sup>2</sup> values averaged over all the 7 available datasets. Numbers below marks (circles and crosses) indicate the percentile of the distribution of R<sup>2</sup> values used for selecting the lower and upper thresholds. Panel B shows ROC curve for each dataset obtained using as a threshold for co-regulation prediction only the 90<sup>th</sup> percentile of the corresponding distribution.</p></caption></fig> <fig id="pcbi-1000141-g004" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g004</object-id><label>Figure 4</label><caption><title>Sample gene pairs distribution of simultaneous R<sup>2</sup> and lead-lag R<sup>2</sup> minus simultaneous R<sup>2</sup> versus half-life ratios.</title><p>Panel A shows the bar plot of the mean values of the R<sup>2</sup> distributions for co-regulated pairs. On the top of each bar is indicated the standard deviation as a percentage of the mean value. The difference between the means of the first and the last bin of the lead-lag R<sup>2</sup> minus simultaneous R<sup>2</sup> is significant according to a <italic>t</italic>-test with a confidence level of 95%, (<italic>p</italic>-value 10<sup>−9</sup>) and the difference between the means of the first and the last bin of the simultaneous R<sup>2</sup> is mildly significant according to a <italic>t</italic>-test with a confidence level of 95%, (<italic>p</italic>-value = 0.03). According to the <italic>t</italic> test, no other difference is significant. The number of gene pairs contained in each bin are: 9128, 3726, 1707, 919 and 1260, respectively. Blue bars indicate the distribution of simultaneous R<sup>2</sup> whereas green bars indicate the distribution of lead-lag R<sup>2</sup> minus simultaneous R<sup>2</sup>. Panel B shows the scatterplot of the <italic>Z</italic>-scores corresponding to 100000 randomizations of the half-life measurements. Whereas the <italic>Z</italic> score of the simultaneous R<sup>2</sup> does not show any significant change after half-life ratios randomization for each bin, the <italic>Z</italic> score of the lead-lag R<sup>2</sup> minus the simultaneous R<sup>2</sup>, does show a significant change in the first and the last bin.</p></caption></fig> <fig id="pcbi-1000141-g005" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g005</object-id><label>Figure 5</label><caption><title>Performance comparisons of similarity measures for co-regulation prediction.</title><p>The similarity measures considered for comparisons are: the square of Pearson, Spearman's rank, Kendall's tau, time-warped, time-delayed and lead-lag correlation. Sensitivity (Panel A) and specificity (Panel B) are shown as a function of the chosen percentile threshold. The results obtained using the for cosine correlation are not reported since it performed as a random choice.</p></caption></fig> <fig id="pcbi-1000141-g006" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g006</object-id><label>Figure 6</label><caption><title>Examples of lead-lag relationships of key cell cycle regulators.</title><p>The gene pairs reported in the figure are: CLB6/GIN4 (Panel A), CDC6/ASH1 (Panel B), SWI5/NCE102 (Panel C) and YOX1/MNN1 (Panel D). Gene expression time profiles are taken from Pramila <italic>et al.</italic>
##REF##16912276##[43]##, alpha_38 dataset and the expression values are normalized with respect to peak-to-peak amplitude. Each half-life dataset is indicated in brackets: “W” denotes the Wang dataset ##REF##11972065##[26]##, “K” denotes the Kuai dataset ##REF##16166263##[45]##.</p></caption></fig> <fig id="pcbi-1000141-g007" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000141.g007</object-id><label>Figure 7</label><caption><title>Lead-lag analysis of gene expression profiles of two components of the replication complex.</title><p>The MCM cluster is composed by the time course of genes belonging to the MCM protein complex and the POL cluster by those belonging to the DNA polymerases δ and ε complexes. Panel A shows gene expression time profiles during one cell cycle (##REF##16912276##[43]##, alpha38 dataset) of the two groups of genes. In Panel B the scatterplots are computed for simultaneous and lead-lag R<sup>2</sup> by considering within pairs for each cluster (MCM and POL) and between pairs (MCM/POL). The between pairs MCM/POL display a low value of simultaneous R<sup>2</sup> and an high value of the lead-lag R<sup>2</sup> consistently with the formation of the MCM/POL complex. Panel B also shows the results on a negative control composed by a group of simultaneously expressed genes encoding proteins that belong to the cytoplasmic ribosomal large subunit (RIB cluster). The between pairs MCM/RIB and POL/RIB display low values of both simultaneous and lead-lag R<sup>2</sup> consistently with the absence of any interaction or co-regulation between the ribosomal proteins and the replication complex.</p></caption></fig>
<disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula><label>(1)</label></disp-formula> <disp-formula></disp-formula> <disp-formula><label>(2)</label></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula><label>(3)</label></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula> <disp-formula></disp-formula>
<supplementary-material content-type="local-data" id="pcbi.1000141.s001"><label>Text S1</label><caption><p>Supporting Information file</p><p>(0.10 MB DOC)</p></caption></supplementary-material>
<fn-group><fn fn-type="COI-statement"><p>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p>This work was partially supported by a grant from ASI, Biotechnology Program.</p></fn></fn-group>
<graphic xlink:href="pcbi.1000141.g001"/> <graphic xlink:href="pcbi.1000141.g002"/> <graphic xlink:href="pcbi.1000141.e001.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.g003"/> <graphic xlink:href="pcbi.1000141.g004"/> <graphic xlink:href="pcbi.1000141.g005"/> <graphic xlink:href="pcbi.1000141.g006"/> <graphic xlink:href="pcbi.1000141.g007"/> <graphic xlink:href="pcbi.1000141.e002.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e003.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e004.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e005.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e006.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e007.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e008.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e009.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e010.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e011.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e012.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e013.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e014.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e015.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e016.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e017.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000141.e018.jpg" mimetype="image" position="float"/>
<media xlink:href="pcbi.1000141.s001.doc"><caption><p>Click here for additional data file.</p></caption></media>
[{"label": ["3"], "element-citation": ["\n"], "surname": ["Guan", "Zheng", "Tang", "Liu", "Zinkel"], "given-names": ["Q", "W", "S", "X", "RA"], "year": ["2006"], "article-title": ["Impact of nonsense-mediated mRNA decay of the global expression profile of budding yeast."], "source": ["Plos Genetics"], "volume": ["2"], "fpage": ["1924"], "lpage": ["1943"], "comment": ["doi:10.1371/journal.pgen.0020203"]}, {"label": ["5"], "element-citation": ["\n"], "surname": ["Tavazoie", "Hughes", "Campbell", "Cho", "Church"], "given-names": ["S", "JD", "MJ", "RJ", "GM"], "year": ["1999"], "article-title": ["Systematic determination of genetic network architecture."], "source": ["Nature Gen"], "volume": ["22"], "fpage": ["281"], "lpage": ["285"]}, {"label": ["11"], "element-citation": ["\n"], "surname": ["Brazhnik", "de la Fuente", "Mendes"], "given-names": ["P", "A", "P"], "year": ["2002"], "article-title": ["Gene networks: how to put the function in genomics."], "source": ["Trends Biotech"], "volume": ["20"], "fpage": ["467"], "lpage": ["472"]}, {"label": ["13"], "element-citation": ["\n"], "surname": ["Basso", "Margolin", "Stolovitzky", "Klein", "Dalla-Favera"], "given-names": ["K", "AA", "G", "U", "R"], "year": ["2005"], "article-title": ["Reverse engineering of regulatory networks in human B cells."], "source": ["Nature Gen"], "volume": ["37"], "fpage": ["382"], "lpage": ["390"]}, {"label": ["14"], "element-citation": ["\n"], "surname": ["Wit", "McClure"], "given-names": ["E", "J"], "year": ["2004"], "source": ["Statistics for Microarrays Design, Analysis and Inference"], "publisher-loc": ["Chichester, UK"], "publisher-name": ["John Wiley & Sons Ltd"]}, {"label": ["18"], "element-citation": ["\n"], "surname": ["Pilpel", "Sudarsanam1", "Church"], "given-names": ["Y", "P", "GM"], "year": ["2001"], "article-title": ["Identifying regulatory networks by combinatorial analysis of promoter elements."], "source": ["Nature Gen"], "volume": ["29"], "fpage": ["153"], "lpage": ["159"]}, {"label": ["20"], "element-citation": ["\n"], "surname": ["Yu", "Luscombe", "Qian", "Gerstein"], "given-names": ["H", "NM", "J", "M"], "year": ["2003"], "article-title": ["Genomic analysis of gene expression relationships in transcriptional regulatory networks."], "source": ["Trends Gen"], "volume": ["19"], "fpage": ["422"], "lpage": ["427"]}, {"label": ["29"], "element-citation": ["\n"], "surname": ["Franklin", "Powell", "Abbas"], "given-names": ["G", "JD", "Emami-Naeini"], "year": ["2002"], "article-title": ["Feedback Control of Dynamic Systems."], "source": ["Prentice-Hall, 4"], "sup": ["th"]}, {"label": ["30"], "element-citation": ["\n"], "surname": ["Ge", "Liu", "Church", "Vidal"], "given-names": ["H", "Z", "GM", "M"], "year": ["2001"], "article-title": ["Correlation between transcriptome and interactome mapping data from "], "italic": ["Saccharomyces cerevisiae"], "source": ["Nature Gen"], "volume": ["29"], "fpage": ["482"], "lpage": ["486"]}, {"label": ["36"], "element-citation": ["\n"], "surname": ["Arkin", "Shen", "Ross"], "given-names": ["A", "PD", "J"], "year": ["1997"], "article-title": ["A Test Case of Correlation Metric Construction of a Reaction Pathway from Measurements."], "source": ["Science"], "volume": ["277"], "fpage": ["1275"], "lpage": ["1279"]}, {"label": ["42"], "element-citation": ["\n"], "surname": ["Sankoff", "Kruskal"], "given-names": ["D", "J"], "year": ["1983"], "source": ["Time warps, string edits, and macromolecules: the theory and practice of sequence comparison"], "publisher-loc": ["Reading, MA"], "publisher-name": ["Addison Wesley"]}, {"label": ["44"], "element-citation": ["\n"], "surname": ["Breeden"], "given-names": ["LL"], "year": ["2000"], "article-title": ["Cyclin transcription: timing is everything."], "source": ["Curr Biol"], "volume": ["10"], "fpage": ["586"], "lpage": ["588"]}, {"label": ["46"], "element-citation": ["\n"], "surname": ["Breeden"], "given-names": ["LL"], "year": ["2003"], "article-title": ["Periodic transcription: a cycle within a cycle."], "source": ["Curr Biol"], "volume": ["13"], "fpage": ["31"], "lpage": ["38"]}, {"label": ["49"], "element-citation": ["\n"], "surname": ["Desmyter", "Verstraelen", "Dewaele", "Libert", "Contreras"], "given-names": ["L", "J", "S", "C", "R"], "year": ["2007"], "article-title": ["Nonclassical export pathway: overexpression of NCE102 reduces protein and DNA damage and prolongs lifespan in an SGS1 deficient "], "italic": ["Saccharomyces cerevisiae"], "source": ["Biogerontol"], "volume": ["8"], "fpage": ["527"], "lpage": ["535"]}]
PLoS Comput Biol. 2008 Aug 1; 4(8):e1000141