Simulation example

where + and − indicate the existence and regulation types of the edge. We randomly generate m stable regulatory matrices A(k), k = 1, . . . , m, according to the template A₀, such that sign(A(k)) = sign(A₀). Then, m simulated time-course gene expression datasets, each with the number of time points, n _k , are generated from (5) with randomly chosen expression values at the first time point. The simulated error follows a mixed Gaussian distribution: with probability of 0.8, it has the distribution N (0, 1) and with probability of 0.2, it has the distribution N (0, 102). In this way, the simulated data contain large errors and outliers. To investigate the performances of our methods in different situations, we vary the values of m and n _k and apply the group LASSO and Huber group LASSO respectively to these generated datasets and compare the results from these two methods.

Data are generated under three situations (m = 8, n _k = 15), (m = 4, n _k = 15) and (m = 4, n _k = 8). Using the stability selection procedure that is introduced in the Method Section, network typologies consisting of edges with scores or probabilities are inferred by Huber group LASSO and group LASSO. For the first two situations, we set the number of bootstrap samples as 30 and the moving block length as 10. For the third situation, we set the number of bootstrap samples as 30 and the moving block length as 5. Varying the threshold, the ROC plots and precision-recall plots of each method for different situations are obtained and illustrated in Figure 1. The areas under the ROCs (AUROCs) and precision-recall curves (AUPRs) are calculated and reported in Table 1. From Figure 1 and Table 1 we can see that for each situation, the Huber group LASSO outperforms the group LASSO, i.e. the AUROC and AUPR of Huber group LASSO are larger than those of group LASSO. ROC plots in Figure 1 also show that both methods have better performances than the random guess. For the case of m = 8 and n _k = 15, the Huber group LASSO even achieves the maximum value of AUROC and AUPR. It can also be seen that for each method, the more the observations or the more the datasets, the larger AUROC and AUPR can be obtained. This is in accord with the intuition because, in this example, more observations or datasets indicate more information as these simulated data are generated under quite similar circumstances. All the simulation results have demonstrated the effectiveness of our proposed method.

In vivo reverse engineering and modeling assessment (IRMA) data

The data used in this example come from the In vivo Reverse Engineering and Modeling Assessment (IRMA) experiment [17], where a network composed of five genes (GAL80, GAL4, CBF1, ASH1 and SWI5) was synthesized in yeast Saccharomyces cerevisiae, in which genes regulate each other through a variety of regulatory interactions. The network is negligibly affected by endogenous genes and it is responsive to small molecules. Galactose and glucose are respectively used to switch on and off the network. In this study, we use the IRMA time-course data consisting of four switch off datasets (with the number of time points varying from 19 to 21) and five switch on datasets (with the number of time points varying from 11 to 16).

The Huber group LASSO and the group LASSO are applied to these data in three cases: (1) switch on datasets, (2) switch off datasets and (3) all datasets, i.e., combing switch on and switch off datasets. In the stability selection procedure, the number of bootstrap samples is 30 for all cases and the moving block length is 14 for the second case and 8 for the other cases. The ROC plots and precision-recall plots for the Huber group LASSO and the group LASSO for each case are illustrated in Figure 2 and the corresponding AUROCs and AUPRs are summarized in Table 2. It can be seen that except the group LASSO for the switch on datasets, the performances of the methods are better than random guesses. The Huber group LASSO outperforms the group LASSO in both AUROCs and AUPRs. All methods for the switch off datasets perform better than for the switch on datasets. The group LASSO for all datasets has better performance than for the switch on datasets but is not as good as for the switch off datasets. The Huber group LASSO for all datasets has the best performance among all cases. This indicates that combing multiple datasets may lead to either the best result or a robust result which is better than the worst case. The network topology with false positive rate (FPR) 0.08 of the Huber group LASSO for all datasets is shown in Figure 3 and the corresponding true positive rate (TPR) is 0.75 with precision 0.86, in which the red edges represent true positives while black edges are false positives. The results show the effectiveness of our method for the IRMA data.

E. coli SOS network

In this example, we apply the proposed method to identify the real GRN, E. coli SOS DNA repair system as shown in Figure 4. This network is responsible for repairing the DNA after some damage happens. LexA acts as the master repressor of many genes in the normal states. When a damage occurs, RecA acts as a sensor and binds to single-stranded DNA to sense the damage and mediates the autocleavage of LexA. The repressions of the SOS genes are halted by the drop in LexA levels. The SOS genes are activated and start to repair the damages. When the repair is done, RecA level drops and stops mediating the autocleavage of LexA. Then, LexA accumulates and represses the SOS genes to make the cell go back to the normal state.

