 Methodology article
 Open Access
 Published:
Recursive regularization for inferring gene networks from timecourse gene expression profiles
BMC Systems Biology volume 3, Article number: 41 (2009)
Abstract
Background
Inferring gene networks from timecourse microarray experiments with vector autoregressive (VAR) model is the process of identifying functional associations between genes through multivariate time series. This problem can be cast as a variable selection problem in Statistics. One of the promising methods for variable selection is the elastic net proposed by Zou and Hastie (2005). However, VAR modeling with the elastic net succeeds in increasing the number of true positives while it also results in increasing the number of false positives.
Results
By incorporating relative importance of the VAR coefficients into the elastic net, we propose a new class of regularization, called recursive elastic net, to increase the capability of the elastic net and estimate gene networks based on the VAR model. The recursive elastic net can reduce the number of false positives gradually by updating the importance. Numerical simulations and comparisons demonstrate that the proposed method succeeds in reducing the number of false positives drastically while keeping the high number of true positives in the network inference and achieves two or more times higher true discovery rate (the proportion of true positives among the selected edges) than the competing methods even when the number of time points is small. We also compared our method with various reverseengineering algorithms on experimental data of MCF7 breast cancer cells stimulated with two ErbB ligands, EGF and HRG.
Conclusion
The recursive elastic net is a powerful tool for inferring gene networks from timecourse gene expression profiles.
Background
The inference of gene networks from timecourse microarray data can be defined as the process of identifying functional interactions between genes over time. Typically, a gene network is represented by a directed or undirected graph where nodes indicate genes encoded in a given organism of interest and edges represent various functional properties. The elucidation of gene networks has been expected to having an essential role for better understanding of molecular mechanisms and can be useful in the identification of new drug targets [1–5].
In this article, we use vector autoregressive (VAR) model [6, 7] to estimate gene networks from timecourse microarray data. The process of inferring gene networks based on the VAR model is to choose nonzero coefficients in the coefficient matrix, which can be considered as a problem of statistical model selection, especially as a variable selection problem [8]. Although a variety of variable selection methods have been developed, e.g., bestsubset selection [9], subset selection [9] and the lasso [10], these methods often suffer from the following crucial problems due to the limited number of samples (time points) compared with the large number of variables (genes) in timecourse microarray data.

1.
High computational cost for model selection: When the number of variables is m, there are m × 2^{m}candidate models in model selection. The bestsubset selection is computationally prohibitive when the number of variables is large.

