Recursive regularization for inferring gene networks from time-course gene expression profiles

Vector Autoregressive Model

We now consider gene expression data in time-course 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_t-1at the previous time step t - 1 with a first-order VAR model:

(1)

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 m-dimensional 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

(2)

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_t-1is 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_t-1,1, ..., y_{t-1, m}. For notational simplicity, we omit the suffix j of z _j , β _j , λ_1jand λ_2j, and use the notation z, β, λ₁ and λ₂.

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₁-regularization and l₂-regularization together. Similar to the lasso, due to the property of l₁-penalty, the elastic net performs automatic variable selection. It can also choose highly correlated variables which is due to the benefit of l₂-penalty. The loss function of the elastic net for parameter estimation can be represented by

(3)

where λ₁ and λ₂ are regularization parameters. The naive elastic net estimator is then the minimizer of (3):

(4)

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₁- and l₂-penalties. To solve the problem, they proposed the elastic net estimator as a bias-corrected version of the naive elastic net estimator given by:

(5)

However, the elastic net suffers from too many false positives in the variable selection problem for extremely high-dimensional data, which will be illustrated in Simulation Section. This is due to the same amount of shrinkage for all the coefficients in the l₁-penalty. That is, the same regularization parameter λ₁ 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₁-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:

(6)

where λ₁ and λ₂ are regularization parameters, and w₁, ..., w _m are positive coefficient weights. The weighted elastic net estimator and its bias-corrected estimator are defined as:

(7)

and

(8)

The modification of the l₁-penalty was first proposed by Zou [18]. This type of l₁-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,

(9)

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 plug-in estimator instead of the true coefficients. In this article, we propose a new multi-step 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. 1.

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

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

   
    
3. 3.
   Solve the weighted elastic net problem and estimate the coefficient vector:

   

   (10)
    
4. 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 zero-valued 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 bias-corrected version of the recursive elastic net by replacing (10) with

(11)

We call the iterative procedure using the elastic net estimator (5) as the initial variable weights and updating the coefficients by the bias-corrected 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 l-th 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 λ₂ > 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 λ₂ > 0.

Furthermore, the following theorem provides an upper bound on the difference of the i-th and j-th coefficients for the l-th iteration with the proposed method:

Theorem 1 For given dataset (z, X), regularization parameters (λ₁, λ₂) and variable weights of the l-th iteration , suppose that the response z is centered and the covariate matrix X is standardized. Let be the recursive elastic net estimator of the l-th iteration (10). Suppose that two coefficients satisfy and define

(12)

then

(13)

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 i-th and j-th 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 LARS-EN algorithm [13]. The modified algorithm called LARS-WEN can be described as:

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 λ₁and λ₂ in the recursive elastic net. Traditionally, cross-validation is often used for estimating the predictive error and evaluating different models. For example, leave-one-out cross-validation 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, cross-validation can be time-consuming 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 l-th iteration is given by

(14)

where denotes the active set of and the corresponding design matrix is given by

Note that does not include λ₁ and depends on and λ₂. Let and be the coefficient vector and the active set of the LARS-WEN step α, respectively, and λ_α,1be the corresponding l₁-regularization parameter of the step α. Hence, we can use in place of λ_α,1as a tuning parameter.

In our numerical studies, two model selection criteria, Bayesian information criterion (BIC) [23] and bias-corrected 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

(15)

(16)

where

(17)

For a given λ₂, 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 λ₂. The algorithm is as follows:

Algorithm:RENET-VAR

1. 1.

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

    
2. 2.
   For a fixed λ₂, start with the initial coefficient weights

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

   1. (a)

      For the l-th iteration, given , the LARS-WEN algorithm after u steps produces u solutions,

       

   where

   1. (b)

      Calculate the values of model selection criterion for u candidate models, say

       

   where denotes the active set of , and MC = AICc or BIC.

   1. (c)

      Select from among which satisfies

       
   2. (d)

      Update the variable weights by

       

   

   

   (18)

   

   

   
    
4. 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:

(19)

The algorithm stops if SMC^(l)(λ₂) - SMC^{(l - 1)}(λ₂) > 0 (l = 2, 3, ..., M) holds.

Typically, we first pick a grid of values for λ₂, say

For each λ_{γ, 2}, the RENET-VAR algorithm produces . Finally, we select which satisfies

We can also estimate with the bias-corrected 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.

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₁-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

(20)

where β* was a 1000-dimensional vector of coefficients, and x followed a 1000-dimensional 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 signal-to-noise ratio was 10, and 50 observations were generated from this model. We set the grid of the values of the l₂-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₁-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.

Method	TP	FP	TN	FN	TDR	SE
LA	4.60	25.38	969.62	0.40	0.22	0.92
NEN	4.38	11.90	983.10	0.62	0.30	0.88
EN	4.22	7.89	987.11	0.78	0.40	0.84
REN	4.07	1.01	993.99	0.93	0.84	0.81
CREN	3.78	1.05	993.95	1.22	0.84	0.76

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 0-th 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.

Experimental Results

In this section we used experimental data of MCF-7 breast cancer cells stimulated with two ErbB ligands, epidermal growth factor (EGF) and heregulin (HRG), to compare the proposed method with various reverse-engineering 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 U133A-2. 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 James-Stein 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 time-course 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 gene-gene 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 MRN-L, CLR-L and ARA-L, 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 reverse-engineering algorithms in terms of how to identify pairwise gene-gene 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 time-course expression profiles of the 110 Entrez Gene IDs for EGF-stimulation and HRG-stimulation, respectively.

Next we applied our method to time-course expression profiles of the 228 Entrez Gene IDs stimulated with EGF and HRG, respectively, in order to investigate the difference between the EGF- and HRG-induced 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 MCF-7 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 HRG-induced 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 in-degree, the number of directed edges that point toward the node, and the out-degree, the number of directed edges that start at the node. The in-degree distributions and the out-degree distributions of the EGF- and HRG-induced networks are illustrated in Additional File 3. We found that each of the in-degree distributions concentrated around its average value, whereas the out-degree distributions had long tails. In particular, the out-degree distribution of the HRG-induced network was more heavy tailed than that of the EGF-induced network. Furthermore, in the HRG-induced network, we found that three hub genes, FOS, MAP3K8 and PDLIM5, had larger out-degrees (30, 68 and 75, respectively) than the other genes. While, in the EGF-induced network, each of FOS, MAP3K8 and PDLIM5 had only one outgoing edge (one out-degree). Thus, these hub genes are thought to be essential in only the HRG-induced network since their disruptions do not lead to a major loss of connectivity in the EGF-induced network, but the loss of them causes the breakdown of the HRG-induced 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 MCF-7 cells to be differentiated.