As the conditions for the first two experiments are different for the last two experiments, we consider applying the method to three cases: (1) datasets of experiment 1 and 2, (2) datasets of experiment 3 and 4 and (3) all experiment datasets. In the stability selection procedure, the number of bootstrap samples is 30 and the moving block length is 25 for all cases. The ROC plots and precision-recall plots for the Huber group LASSO and the group LASSO for each case are illustrated in Figure 5 and the corresponding AUROCs and AUPRs are illustrated in Table 4. From the ROC plots and AUROCs, it can be seen that the Huber group LASSO performs significantly better than random guess while the group LASSO method is only a little bit better than random guess. Obviously, the Huber group LASSO outperforms the group LASSO both in AUROCs and AUPR for all cases. The Huber group LASSO using experiment 3 and 4 datasets has the best performance. Performance of the Huber group LASSO using all datasets is between that in the first case and that in the second case. It can be considered as a robust result because of the using of multiple datasets. The network topology with FPR 0 and TPR 0.59 of the Huber group LASSO for all datasets is shown in Figure 6, in which all inferred edges are correct. These results demonstrate the effectiveness of our method for the E. coli SOS data.

S. cerevisae cell cycle subnetwork

A cell cycle regulatory subnetwork in S. cerevisae is inferred by the proposed method from 5 experimental microarray datasets. As in [22], the subnetwork consists of 27 genes including 10 genes for producing transcription factors (ace2, fkh1, swi4, swi5, mbp1, swi6, mcm1, fkh2, ndd1, yox1) and 17 genes for producing cyclin and cyclin/CDK regulatory proteins (cln1, cln2, cln3, cdc20, clb1, clb2, clb4, clb5, clb6, sic1, far1, spo12, apc1, tem1, gin4, swe1 and whi5). The time-course datasets we use include cell-cycle alpha factor release, cdc15, alpha factor fkh1 fkh2, fkh1,2 alpha factor and Elutriation, which are all downloaded from Stanford Microarray Database (SMD). We apply the Huber group LASSO and the group LASSO respectively to infer the network from the datasets.

In order to demonstrate the effectiveness of the proposed method, the inferred results are compared with the interaction network of the chosen 27 genes, drawn from BioGRID database [23]. The network in the database has 112 interactions, not including the self-regulations, and we take it as the gold standard regulatory network. In the stability selection procedure, the number of bootstrap samples is 30 and the moving block length is 9. The ROC plots and precision-recall plots are illustrated in Figure 7 and the AUROCs and AUPRs are shown in Table 5. We can see that both methods have better performances than random guess and the Huber group LASSO outperforms the group LASSO. One network from Huber group LASSO with FPR 0.43 and TPR 0.59 is shown in Figure 8, in which red edges are those inferred edges having been identified in the database and grey edges might be either false positives or novel discovered regulatory relations. Although the gold standard network extracted from the database may contain false edges or not be complete, this shows the effectiveness of our method to some extent.

Huber group LASSO

where A _i (k) ^T is the i th row of the matrix A(k) and x _j (k) is the j th column of the matrix X(k). w _k is the weight for the k th dataset, which can be assigned by experience. In this study, we choose w _k = n _k / Σn _k i.e., the more observations the dataset has, the higher weight it is assigned with. The penalty term in (6) takes advantage of the sparse nature of GRNs and has the effect making the grouped parameters to be estimated either all zeros or all non-zeros [10], i.e., a _iℓ (k)'s, k = 1, . . . , m, become either all zeros or all non-zeros. Therefore, a consistent network topology can be obtained from the group LASSO method. λ is a tuning parameter which controls the degree of sparseness of the inferred network. The larger the value of λ, the more grouped parameters become zeros.

The squared error and Huber loss function are illustrated in Figure 9. It can be seen that for small errors, these two loss functions are exactly the same while for large errors, Huber loss which increases linearly is less than the squared error loss which increases quadratically. Because Huber loss penalizes much less than the squared error loss for large errors, the Huber group LASSO is more robust than group LASSO when there exists large noise or outliers in the data. It is also known that the Huber loss is nearly as efficient as squared error loss for Gaussian errors [24].

where x _j is the j th column of X, y _ij is the j th element of Y _i , $n={\sum }_{k=1}^m{n}_k$ and ${\omega }_i=w_{1}I\left(i\le n_1\right)+{\sum }_{k=2}^m{w}_{k}I\left({\sum }_{l=1}^{k=1}{n}_l<i\le {\sum }_{l=1}^k{n}_l\right)$. (6) can be rewritten similarly.

Optimization algorithm

where for notational convenience, we omit the subscript i here and b _ℓ is a block of parameters of b, i.e. $b={\left[b_1^T,\dots ,b_p^T\right]}^T$.

where γ is the largest eigenvalue of ${\sum }_{j=1}^n{\omega }_j{x}_j{x}_j^T$. It can be easily shown that this auxiliary function satisfies (11).

Note that the algorithm can be adapted to solve (6) with quite similar derivations. In the following section, we show that the sequence {b^(k)} generated from the algorithm guarantees the objective function J (b) keep decreasing. We also show that the limit point of the sequence generated is indeed the minimum point of J (b).

Convergence analysis

The convergence of the optimization algorithm for the minimization of (10) is analyzed in the way similar to [25]. We first show the descent property of the algorithm.

Lemma 1 The sequence {b(k)} generated from the optimization algorithm keeps the objective function J(b) decreasing, i.e., J(b^(k)) ≥ J (b^(k+1)).

Proof By (11) and (13), we have

