The Gene Expression Prediction Challenge of DREAM3

The challenge for predicting gene expression provided by the DREAM project is of great importance to explore the benefits and bottlenecks of the state-of-the-art algorithms in a fair competition. It represents a solution to the non-trivial problem of designing relevant challenges which at same time addresses biological and computational interesting problems. From the DREAM web-site [wiki.c2b2.columbia.edu/dream/index.php/The_DREAM_Project, accessed October 10, 2008] we quote for the gene expression prediction challenge within DREAM3:

Gene expression time course data is provided for four different strains of yeast (S. Cerevisiae), after perturbation of the cells. The challenge is to predict the rank order of induction/repression of a small subset of genes (the “prediction targets” in one of the four strains, given complete data for three of the strains, and data for all genes except the prediction targets in the other strain. Predictors are also allowed to use any information that is in the public domain but are expected to be forthcoming about what information was used.

Background. GAT1, GCN4, and LEU3 are yeast transcription factors. Each of these transcription factors has something to do with controlling genes involved in nitrogen or amino acid metabolism. The genes are not essential because strains that have perfect deletions of any of these genes are viable. In this challenge, we provide gene expression data from four strains: (i) a strain that is wild-type for all three transcription factors (wt, or parental), (ii) a strain that is identical to the parental strain except that it has a deletion of the GAT1 gene (gat1Δ), (iii) a strain that is identical to the parental strain except that it has a deletion of the GCN4 gene (gcn4Δ), and (iv) a strain that is identical to the parental strain except that it has a deletion of the LEU3 gene (leu3Δ).

Expression levels were assayed separately in all four strains following the addition of 3-aminotriazole (3AT). 3AT is an inhibitor of an enzyme in the histidine biosynthesis pathway and, in the appropriate media (which is the case in these experiments) inhibition of the histidine biosynthetic pathway has the effect of starving the cells for this essential amino acid.

Data from eight time points was obtained from 0 to 120 minutes. Time t = 0 means the absence of 3AT.

The challenge. Predict, for a set of 50 genes, the expression levels in the gat1Δ strain in the absence of 3-aminotriazole (t = 0) and at 7 time points (t = 10, 20, 30, 45, 60, 90 and 120 minutes) following the addition of 3AT. Absolute expression levels are not required or desired; instead, the fifty genes should be ranked according to relative induction or repression relative to the expression levels observed in the wild-type parental strain in the absence of 3AT.

This challenge is biologically relevant, and the fact a gold standard exists but is hidden makes the challenge objective and fair. Further, the probe names were given, which allows for data integration of publicly available experiments and a priori knowledge, making the challenge even more realistic in describing a situation which can occur in one's laboratory. However, the problem is somewhat different from the normal setting in systems biology where the aim is not only to predict future experiments but also to obtain interpretable models from which we can gain an increased biological understanding [8]–[10]. The data for this DREAM challenge was kindly delivered by Neil Clarke and co-workers, a fact which was revealed first after the submission period for predictions had closed. We will henceforth refer to this data as the “DREAM data”.

Modelling Assumptions

The quest for modelling gene networks has taken many different forms during the last decade [5], [9], [10], [13]. One of the major modelling frameworks is provided by ODEs, ordinary differential equations, which is the form we exploit here for the development of our algorithm.

An often utilized approach for large-scale modelling of gene regulatory networks is to only consider the transcripts, and thereby letting all interactions be projected onto the space of genes only [14]. By this, one obtains a gene-to-gene network, sometimes referred to as “‘influential’ gene regulatory network” [5]. Further, the amount of data available makes most models except linear ones behind reach, and one therefore assumes a linear relation between the expression rates and the expression levels [15]–[17] (or some non-linear transformation of them, the equations still being linear, though [18]). The rationale for this assumption of linearity is normally expressed as a belief in the system being close to some equilibrium or working point, and thus an expansion keeping only the linear terms would be appropriate. The basic dynamical equation then looks like(1)Here denotes the expression level of gene at time , and the dot on the left hand side denotes its time-derivative. The coefficients and the elements are to be inferred and describe the gene regulatory system on a gene-to-gene level. The interactions, i.e., the non-zero , can both correspond to (semi)-direct interaction, for instance when denotes the RNA-level of a transcription factor which binds to the promoter of gene , and indirect interactions, for instance when gene stands for an inhibitor of a transcription factor which in its turn regulate gene . Note that the coefficient can be both positive, describing activation, and negative, describing inhibition. Since this is a dynamical equation, it can be used for studying the time-development of the model, including aspects as stability and flexibility [19].

Here, however, we are not primarily interested in the dynamics or in the derivatives of the 50 genes in the gat1Δ strain which were removed from the data file, but in the prediction of their expression levels. By denoting this set of 50 genes as , we reformulate and generalize Eq. (1) to(2)for all .

Model Selection

In order to obtain a ranking list based on the expressions of the 50 genes in , we must predict the values of for all from (2). For this, we need explicit values of and . Given the values of and for all and for all , this can be formulated as a minimization problem with the objective function(3)We will here stick to the cases where , least absolute deviations (LAD), and , least squares (LS), respectively. The penalty term is utilized for discriminating among models, and is to be determined later. For the moment, it is set to zero. The experiments are assumed to have been performed at time , where we for ease of notation let the index run over all three complete time-series (the gat1Δ-series, i.e., the series which should be predicted, is never utilized for the inference, since it is from this series we determine the output which is submitted for the challenge – hence there is a risk of over fitting if we utilize it twice).

The DREAM data are measured by Affymetrix chips of 9335 probes, and obtained from two biological and two technical replicates. We map the probe-names onto unique gene names, which leave 7804 units, where we use mean values when more than one probe corresponded to one gene. Furthermore, we approximate the derivatives from central differences, except at the end points of each series. At time we set the derivative to zero, since it was right before the addition of 3AT, and at the final time we approximate the derivatives from backward differences.

Throughout the article, we utilize cross-validation (CV) to discriminate among models. We hold one of the three time-series provided by DREAM out from the inference, and utilize the other two, and occasionally also other data sets, for finding the searched parameters. We then use data from the left out strain to predict the expression values of the 50 searched genes for each time point in the this series, rank them according to the predicted levels such that the highest expressed gene obtain rank number one, second highest rank number two, etc., and calculate the Spearman rank correlation with the observed ranking of the same series. This is repeated three times, holding each of the provided time-series out a time. We end up with 24 different ranking lists for the 50 genes in , eight for each time-series. Henceforth, we will refer to this procedure simply as “cross-validation”. Note that this approach to only hold one of the given time-series out still holds also when external datasets are introduced.

Before we start exploring various versions of the penalty term in (3), we try to simplify the model (2). The strategy is to primarily work with the DREAM data, in order to reduce the model. When this first reduction is obtained, we will utilize also other publicly available data in order to further strengthen the predictive power of our mathematical model. This first model selection is performed among the models with perfect fits, i.e., the ones where the terms for the first sum of absolute values in (3) all are zero, making the exponent irrelevant, and we get an idea on what to include. This perfect fit is possible since for each gene value we should predict, that is, for each i, we have 15 509 free parameters to determine from 16 linear equations (one per time point). We explore the following three scenarios:

• Only expression values are included
• Both expression values and expression rates (derivatives) are included
• Only expression rates (derivatives) are included

By picking the solution with zero value for the objective function (without penalty term) and choosing the coefficients and such that their L1-norm and L2-norm , respectively, are minimized, we obtain the values of the Spearman rank correlations presented in Table 1. The rank correlations are for each series comprising eight time points as well as for the full 24 lists simultaneously (referred to as “overall”) obtained by cross-validation. Note that this is not the same as the mean value of each time-series or time point.

We see that the highest values for the correlations are obtained when we only include the expression levels. Inclusion of the expression rates makes the result slightly worse, except for the least squares where the correlations are equal. However, with the same predictive power, we apply Occam's razor and prefer the simplest model. To only use the rates gives the least satisfying result of them all. Therefore, in the sequel, we choose to discard all derivative terms and determine the parameters according to(4)Here we use the short-hand notation for and drop the intercept term since we also centre the data. What remains is to determine the value of m and the explicit form and value of the penalty term.

By choosing the exponent we turn the problem into Least Squares (LS), while gives Least Absolute Deviations (LAD). Least squares are known to be sensitive for outliers, while LAD is not affected by such at all. However, the solutions obtained from LAD are unstable, in the sense that small variations in data might result in large variations of the inferred parameters, that is, least squares are a more stable inference method. Indeed, LAD can be compared with picking the median, which may be devastating when we are already short of data [20]. Both methods will here be explored.

The penalty term can take many different forms. A review for least squares of more classical forms as Mallow's Cp, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Minimum Description Length (MDL) etc in the context of gene networks can be found in [21]. We will here stick to forms which are computationally attractive, which is important since the size of the present problem makes methods based on exhaustive search practically impossible. We choose for least squares the form(5)This is known as the elastic net [22] and is a convex combination of the well established methods of ridge regression [23] and of the lasso , [24]. For LAD, least absolute deviations, only the “lasso-form”, i.e., , is of practical interest due to implementation considerations. This case is called regularized least absolute deviation (RLAD) [25].

This convex combination is a compromise between the two goals of a good solution, to have good predictive power and to be interpretable [5]. The L1 penalty in the lasso favours sparse solutions, i.e., it performs effectively a subset selection, but is greedy in the case of correlated predictors since it only picks the most correlated predictor. On the other hand, the L2 penalty in ridge regression keeps all predictors as non-zero with probability one, but is not greedy as the lasso. It has been argued that this compromise can overcome some of the limitations of the lasso and ridge regression, with the preservation of the benefits from each of the pure methods [22]. In Table 2 we investigate what kind of penalty is most predictive for our purpose. We investigate the elastic net for least squares and RLAD for least absolute deviations on the DREAM data set. The parameters and are obtained from an exhaustive grid search over the parameter space and 1000 equidistant values of between zero and the upper limit (defined as the limit where the found solution coincides with the solution without penalty term). The parameters are eventually chosen with cross-validation. In practice, we utilize the R-package glmnet by Friedman et al. [26]. It turns out that the values of vary with in the whole interval between zero and one, with a median of 0.75. That is, the solution is “lasso-like”, but still the maximum in-degree turns out to be 488 (with median 23.5) which is not possible for a pure lasso where the number of predictors cannot exceed the number of experiments.

We show in Table 2 the results of such an optimization. The entries are the Spearman rank order correlations obtained from cross-validation, i.e., they are the best measure of performance we have without adding any prior information. The general conclusion from Table 2 is that least squares regression combined with the elastic net is to prefer among the methods. From a correlation of 0.794 it is worthy to proceed with further improvements by integrating more data into the inference process. Worth noting is that in this process, the parameters and are not fixed, but will be recalculated by new cross-validation procedures.

Data Integration

One way to improve the performance of the algorithm is to include more data. This is a challenging problem which is crucial for all kinds of large-scale inference problems [5]. Previous work has shown, however, that uncritical inclusion of more data can even decrease the performance of an algorithm [27]. For yeast, huge amounts of public data are available, which potentially can improve the quality of the inferred models. The main issue is to select and process only the relevant data sets. This must be done in a careful, supervised way to avoid over fitting, because validation data are on the other hand severely limited. The here presented framework introduces two main possibilities, integration of more expression data and introduction of prior knowledge of regulatory interactions.

Prior knowledge of TF-DNA regulations.

The TF-DNA binding data are taken from Yeastract [32], [33], and give clues on which genes are likely to act as regulators. To take this knowledge into account, we utilize a modified form of the penalty in (5). Now, we choose a different penalty for each predictor such that we increase the probability for a non-zero if gene j somehow, to be specified below, is co-regulated with gene i. By utilizing the cross-validation scheme, we integrate this prior knowledge in a soft way, i.e., we bias the search to known relations, but allow also the possibility of novel links. Explicitly, expression (5) is modified by the introduction of parameters , called priors, reflecting the likeliness of gene not to be co-regulated with gene i, to look like(7)That is, the prior corresponds to our prior belief that there is no correspondence between gene j and gene i, since a high value of implies a high penalty, and thus there has to be a high correlation between the genes to include gene j as a predictor for gene i. We let the parameters be between zero and unity, with no loss of generality, since it is only the relative difference which is of relevance, while takes care of the global magnitude of the penalty.

The reason why we focus on co-regulation rather than regulatory interactions is that the values of the inference are based on transcript levels, and TFs are known to be expressed on a low level. Also, their activity is often determined by phosphorylation and other effects rather than their amount [34], which makes the levels even more dubious. Therefore, we concentrate not on the TFs themselves but on genes which are controlled by the same TFs and thus might act as predictors for each others. The rationale can also be formulated such that if two genes are regulated by the same set of TFs they are likely to have large (anti-) correlation, thus the level for one of the genes should be a good predictor for the other, see Figure 3. The explicit calculation is carried out through the TF interaction matrices obtained from Yeastract [32], [33]. We set if there is documented experimental evidence for TF-j binding upstream of gene i and if there is either experimental evidence or putative evidence (or both) for such a binding. Otherwise, the elements are set to zero. From these TF interaction matrices, we calculate the weighted shared fraction of TFs between gene j and gene i, , as(8)Note that is a symmetric matrix, , just as it should from its definition of a weighted shared fraction of TFs. The prior finally is chosen as:(9)where is a free global parameter to be determined later by cross-validation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Schematic view of how to obtain the weighed shared fraction of TFs.

We utilize as prior information that genes which are co-regulated are likely to be effective predictors of each other; the more TFs in common, the more likely to be co-regulated.

https://doi.org/10.1371/journal.pone.0009134.g003

To summarize the discussion above, the objective function takes the form(10)where the penalty term is(11)Thus, the data integration leaves us with 207,103 parameters () in total to determine. The global parameters result in the fractional value of 4142.06 parameters for each gene to predict. That is, the cross-validation procedure we utilize has to minimize expression (10) with a penalty term of the form presented in (11) and detailed in (8). We minimize it by a greedy stepwise procedure, again using the R-package glmnet [26].

Effectively, for each target gene , we start from the value of found without any prior and explore and individually. For each sets of values for these parameters, we run glmnet with 1000 values of (described above) and perform local searches in around its present value, changing it if the Spearman rank correlation increases. The parameters and are introduced in the order they increase the performance of the algorithm, but when a value is determined, we do not change it again. It turns out they are introduced in the order and , with the values and , respectively. They are obtained by in an iteratively refinement process, where the initial searches for are in the range between zero and one, and values for from the set . The refined search utilizes steps of the size for and 0.01 for . The result is shown in the second line of Table 3. We can see how this final step further increases the quality of the performance of the algorithm, although the increase this time is not as drastic as before. Worth noting is also how the Rosetta data contribute ten times more than the data from the ncbi omnibus. The rather low value of is also an indication that the positive contribution of known TF-DNA bindings is limited.

Finally, from a computational point of view, we remark that all implementations and calculations have been performed on an ordinary laptop in the languages R and Matlab. That is, the complexity of the problem is not worse than it can be handled in any laboratory.