2.
High correlation between variables: When the number of variables is much larger than the number of samples, two or more variables tend to be highly correlated [11]. In this situation, the coefficient estimates of the subset selection or the lasso may change erratically in response to small changes in the observed data, and thus the resulting models have poor performances [12, 13]. What is worse that these methods tend to select only one variable from the highly correlated variables [13] which can lead to reducing the number of true positives in gene network inference.
One solution for the above problems is to use a regularization method called elastic net [13] which minimizes a penalized loss function with l_{1} and l_{2}penalties of the coefficients. Applying an l_{1}penalty regularizes the least squares fit and shrinks some coefficients exactly to zero, i.e., achieves automatic variable selection, as the lasso does. Adding of an l_{2}penalty to an l_{1}penalty encourages a grouping effect so that highly correlated variables will be in or out of the model together. The elastic net is also capable of selecting a set of relevant variables with low computational effort even when the number of variables is much larger than the number of observations with LARSEN algorithm [13]. However, although VAR modeling with the elastic net succeeds in increasing the number of true positives, it also results in increasing the number of false positives. This is because the elastic net shrinks the same amount of the l_{1}penalty to each coefficient without assessing their relative importance.
To increase the capability of the elastic net, i.e., to decrease the number of false positives while keeping the high number of true positives in inferring gene networks based on VAR model, we propose a new class of regularization, called recursive elastic net, by incorporating relative importance of the coefficients. The recursive elastic net is a kind of iterative procedures that can reduce the number of false positives gradually by updating the relative importance of each VAR coefficient. The performance of variable selection strongly depends on the regularization parameters of the l_{1} and l_{2}penalty terms in the recursive elastic net. For selecting these regularization parameters automatically, we derive a modified Bayesian information criterion and a modified biascorrected Akaike information criterion.
Methods
Vector Autoregressive Model
We now consider gene expression data in timecourse microarray experiments at n time steps t ∈ {1, 2, ..., n}. Let y_{ t }= (y_{t,1}, ..., y_{t, m})' be the vector of m gene expressions at time step t. We assume that y_{ t }at time step t is generated from y_{t1}at the previous time step t  1 with a firstorder VAR model:
where c is a constant vector, B= [B_{ ij }]_{1≤i, j≤m}is an m × m autoregressive coefficient matrix, ϵ_{ t }is assumed to follow an mdimensional multivariate normal distribution N(0, Σ), and the notation denotes the transpose of the vector y_{ t }. Throughout the paper, we assume that Σ is a diagonal matrix, that is, . Note that the coefficient B_{ ij }measures the influence that node i exerts on node j and the nonzero B_{ ij }provides a functional connectivity which is related to Granger causality [14]. The Granger causality is a concept widely used in econometrics to analyze the relationships between variables and is based on the intuitive idea that a cause always occurs before its effect. Thus, we can describe a gene network by a directed graph based on the coefficient matrix B where nodes represent genes and edges show functional connectivities. In the estimated graph, node i is linked to node j by the edge i → j if and only if B_{ ij }≠ 0.
For simplicity of explanation, we introduce the following notations:
Since the covariance matrix Σ is diagonal, the first order VAR model of these observations can then be considered as m linear regression models
Here, we center each input variable so that there is no intercept in (2) and scale each input variable so that the observed standard deviation is equal to one.
The major objective of our analysis is to estimate the coefficient matrix B from observations. Especially, our interest is its structures, that is, we want to estimate which elements of each column of B are zero or nonzero. Let denote the set of indices corresponding the nonzero components of by . The graph inference problem can be formally stated as the problem of predicting from the observations. This formulation is equivalent to m variable selection problems in linear regression where each component in y_{ t }is the response and y_{t1}is the covariates.
Recursive Elastic Net
From this section, we focus on a variable selection problem in linear regression model with response y_{t, j}and covariates y_{t1,1}, ..., y_{t1, m}. For notational simplicity, we omit the suffix j of z_{ j }, β_{ j }, λ_{1j}and λ_{2j}, and use the notation z, β, λ_{1} and λ_{2}.
To obtain a sparse solution for parameter estimation with high prediction accuracy and encourage a group effect between variables in linear regression, Zou and Hastie [13] proposed the elastic net, which combines both l_{1}regularization and l_{2}regularization together. Similar to the lasso, due to the property of l_{1}penalty, the elastic net performs automatic variable selection. It can also choose highly correlated variables which is due to the benefit of l_{2}penalty. The loss function of the elastic net for parameter estimation can be represented by
where λ_{1} and λ_{2} are regularization parameters. The naive elastic net estimator is then the minimizer of (3):
Although the naive elastic net overcomes the limitations of the lasso, Zou and Hastie [13] showed the naive elastic net does not perform satisfactorily in numerical simulations and this deficiency in performance is due to imposing a double shrinkage with l_{1} and l_{2}penalties. To solve the problem, they proposed the elastic net estimator as a biascorrected version of the naive elastic net estimator given by:
However, the elastic net suffers from too many false positives in the variable selection problem for extremely highdimensional data, which will be illustrated in Simulation Section. This is due to the same amount of shrinkage for all the coefficients in the l_{1}penalty. That is, the same regularization parameter λ_{1} is used for each coefficient without assessing their relative importance. In a typical setting of linear regression, it has been shown that such an excessive l_{1}penalty can affect the consistency of model selection [15–18].
To enable different amounts of shrinkage for each coefficient, we define the following weighted loss function:
where λ_{1} and λ_{2} are regularization parameters, and w_{1}, ..., w_{ m }are positive coefficient weights. The weighted elastic net estimator and its biascorrected estimator are defined as:
and
The modification of the l_{1}penalty was first proposed by Zou [18]. This type of l_{1}regularization was applied to parameter estimation for several statistical models, such as Cox's proportional hazard model [19], least absolute deviation regression model [20], and graphical Gaussian model [21].
We next consider how to set the values of w_{ k }. One possible way is to select the weights so that they are inversely proportional to the true coefficients that is,
where L is some specified large value. If w_{ k }is small, tends to be nonzero. While tends to be zero if w_{ k }is large. By using the weight (9) and appropriate values of the regularization parameters, we can choose the true set of variables, i.e., . However, in many situations, one does not know true coefficients, and thus we need some plugin estimator instead of the true coefficients. In this article, we propose a new multistep regularization, called recursive elastic net. The proposed method tries to obtain a better estimator than the elastic net estimator by solving the weighted elastic net problem and updating the weights alternatively. The recursive elastic net can be described as:
Recursive Elastic Net

1.
Set a maximum number of iterations to be M and choose an initial estimator by the naive elastic net (4).

2.
For the lth iteration (l = 1, ..., M), define the coefficient weights by using the (l  1)th weighted elastic net estimator as follows:

3.
Solve the weighted elastic net problem and estimate the coefficient vector:
(10) 
4.
Repeat Step 2 and Step 3 until a stopping criterion is fulfilled or when the number of iterations attains M.
The parameter δ > 0 is introduced to avoid that the weights of zerovalued components in take the value of infinity. From numerical simulations which are not demonstrated, δ should be taken near 0. In this article, we set δ = 10^{5}. Following Zou and Hastie [13], we can also consider a biascorrected version of the recursive elastic net by replacing (10) with
We call the iterative procedure using the elastic net estimator (5) as the initial variable weights and updating the coefficients by the biascorrected weighted elastic net (11) corrected recursive elastic net. Note that our solution with the iterative procedure is a local optimum when initializing some variable weights and we cannot guarantee the convergence of a global minimum. Thus, it is important to choose a suitable starting point for the variable weights. In the recursive elastic net and the corrected recursive elastic net, we propose to initialize with the naive elastic net and the elastic net estimators, (4) and (5), that is, the unweighted solutions of (10) and (11), respectively. In the recursive elastic net and the corrected recursive elastic net, the regularization parameters determine whether or not the lth estimator is better than the (l  1)th estimator. The selection problem of the regularization parameters will be discussed in the latter section.
Grouping Effect
In this section, we show that the recursive elastic net estimator can lead to a desirable grouping effect for correlated variables. In the framework of linear regression, the grouping effect can be considered as an effect that the coefficients of highly correlated variables should be similar each other. As discussed in Zou and Hastie [13], the lasso does not have the grouping effect and tends to select one among a correlated set. The following Lemma which is quoted from Zou and Hastie [13] guarantees the grouping effect for the recursive elastic net in the situation with identical variables, since the loss function of the recursive elastic net is strictly convex when λ_{2} > 0. The Lemma is given by:
Lemma 1 Assume that x_{ i }= x_{ j }for i, j ∈ {1, ..., m} and is estimated by (10), then for any λ_{2} > 0.
Furthermore, the following theorem provides an upper bound on the difference of the ith and jth coefficients for the lth iteration with the proposed method:
Theorem 1 For given dataset (z, X), regularization parameters (λ_{1}, λ_{2}) and variable weights of the lth iteration , suppose that the response z is centered and the covariate matrix X is standardized. Let be the recursive elastic net estimator of the lth iteration (10). Suppose that two coefficients satisfy and define
then
where and is the sample correlation.
The proof of Theorem 1 is given by Appendix. We notice that the grouping effect of the recursive elastic net contributes not only the sample correlation ρ but also the difference between the ith and jth coefficients of the (l  1)th iterations . If x_{ i }and x_{ j }are highly correlated, i.e., ρ ≐ 1, and the difference is almost zero, then Theorem 1 says that the difference is also almost zero.
Computational Algorithm
In this section, we provide an algorithm for solving the weighted elastic net problem in (6) by modifying the LARSEN algorithm [13]. The modified algorithm called LARSWEN can be described as:
Algorithm: LARSWEN

1.
Define the (n  1) × m design matrix by

2.
Solve the elastic net problem by the LARSEN algorithm:

3.
Output the weighted lasso estimator by
For a fixed λ_{2}, the LARSWEN algorithm sequentially updates the locations of the nonzero coefficients and its values, and produces the entire solution path of the weighted elastic net estimator . The R source code is available from the authors upon request. The computational cost of the LARSWEN algorithm is the same as that of the LARSEN algorithm except for the treatment of Step 1 and Step 3. If the algorithm for calculating is stopped after u steps, it requires O(u^{3} + mu^{2}) operations.
Selection of Regularization Parameters
The regularization parameters play a critical role for determining the performance of variable selection with the recursive elastic net. We now discuss how to choose the values of the regularization parameters λ_{1}and λ_{2} in the recursive elastic net. Traditionally, crossvalidation is often used for estimating the predictive error and evaluating different models. For example, leaveoneout crossvalidation involves using a single observation from the original sample as the validation data and the remaining observations as the training data, and requires additional calculations for repeating such that each observation in the sample is used once as the validation data. Note that we need to choose two regularization parameters for each variable j in the VAR model. Thus, crossvalidation can be timeconsuming and hardly be applied for the regularization parameter selection in the estimation of the VAR model with the recursive elastic net. In order to choose the two regularization parameters, we use the degrees of freedom as a measure of model complexity. The profile of degrees of freedom clearly shows that how the model complexity is controlled by shrinkage, which helps us to pick an optimal model among all the possible candidates [22]. In recent years, an unbiased estimate of the degrees of freedom for the lasso [22] was proposed. Using the result of Zou et al. [22], an unbiased estimate of the degrees of freedom for the recursive elastic net in the lth iteration is given by
where denotes the active set of and the corresponding design matrix is given by
Note that does not include λ_{1} and depends on and λ_{2}. Let and be the coefficient vector and the active set of the LARSWEN step α, respectively, and λ_{α,1}be the corresponding l_{1}regularization parameter of the step α. Hence, we can use in place of λ_{α,1}as a tuning parameter.
In our numerical studies, two model selection criteria, Bayesian information criterion (BIC) [23] and biascorrected Akaike information criterion (AICc) [24, 25], are considered. Substituting the number of parameters in BIC and AICc by the degrees of the freedom we can define a modified BIC and a modified AICc for selecting the regularization parameters as
where
For a given λ_{2}, we choose the optimal step that gives the smallest BIC or the smallest AICc from the candidate models.
VAR Modeling with Recursive Elastic Net
As shown in the previous section, the inference of nonzero components in the VAR model can be formalized as m variable selection problems, i.e., our objective is to estimate . Therefore, we can apply the recursive elastic net directly to estimate the coefficient matrix B in the VAR model. In this section, we provide an algorithm for estimating B in the VAR model with the recursive elastic net and the derived model selection criteria for a fixed λ_{2}. The algorithm is as follows:
Algorithm:RENETVAR