$\begin{array}{cc}J\left(b^{\left(k\right)}\right)\hfill & =Q\left(b^{\left(k\right)}|b^{\left(k\right)}\right)\ge Q\left(b^{\left(k\right)}\left(1\right)|b^{\left(k\right)}\right)\hfill \\ \ge Q\left(b^{\left(k\right)}\left(2\right)|b^{\left(k\right)}\left(1\right)\right)\ge \cdots \ge Q\left(b^{\left(k\right)}\left(p\right)|b^{\left(k\right)}\left(p-1\right)\right)\ge J\left(b^{\left(k\right)}\left(p\right)\right)\hfill \\ =J\left(b^{\left(k+1\right)}\right).\hfill \end{array}$

Next, we show that if the generated sequence satisfies some conditions, it converges to the optimal solution.

Lemma 2 Assume the data (y, X) lies on a compact set and the following conditions are satisfied:

1 The sequence {b^(k)} is bounded.

2 For every convergent subsequence $\left\{b^{\left(n_k\right)}\right\}\subset \left\{b^{\left(k\right)}\right\}$, the successive differences converge to zeros, $b^{\left(n_k\right)}-b^{\left(n_k-1\right)}\to 0$.

Then, every limit point b^∞ of the sequence {b^(k)} is a minimum for the function J(b), i.e., for any $\delta ={\left({\delta }_1^T,\dots ,{\delta }_p^T\right)}^T\in {ℝ}^{mp}$,

since $b_j^T{\delta }_j\le {∥b_j∥}_2{∥{\delta }_j∥}_2$.

for any 1 ≤ j ≤ p.

For $\delta ={\left({\delta }_1^T,\dots ,{\delta }_p^T\right)}^T\in {ℝ}^{mp}$, due to the differentiability of f (b),

$\begin{array}{cc}{\displaystyle \underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}}\left\{\frac{J\left(b^{\infty }+\alpha \delta \right)-J\left(b^{\infty }\right)}{\alpha }\right\}\hfill & ={\displaystyle \sum _{j=1}^p}{\nabla }_{j}f{\left(b^{\infty }\right)}^T{\delta }_j+{\displaystyle \sum _{j=1}^p}\underset{\alpha \downarrow 0+}{\text{lim}\text{inf}}\left\{\frac{\lambda {∥b_j^{\infty }+\alpha {\delta }_j∥}_2-{∥b_j^{\infty }∥}_2}{\alpha }\right\}\hfill \\ ={\displaystyle \sum _{j=1}^p}\left\{{\nabla }_{j}f{\left(b^{\infty }\right)}^T{\delta }_j+\partial P\left(b^{\infty };{\delta }_j\right)\right\}\ge 0.\hfill \end{array}$

Finally, we show that the sequence generated from the proposed algorithm satisfies these two conditions.

Theorem 3 Assuming the data (y, X) lies on a compact set and no column of X is identically 0, the sequence {b^(k)} generated from the algorithm converges to the minimum point of the objective function J (b).

Proof We only need to show that the generated sequence meets the conditions in Lemma 2.

Let b(u) be the vector containing u as its j th block of parameters with other blocks being the fixed values.

where $b^{\left(k\right)}\left(j\right)={\left[b_1^{{\left(k+1\right)}^T},\dots ,b_j^{{\left(k+1\right)}^T},b_{j+1}^{{\left(k\right)}^T},\dots ,b_p^{{\left(k\right)}^T}\right]}^T.$

Note that by Lemma 1, {J (b^(k))} converges as it keeps decreasing and is bounded from below. The convergence of {J (b^(k))} yield the convergence of {b^(k)}. Hence, conditions of Lemma 2 hold which imply that the limit of {b^(k)} is the minimum point of J (b).

Implementation

where $A_{\lambda }^{*\left(b\right)}\left(i,j\right)$ is the (i, j)'s entry of $A_{\lambda }^{*\left(b\right)}$ and #{·} is the number of elements in that set.

The final network topology can be obtained by setting a threshold, edges with probabilities or scores less than which are considered nonexistent. This study only focus on giving a list of edges with scores. The selection of threshold is not discussed here. The stability selection procedure can also be applied with the group LASSO method (6).

Since the data used are time series data, the moving block bootstrap method is employed in the first step to draw bootstrap samples from each dataset. For a dataset with n observations, in the moving block bootstrap with block length l, the data is split into n − l + 1 blocks: block j consists of observations j to j + l − 1, j = 1, . . . , n − l + 1. [n/b] blocks are randomly drawn from n − l + 1 blocks with replacement and are aligned in the order they are picked to form a bootstrap sample.

Another tuning parameter δ controls the degree of robustness. Generally, it picks $\delta =1.345\widehat{\sigma }$ where $\widehat{\sigma }$ is the estimated standard deviation of the error and $\widehat{\sigma }=MAD/0.6745$, where MAD is the median absolute deviation of the residuals. In this study, we use the least absolute deviations (LAD) regression to obtain the residuals. To avoid the overfitting of LAD which leads to a very small δ, we choose by $\delta =\text{max}\left(1.345\widehat{\sigma },1\right)$.