1.
Set a maximum number of iterations to be M and a maximum step size to be u.

2.
For a fixed λ_{2}, start with the initial coefficient weights

3.
For j = 1, ..., m

(a)
For the lth iteration, given , the LARSWEN algorithm after u steps produces u solutions,
where

(b)
Calculate the values of model selection criterion for u candidate models, say
where denotes the active set of , and MC = AICc or BIC.

(c)
Select from among which satisfies

(d)
Update the variable weights by
(18) 
(a)

4.
Update l = l + 1 and repeat Step 3 until a stopping criterion is satisfied or the number of iterations attains M.
As a stopping criterion, we calculate the sum of the values of model selection criterion for j = 1, ..., m:
The algorithm stops if SMC^{(l)}(λ_{2})  SMC^{(l  1)}(λ_{2}) > 0 (l = 2, 3, ..., M) holds.
Typically, we first pick a grid of values for λ_{2}, say
For each λ_{γ, 2}, the RENETVAR algorithm produces . Finally, we select which satisfies
We can also estimate with the biascorrected recursive elastic net by replacing the solution (18) with
We note that our iterative algorithm guarantees the decrease of the sum of the values of model selection criterion and searches the model minimizing (19), although it cannot guarantee the decreasing of the loss function.
Results and Discussion
We tested the performance of the proposed method on both simulated and real gene expression timecourse data. Using simulated data, we showed the differences between the proposed method and other l_{1}regularization approaches in the variable selection problem for linear regression model and evaluated how well the proposed method estimates the coefficient matrix of the VAR model compared with other estimation methods. Using real data, we also compared the proposed method with various reverseengineering algorithms for inferring gene networks.
Simulation Results
Simulation 1
We compared our proposed methods, the recursive elastic net (REN) and the corrected recursive elastic net (CREN), with the other l_{1}regularization methods, the lasso (LA), the naive elastic net (NEN) and the elastic net (EN) on simulated data in linear regression. We generated data from a linear regression model
where β* was a 1000dimensional vector of coefficients, and x followed a 1000dimensional multivariate normal distribution N(0, S) with the zero mean vector 0 and the covariance matrix S whose (i, j)th entry was S_{i, j}= 0.3 if i ≠ j and S_{i, i}= 1 otherwise. The true coefficient vector was , i.e., the first five coefficients were related with the response z. The first five entries were generated from a uniform distribution on the interval [1.2, 0.8] and [0.8, 1.2]. The standard deviation σ was chosen so that the signaltonoise ratio was 10, and 50 observations were generated from this model. We set the grid of the values of the l_{2}regularization parameters to {0.5, 1, 2, 5}. The regularization parameters of each method were selected by a model selection criterion AICc. We used AICc as a stopping criterion of REN and CREN and the algorithm stopped if AICc was not decreasing in the next iteration. All the methods were run on each dataset through 100 simulations. Results of Simulation 1 are described in Table 1. TDR stands for true discovery rate (or accuracy) defined as TDR = TP/(TP + FP) and sensitivity (SE) is TP/(TP+FN), where TP, FP, TN and FN are the number of true positives, false positives, true negatives and false negatives, respectively. Compared with REN and CREN with the other l_{1}regularization methods, LA, NEN and EN had larger numbers of true positives than REN and CREN. However, LA, NEN and EN also had larger numbers of false positives, and they suffered from their low true discovery rates. While REN and CREN had the higher true discovery rates than LA, NEN and EN, since they reduced the numbers of false positives drastically while keeping the numbers of true positives as large as possible. Especially, REN reduced the number of false positives from 11.90 to 1.01 and had the highest true discovery rate between all the algorithms.
Next we illustrates the difference between the naive elastic net and the recursive elastic net on a dataset from the model (20). Figure 1 illustrates the coefficient profiles through 10 iterations with the application of the recursive elastic net on a dataset of 50 observations. The red lines indicate the first 5 coefficient profiles related with the response z, that is, the five true coefficient profiles. While the black lines indicate the noisy coefficient profiles. The dot line stands for the iteration when the stopping criterion was fulfilled.
The coefficients of the 0th iteration was equivalent to the naive elastic net estimators where 19 coefficients were identified as nonzero including the 5 true components (true positives) and 14 noisy components (false positives). In the first iteration, among the 19 components selected in the naive elastic net, 7 coefficients including the 5 true components were estimated as nonzero by the recursive elastic net. That is, we succeeded in reducing the number of false positives from 14 to 4 while keeping all the 5 true positives. After 5 iterations, the stopping criterion was fulfilled and the recursive elastic net yielded 5 true positives and only 1 false positive. This reduction of the number of false positives in the linear regression has much effect on the performance of the VAR model which is equivalent to m linear regressions. The performances of the proposed methods in the VAR model will be illustrated in the next section.
Simulation 2
For the second simulation, we used the VAR model with the structure of scalefree networks. We considered five models with (n, m) = (20, 100), (20, 200), (20, 500), (20, 1000) and (20, 2000) where n was the number of time points and m was the number of variables. The topologies of the scalefree networks were generated according to the BarabasiAlbert model [26] by using the R function barabasi.game in the R package igraph [27]. The parameter power was set to 1.2 and the other parameters were set to their default values, and then the numbers of edges were chosen to be 108, 213, 525, 1045 and 2079 for m = 100, 200, 500, 1000 and 2000, respectively. The structures of the simulated networks are illustrated in Additional File 1. The dataset from the VAR model with a given structure was drawn as follows:

1.
The nonzero entities for the true VAR coefficients B* were generated from a uniform distribution on the interval [0.8, 0.4] and [0.4, 0.8].

2.
The initial value y_{0} was drawn from a uniform distribution on the interval [10, 100].

3.
The data at time points t ∈ {1, 2, ..., 20} was generated recursively according to the VAR model with identity random noise matrix Σ = I.
We simulated 100 different datasets from the above procedure for each combination of (n, m). We compared the recursive elastic net (REN) and the corrected recursive elastic net (CREN) with the lasso (LA), the naive elastic net (NEN), the elastic net (EN), and the JameStein shrinkage (JS) [7] on simulated datasets. OpgenRhein and Strimmer [7] proposed an improved estimator of the VAR coefficients by using a JamesSteintype shrinkage approach. All the algorithms were run out on each drawn dataset {y_{1}, ..., y_{20}}. The setting of the l_{2}regularization parameters was same as given in Simulation 1. The regularization parameters of LA, NEN, EN, REN and CREN were selected by a model selection criterion, BIC or AICc. JS is a testbased approach and a cutoff value for the local false discovery was set to 0.2 (JSA), which was used in OpgenRhein and Strimmer [7]. We also used a cutoff value such that the number of significant edges detected by JS was as same as that of REN with AICc (JSB).
For each algorithm, its network inference performance was evaluated by true discovery rate (TDR = TP/(TP + FP) and sensitivity (SE = TP/(TP + FN)) where TP, FP, TN and FN are the number of true positives, false positives, true negatives and false negatives, respectively. Table 2 presents the results of Simulation 2. These results can be summarized as follows:

1.
The sensitivity and true discovery rate of each algorithm were gradually decreasing as increasing the number of variables.

2.
NEN and EN had higher sensitivities than the other methods, but their true discovery rates were much lower than those of REN and CREN.

3.
The sensitivity of LA was slightly higher than that of REN, while the true discovery rate of REN was over ten times higher than that of LA.

4.
Compared with the same significant edges of REN and JSB, REN outperformed JSB in both of sensitivity and true discovery rate.

5.
CREN succeeded in higher sensitivity than REN, but it suffered from its low true discovery rate compared with REN when the number of variables was large.

6.
REN had the highest true discovery rate between all the algorithms in all cases while its sensitivity was almost same as that of LA, JSA and CREN.

7.
REN with AICc was slightly better than REN with BIC.
We also compared our proposed method with other competing methods on simulated dataset by using the VAR model with hierarchical structures. These results are given in Additional File 2.
Experimental Results
In this section we used experimental data of MCF7 breast cancer cells stimulated with two ErbB ligands, epidermal growth factor (EGF) and heregulin (HRG), to compare the proposed method with various reverseengineering algorithms. This dataset is available from Gene Expression Omnibus (GEO, Accession Number GSE6462) [28]. In GSE6462, cells were stimulated with 0.1, 0.5, 1, or 10 nM of either EGF or HRG and the expression values were measured at eight time points (0, 5, 10, 15, 30, 45, 60, and 90 minutes) by Affymetrix GeneChip U133A2. The expressions with 10 nM of EGF at 60 minutes were not observed, and thus its values were estimated with linear interpolation. These microarray datasets were also normalized by faster cyclic loess [29].
In addition to the lasso, the naive elastic net, the elastic net, and the JamesStein shrinkage which were used in Simulation Section, we considered MRNET [30], CLR [2], and ARACNE [31] as competing methods. These methods are available from the R package minet [32]. Note that they do not intend to analyze timecourse data. To increase their capabilities, we considered a lagged version of a time series matrix, denoted by D, whose (i, j)element indicates D_{i, j}= y_{i+1, j} y_{i, j}and used it as an input of MRNET, CLR and ARACNE. This modification enables us to extract genegene interactions between time steps i + 1 and i. To distinguish between the original version and the modified version, we refer the original MRNET, CLR and ARACNE as MRN, CLR and ARA, and the modified MRNET, CLR and ARACNE as MRNL, CLRL and ARAL, respectively. In our analysis, we assumed that time points were equidistant and that the structure of gene networks did not change with different dosages of EGF or HRG. Thus, the datasets from different dose conditions of EGF or HRG were considered as replicated observations. We validated the performances of these reverseengineering algorithms in terms of how to identify pairwise genegene interaction edges by using TRANSPATH database [33]. Although we do not have perfect knowledge of true network, we used biological interactions obtained from TRANSPATH database as true edges. We started by focusing on 254 probes by Nagashima et al. [34]. Among these probes, 233 probes were mapped to unique Entrez Gene IDs, 7 probes were mapped to more than one Entrez Gene IDs, and 14 probes were not mapped, respectively. We only used 233 probes with 228 unique Entrez Gene IDs. To correspond each Entrez Gene ID to unique probe, we selected a probe for each Entrez Gene ID with the largest mean value of expressions through all the experiments. Among these 228 Entrez Gene IDs, TRANSPATH database identified a biological pathway which includes 110 Entrez Gene IDs. While the other 128 Entrez Gene IDs do not connect each other in TRANSPATH database. Thus, all the algorithms were run on timecourse expression profiles of the 110 Entrez Gene IDs for EGFstimulation and HRGstimulation, respectively.
The setting of the regularization parameters in the lasso (LA), the naive elastic net (NEN), the elastic net (EN), the recursive elastic net (REN) and the corrected recursive elastic net (CREN) was same as given in Simulation and these parameters were selected by BIC. As a result, LA, NEN, EN, REN and CREN identified 2202, 1721, 870, 162 and 436 significant edges from the EGF dataset, respectively. While they also identified 2904, 2335, 411, 154 and 171 significant edges from the HRG dataset, respectively. In order to facilitate comparison among them, we controlled the numbers of significant edges detected by LA, NEN, EN and CREN same as that of REN by using a threshold value, that is, set to 0 all edges whose coefficients are lower than the threshold. We also chose the same number of significant edges detected by the JamesStein shrinkage (JS) as that of REN in a similar way as in the previous section. MRN, CLR, ARA, MRNL, CLRL and ARAL require the calculation of a mutual information matrix. We used a MillerMadow asymptotic bias corrected empirical estimator as each mutual information of them. We set the numbers of significant edges detected by them as same as that of REN by using a threshold. Notice that each of MRN, CLR, ARA, MRNL, CLRL and ARAL infers the network just as an undirected graph, while each of the others infers the network as a directed graph. In our analysis, we first transformed the real and the inferred directed networks into undirected graphs and then performed a comparison. Results from experimental data are described in Table 3. TDR stands for true discovery rate defined as TP = (TP + FP) where TP and FP are the numbers of true positives and false positives, respectively. The method RAND is an algorithm which selects the same number of significant edges same as that of REN randomly. Each value of RAND indicates the mean value through 1,000,000 simulations. The pvalue of each algorithm was calculated using a Binomial distribution based on the random model through 1,000,000 simulations. Notice that the overall low performances on the EGF and HRG datasets are owing to the limited number of data and the imperfect knowledge of the real network. We observed that, in the EGF dataset, ARA, MRNL, LA, NEN, REN and CREN had higher numbers of true positives than RAND, however only the performance of REN was significantly better than that of RAND. While, in the HRG dataset, MRNL, CLRL, ARAL, NEN, EN, REN and CREN performed well than RAND. Between them, the performances of EN, REN and CREN were significantly better than that of RAND. As a result, REN had the highest number of true positives on both of the EGF and HRG datasets.
Next we applied our method to timecourse expression profiles of the 228 Entrez Gene IDs stimulated with EGF and HRG, respectively, in order to investigate the difference between the EGF and HRGinduced gene networks. In the ErbB signaling network, two different ErbB receptor ligands, EGF and HRG, play a key role in controlling diverse cell fates. In MCF7 cells, stimulation with EGF causes sustained network activation and leads to cell differentiation, while stimulation with HRG causes transient network activation and leads to cell proliferation.
Figures 2 and 3 show the inferred EGF and HRGinduced gene networks with the recursive elastic net.
These figures were illustrated by Cell Illustrator [35, 36]. We observed that the two directed graphs had different network topologies each other. We compared with the two graphs from characteristics of node degree, that is, the number of edges that they had. In directed networks, we distinguish between the indegree, the number of directed edges that point toward the node, and the outdegree, the number of directed edges that start at the node. The indegree distributions and the outdegree distributions of the EGF and HRGinduced networks are illustrated in Additional File 3. We found that each of the indegree distributions concentrated around its average value, whereas the outdegree distributions had long tails. In particular, the outdegree distribution of the HRGinduced network was more heavy tailed than that of the EGFinduced network. Furthermore, in the HRGinduced network, we found that three hub genes, FOS, MAP3K8 and PDLIM5, had larger outdegrees (30, 68 and 75, respectively) than the other genes. While, in the EGFinduced network, each of FOS, MAP3K8 and PDLIM5 had only one outgoing edge (one outdegree). Thus, these hub genes are thought to be essential in only the HRGinduced network since their disruptions do not lead to a major loss of connectivity in the EGFinduced network, but the loss of them causes the breakdown of the HRGinduced network into isolated clusters. Our analysis based on the VAR model with the proposed method produced the hypothesis that these hub genes might be implicated why only HRG allows MCF7 cells to be differentiated.
Conclusion
In this article, we have proposed a new class of l_{1}regularization, called recursive elastic net for inferring gene networks from timecourse microarray data based on the VAR model. The recursive elastic net addresses the drawback of the elastic net and has higher true discovery rate than other competing methods in gene network inference. Numerical simulations demonstrated that the proposed method succeeded in reducing the number of false positives drastically while keeping the large number of true positives and achieved two or more times higher true discovery rate than the lasso, the naive elastic net, the elastic net and the JamesStein shrinkage even when the number of time points was small. We also compared our method with various reverseengineering algorithms including MRNET, CLR and ARACNE by using experimental data of MCF7 breast cancer cells and TRANSPATH database. As a result, we found that the proposed method had the best performance between them and provided some differences between the EGF and HRGinduced networks.
An interesting direction for future work would be to incorporate prior information elicited from the known biological networks. In the recursive elastic net, we can incorporate the topology of the known network directly as the initial variable weights (10). For example, if an edge from gene i to gene j exists the known biological pathway, we impose a small weight less than one on the corresponding coefficient. This incorporation would further provide higher sensitivity and higher true discovery rate of the proposed method in gene network inference. Another promising direction for future work would be the extension of the recursive elastic net to other type of mathematical models. A limitation of the VAR model comes from their equidistant assumption which is not true in many real situations. We might solve the problem by applying the proposed method to differential equation models [37, 38].
Appendix: Proof of Theorem 1
We provide proof of Theorem 1 which is similar to that of Zou and Hastie [13]. If , then both and are nonzero, and we have . Because of (10), the estimator satisfies
Hence we have
Subtracting (21) from (22) gives
which is equivalent to
Using the identity A ± B ≤ A + B for all A and B, we have
By equation (6) we must have
that is,
So
Then equation (23) implies
If we assume that X is standardized, where . Thus inequality (24) reduces to
We assume that . Since the function f(x) = x^{1} is Lipshitz continuous for x > 0, then we have
Using (26) gives
Thus inequality (25) further simplifies to
This completes the proof.
References
 1.
Di Bernardo D, Thompson MJ, Gardner TS, Chobot SE, Eastwood EL, Wojtovich AP, Elliott SJ, Schaus SE, Collins JJ: Chemogenomic profiling on a genomewide scale using reverseengineered gene networks. Nat Biotechnol. 2005, 23: 377383.
 2.
Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ, Gardner TS: Largescale mapping and validation of Escherichia coli transcriptional regulation from a Compedium of expression profiles. PLoS Biol. 2007, 5: e8.
 3.
Gardner TS, di Bernardo D, Lorenz D, Collins JJ: Inferring genetic networks and identifying compound mode of action via expression profiling. Science. 2003, 301: 102105.
 4.
Imoto S, Tamada Y, Araki H, Yasuda K, Print CG, CharnockJones SD, Sanders D, Savoie CJ, Tashiro K, Kuhara S, Miyano S: Computational strategy for discovering druggable gene networks from genomewide RNA expression profiles. Pac Symp Biocomput. 2006, 559571.
 5.
Tamada Y, Imoto S, Tashiro K, Kuhara S, Miyano S: Identifying drug active pathways from gene networks estimated by gene expression data. Genome Informatics. 2005, 16: 182191.
 6.
Fujita A, Sato JR, GarayMalpartida HM, Yamaguchi R, Miyano S, Sogayar MC, Ferreira CE: Modeling gene expression regulatory networks with the sparse vector autoregressive model. BMC Syst Biol. 2007, 1: 39.
 7.
OpgenRhein R, Strimmer K: Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process. BMC Bioinformatics. 2007, 8 (Suppl 2): S3.
 8.
George EI: The variable selection problem. J Am Statist Assoc. 2000, 95: 13041308.
 9.
Miller A: Subset Selection in Regression. 2002, Chapman & Hall/CRC, Second.
 10.
Tibshirani R: Regression shrinkage and selection via the lasso. J R Statist Soc B. 1996, 58: 267288.
 11.
Fan J, Lv J: Sure independence screening for ultrahigh dimensional feature space. J R Statist Soc B. 2008, 70: 849911.
 12.
Breiman L: Heuristics of instability and stabilization in model selection. The Annals of Statistics. 1996, 24: 23502383.
 13.
Zou H, Hastie T: Regularization and variable selection via the elastic net. J R Statist Soc B. 2005, 67: 301320.
 14.
Granger CWJ: Investigating causal relation by econometric and crosssectional method. Econometrica. 1969, 37: 424438.
 15.
Leng C, Lin Y, Wahba G: A note on lasso and related procedures in model selection. Statistica Sinica. 2006, 16: 12731284.
 16.
Wang H, Leng C: A note on adaptive group lasso. Comput Statist Data Anal. 2008, 52: 52775286.
 17.
Yuan M, Lin Y: On the nonnegative garrote estimator. J R Statist Soc B. 2007, 69: 143161.
 18.
Zou H: The adaptive lasso and its oracle properties. J Am Statist Assoc. 2006, 101: 14181429.
 19.
Zhang HH, Lu W: Adaptive lasso for Cox's proportional hazards model. Biometrika. 2006, 94: 691703.
 20.
Wang H, Guodong L, Guohua J: Robust regression shrinkage and consistent variable selection through LADlasso. J Bus Econ Statist. 2007, 25: 347355.
 21.
Shimamura T, Imoto S, Yamaguchi R, Miyano S: Weighted lasso in graphical Gaussian modeling for large gene network estimation based on microarray data. Genome Informatics. 2007, 19: 142153.
 22.
Zou H, Hastie T, Tibshirani R: On the "degrees of freedom" of the lasso. Ann Statist. 2007, 35: 21732192.
 23.
Schwarz G: Estimating the dimension of a model. Ann Statist. 1978, 6: 461464.
 24.
Hurvich CM, Tsai CL: Bias of the corrected AIC criterion for underfitted regression and time series models. Biometrika. 1991, 78: 521530.
 25.
Sugiura N: Further analysis of the data by Akaike's information criterion and the finite corrections. Comm Statist A. 1978, 7: 1326.
 26.
Barabási AL, Albert R: Emergence of scaling in random networks. Science. 1999, 286: 509512.
 27.
Csárdi G, Nepusz T: igraph Reference Manual. http://cneurocvs.rmki.kfki.hu/igraph/
 28.
NCBI Gene expression omnibus Accession Number GSE6462. http://www.ncbi.nlm.nih.gov/geo/
 29.
Ballman KV, Grill DE, Oberg AL, Therneau TM: Faster cyclic loess: normalizing RNA arrays via linear models. Bioinformatics. 2004, 20: 27782786.
 30.
Meyer PE, Kontos K, Lafitte F, Bontempi G: Informationtheoretic inference of large transcriptional regulatory networks. EURASIP J Bioinform Syst Biol. 2007, 2007: 79879.
 31.
Basso K, Margolin AA, Stolovitzky G, Klein U, DallaFavera R, Califano A: Reverse engineering of regulatory networks in human B cells. Nat Genet. 2005, 37: 382390.
 32.
Meyer PE, Lafitte F, Bontempi G: minet: A R/Bioconductor package for inferring large transcriptional networks using mutual information. BMC Bioinformatics. 2008, 9: 461.
 33.
TRANSPATH: the pathway databases. http://www.biobase.de/pages/index.php?sid=transpathdatabases
 34.
Nagashima T, Shimodaira H, Ide K, Nakakuki T, Tani Y, Takahashi K, Yumoto N, Hatakeyama M: Quantitative transcriptional control of ErbB receptor signaling undergoes graded to biphasic response for cell differentiation. J Biol Chem. 2007, 282: 40454056.
 35.
Cell Illustrator. http://cionline.hgc.jp/
 36.
Nagasaki M, Doi A, Matsuno H, Miyano S: Genomic Object Net: I. A platform for modeling and simulating biopathways. Applied Bioinformatics. 2003, 2: 181184.
 37.
Yeung MK, Tegne'r J, Collins JJ: Reverse engineering gene networks using singular value decomposition and robust regression. Proc Natl Acad Sci USA. 2002, 99: 61636168.
 38.
Wang Y, Joshi T, Zhang XS, Xu D, Chen L: Inferring gene regulatory networks from multiple microarray datasets. Bioinformatics. 2006, 22: 24132420.
Acknowledgements
We are grateful to Dr. Rainer OpgenRhein for providing the R program code of the JamesStein procedure and to anonymous reviewers for many constructive comments and suggestions. This work was supported by the NextGeneration Supercomputer Project under RIKEN and MEXT, Japan. The computational resource was also provided by the Super Computer System, Human Genome Center, Institute of Medical Science, University of Tokyo.
Author information
Affiliations
Corresponding author
Additional information
Authors' contributions
TS conceived and designed the study, and wrote the manuscript. SI provided statistical expertise and careful manuscript review. RY, AF and MN assisted in preparing the manuscript. SM supervised the whole project. All authors read and approved of the final manuscript.
Electronic supplementary material
Simulation 3
Additional file 2: . This document represents Simulation 3 where we compared our proposed method with other competing methods on simulated dataset by using the VAR model with a hierarchical structure. It also includes Additional Figure 6 that describes the structure of the simulated network and Additional Tables 1 and 2 that show the results of Simulation 3, respectively. (PDF 275 KB)
12918_2008_309_MOESM3_ESM.pdf
Additional file 3: Indegree and outdegree distributions of the EGF and HRGinduced gene networks. This file includes Additional Figures 7 and 8 that describe the indegree distribution and the outdegree distribution of the EGFinduced gene network, and Additional Figures 9 and 10 that describe the indegree distribution and the outdegree distribution of the HRGinduced gene network, respectively. (PDF 72 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Shimamura, T., Imoto, S., Yamaguchi, R. et al. Recursive regularization for inferring gene networks from timecourse gene expression profiles. BMC Syst Biol 3, 41 (2009). https://doi.org/10.1186/17520509341
Received:
Accepted:
Published:
Keywords
 Regularization Parameter
 Lasso
 Model Selection Criterion
 Variable Selection Problem
 True Discovery Rate