Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

This article is part of the supplement: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2011: Bioinformatics

Open Access Proceedings

Inference of gene regulatory subnetworks from time course gene expression data

Xi-Jun Liang1, Zhonghang Xia2*, Li-Wei Zhang1 and Fang-Xiang Wu3

Author Affiliations

1 School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China

2 Dept. of Mathematics and Computer Science, Western Kentucky University, Bowling Green, KY 42101, USA

3 Department of Mechanical Engineering, University of Saskatchewan, 57 Campus Dr., Saskatoon, SK S7N 5A9, Canada

For all author emails, please log on.

BMC Bioinformatics 2012, 13(Suppl 9):S3  doi:10.1186/1471-2105-13-S9-S3


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/13/S9/S3


Published:11 June 2012

© 2012 Liang et al; licensee 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.

Abstract

Background

Identifying gene regulatory network (GRN) from time course gene expression data has attracted more and more attentions. Due to the computational complexity, most approaches for GRN reconstruction are limited on a small number of genes and low connectivity of the underlying networks. These approaches can only identify a single network for a given set of genes. However, for a large-scale gene network, there might exist multiple potential sub-networks, in which genes are only functionally related to others in the sub-networks.

Results

We propose the network and community identification (NCI) method for identifying multiple subnetworks from gene expression data by incorporating community structure information into GRN inference. The proposed algorithm iteratively solves two optimization problems, and can promisingly be applied to large-scale GRNs. Furthermore, we present the efficient Block PCA method for searching communities in GRNs.

Conclusions

The NCI method is effective in identifying multiple subnetworks in a large-scale GRN. With the splitting algorithm, the Block PCA method shows a promosing attempt for exploring communities in a large-scale GRN.

Background

Rapid advances in high-throughput DNA microarray technology generate a huge amount of time course gene expression data which, in turn, calls for efficient computational models to characterize the network of genetic regulatory interactions. A number of methods have been proposed to infer GRNs from gene expression data. Boolean networks [1] use two states, "ON" or "OFF" to represent the state of each gene, and each state at the next time step is determined by Boolean logical rules. Bayesian Networks [2] infer causal relationships between two genes according to conditional probability functions. The stochastic nature makes them more accurate in modeling the dynamics and nonlinearity of gene regulation in large-scale systems. Bayesian Networks, however, usually do not include cycles and, thus, are difficult to deal with feedback motifs. Ordinary differential equations (ODEs) models [3-5] overcome this problem by modeling GRNs as a set of differential equations. Some other models such as signed directed graphs, multiple regression, state space model, etc., are addressed in the survey [6].

Whereas most of the existing work focuses on small-sized GRNs, limited attention has been given to interactions among large scale genes. Conventional approaches are usually designed for the network with connectivity less than a small fixed number [7]. Computational complexity is a major obstacle in reconstruction of large scale GRNs as determining the parameters in such a network is time-consuming. Sparsity is a common assumption used in modeling GRNs to reduce the computational complexity. Typically, in a sparse network, one gene interacts with only a couple of genes [7].

Recently, Yuan et al. [8] proposed a directed partial correlation (DPC) method for regulatory network inference on large-scale gene data. The DPC method combines the directed network inference approach and Granger causality concept for causal inference on time series data to reconstruct large-scale GRNs. Although modular discovery was provided by biclustering in gene expression data, the DPC method cannot present multiple sub-networks simultaneously.

We propose the NCI method for subnetwork identification by detecting community structures from large-scale gene expression data. Usually, GRNs have community structures: genes in the same groups are found with high density of "within-group" interactions and genes in different groups with low density of "between-group" interactions [9]. Many algorithms have been proposed to detect community structures by clustering [9-15]. To accomodate the large-scale GRN inference, we particularly propose a block principal component analysis (Block PCA) method, which explores community structure information for the NCI method.

The NCI method repeats two steps: (1) N-step: identify possible gene regulatory networks; (2) C-step: estimate community structure. At the N-step, a convex quadratic programming, formulated for the community structure, is solved to infer possible GRNs. This quadratic programming can be identically divided into n (the total number of genes) sub-problems, each of which has a much smaller dimension, and, thus, adapt to large-scale networks. At the C-step, the NCI method estimates community structure from the GRNs identified at the first step. When the algorithm terminates, a network with community structures is obtained.

Methods

An ODE model for GRNs

The processes of transcription and translation in a GRN consisting of n genes can be modeled as the following dynamic system:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M1">View MathML</a>

(1)

where vector x = [x1, x2, ..., xn]T ∈ ℝn is the concentration of mRNAs of n genes, C = diag [-c1,-c2, ..., - cn] ∈ Rn × n, ci represents the degradation rate of gene i, r = [r1, r2, ..., rn]T ∈ ℝn is the reaction rates which is a function of concentrations of some mRNAs, and matrix S Rn × n represents the stoichiometric matrix of the biological network. For simplicity, one can assume the reaction rate r is a linear combination of mRNAs concentrations. Let FRn × n be the coefficient matrix. Then,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M2">View MathML</a>

(2)

By substituting (2) into (1), we have

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M3">View MathML</a>

(3)

A standard discretization of system (3) by using the zero-order hold method on m observation points for a given sampling time Δt is

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M4">View MathML</a>

(4)

where A = eCΔt + (eCΔt - I)C-1SF.

Let X ∈ ℝn × m be a matrix of gene expression data, with the columns being the measured gene expression levels at m time points, and n being the number of genes. Let X1 and X2 be the sub-matrices of X made up by the first m - 1 columns and last m - 1 columns of X, respectively. According to [16], the gene regulatory network can be inferred by solving the following optimization problem:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M5">View MathML</a>

(5)

where ||·||2 is the Euclidean norm. Stability is usually used as a criterion to determine the qualification of the inferred GRN. For discrete models, A = (aij)n × m is stable if

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M6">View MathML</a>

(6)

Moreover, since the network is commonly recognized as sparse, l1 regularization is added to Eq. (5) to obtain a sparse matrix A. Hence, with the sparsity and stability conditions, (5) becomes

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M7">View MathML</a>

(7)

where γ is a positive scalar, ||A||1 = ∑i, j |aij| is the l1-norm of matrix A.

The NCI method

Since rows of A are independent in the objective function and constraints, problem (7) can be divided into n sub-problems and solved individually [16]. However, such a solution does not consider the information of community structure which implies multiple subnetworks. In this section, we propose the NCI method to overcome this problem. An observation is that interactions between genes in a community occur more frequently than those between different communities. We introduce a weighted matrix W = (wij)n × m to distinguish genes in the communities with those outside. wij is assigned a small positive value or zero if gene i and j located in the same community; a relatively large value, otherwise.

By adding term 〈W, |A|〉 to (7), we have

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M8">View MathML</a>

(8)

Where μ > 0 is a penalty parameter, |A| = (|aij|)n × m, 〈W, |A|〉 = trace (WT |A|) = ∑i, j wij|aij|. All elements of matrix W are nonnegative.

We propose a clustering method, named Block PCA, to update weight matrix W. With Block PCA, we can obtain matrix L*, reflecting the community structure of its corresponding network. Then, weight matrix W can be updated by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M9">View MathML</a>

(9)

Where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M10">View MathML</a> is the matrix with all 1's.

For example, consider a network with five nodes and

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M11">View MathML</a>

Node 1, 2, and 3 form a community, and node 4 and 5 form another community. Particularly, we apply sparse singular value decomposition (SSVD) [17] on a general L* to identify the communities in GRNs. The NCI method is summarized in Algorithm 1.

Some additional details about Algorithm 1:

1. Stop criteria. The NCI algorithm stops when either of the following two criteria meet. (1) Weighted matrix W converges, that is, ||W(k) - W(k+1)|| ≤ tol for a pre-defined constant tol > 0, where W(k) denotes W at iteration k; (2) The number of iteration reaches the threshold.

2. The efficiency of the algorithm mainly depends on the estimation of the community structure of the underlying GRN. Since matrix A in (8) provides a base for estimation of community structure L* computed by (12), a poor estimation of A may result in an inaccurate W. Hence, instead of using only one estimation of A, we average out the errors by calculating a series of estimations with different arguments γ in (8) and combining them together. More specifically, we choose a set of γ1, ..., γq and compute the corresponding solutions A1, ..., Aq. Then, A in Step 1 of Algorithm 1 is set as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M12">View MathML</a>

After the iteration terminates, model (8) is solved again to compute the matrix A with γ = γτ, where γτ, is a parameter in Algorithm 1.

3. The complexity of the subproblem is our primary concern about the design of the NCI algorithm. Since the subproblems may be called iteratively in Algorithm 1, the complexity of the NCI algorithm is determined by those subproblems. Both sub-problems (8) and the Block PCA model are convex, and can be efficiently solved by CVX [18] and and the proposed splitting algorithm, respectively. As aforementioned, model (8) is dividable: rows of A in the objective function and the constraints are independent. Hence, it is equivalent to n sub-problems:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M13">View MathML</a>

(10)

for i = 1, ..., n, where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M14">View MathML</a> is the i-th row of the matrix <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M15">View MathML</a> is the i-th row of W and A, respectively. Sub-problem (10) can be transformed into a standard (convex) quadratic programming, and solved by software packages such as Mosek or CVX [18].

The Block PCA model

The Block PCA model is motivated by Robust PCA model [19]

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M16">View MathML</a>

(11)

which ||L||* is the nuclear norm of the matrix L, and D is a given matrix.

The block PCA aims to seek a block submatrix in D by solving optimization problem

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M17">View MathML</a>

(12)

where W1 is a weight matrix with all elements nonnegative.

In Block PCA, D ∈ ℝn, n is set to be matrix <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M10">View MathML</a> with all 1's, λ2 is constantly set as <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M18">View MathML</a> and λ1∈(0, λ2). For a network with n nodes, we define weight matrix W1 = (w1ij)n × n, where

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M19">View MathML</a>

(13)

pij is the length of the shortest path between the node i and j, p0 ≥ 0 is a parameter less than the diameter of the network.

As in Robust PCA (11), the nuclear norm || · ||* usually induces a low rank matrix and the l1 norm || · ||1 induces a sparse matrix [19,20]. The constraint D = L + E enforces to split matrix D into a low rank matrix L and a sparse matrix E. Different with Robust PCA, the Block PCA adds an extra term λ1 W1, |L|〉 = λ1 w1ij·|Lij| to (11). The nonnegative weight matrix W1 stands for the prior knowledge about low rank matrix L.

Splitting algorithm for solving Block PCA

Block PCA model (12) can be transformed to a linear semidefinite programming (SDP)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M20">View MathML</a>

(14)

However, this transformation increases the size of the variable matrix from n × n to 2n × 2n. Existing SDP solvers such as CVX [18] can not solve large-scale SDP problems. Instead, we solve Block PCA problem (12) by extending the splitting method [21] for optimization problem

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M21">View MathML</a>

(15)

Where θi: <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M22">View MathML</a> are closed convex functions, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M23">View MathML</a>, b∈ℝl.

Note that Block PCA (12) can be recast as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M24">View MathML</a>

(16)

By letting θ1(·): = ||·||*, θ2(·): = λ1W1, |·|〉, θ3(·): = λ2||·||1, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M25">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M26">View MathML</a> Block PCA (12) can be treated as a generalized case of (15) with matrix variables L, E, U and linear operators <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M27">View MathML</a>

Under the framework of [21], we next present an implementable splitting algorithm for the Block PCA model (12).

Define operator

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M28">View MathML</a>

(17)

t ∈ ℝ and ε > 0. It can be extended to an arbitrary matrix X ∈ ℝn, n by applying element-wise operation, denoted by <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M29">View MathML</a>.

Consider the sigular value decompostion (SVD) of the matrix X

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M30">View MathML</a>

(18)

where U and V are orthogonal matrices consisting of singular vectors, and Σ is the diagnal matrix made up of the singular values. For each τ >0, the soft-thresholding operator <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M54">View MathML</a> is defined as [22]

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M31">View MathML</a>

(19)

More generally, for a matrix W ∈ ℝn, n with all elements nonnegative, we define

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M32">View MathML</a>

(20)

Particularly, if W is the matrix with all elements 1, ||X||w degenerates to ||X||1, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M33">View MathML</a> degenerates to <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M29">View MathML</a>.

Let <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M34">View MathML</a> where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M35">View MathML</a>. Then, for the calculated (Lk, Ek, Uk, Λk), the steps for each iterative (Lk+1, Ek+1, Uk+1k+1) for solving (12) are as follows.

Step 1. Solve Lk+1 by the following problem.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M36">View MathML</a>

(21)

By Theorem 2.1 in [22], the unique solution of (21) is

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M37">View MathML</a>

(22)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M38">View MathML</a>

Step 2. Update the Lagrangian multiplier <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M39">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M40">View MathML</a>

(23)

Step 3. Solve Uk+1, Ek+1 by the following two problem.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M41">View MathML</a>

By the property of the operator Sτ [Y ] shown in [23],

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M42">View MathML</a>

(24)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M43">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M44">View MathML</a>

(25)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M45">View MathML</a>

Step 4. Update the Lagrange multiplier Λk+1 by Lk+1, Ek+1.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M46">View MathML</a>

(26)

The algorithm can be terminated when

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M47">View MathML</a>

(27)

and

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M48">View MathML</a>

(28)

for tolerance ε1 >0, ε2 >0, where ΔLk = Lk+1- Lk, ΔEk = Ek+1 - Ek, ΔUk = Uk+1- Uk.

The splitting algorithm for solving Block PCA model is summarized in Algorithm 2.

In Algorithm 2, arguments β and μ are currently set constant. Adaptive settings of these arguments may speed up the convergence. The discussion of this issue in a simple case can be referred to [24].

Results and discussion

We examine the NCI method based on two synthetic gene regulatory networks with different sizes. The GRN in first test is a small-sized network consisting of 14 genes and 27 interactions. There exist two communities in this GRN. In the second test, the network consists of 50 genes and 100 interactions and the data come from the Artificial Gene Network database [25]. Since the gene network is synthetic, the corresponding matrix A in (5) is known beforehand. We solve the GRN by the NCI method and compare it with A to evaluate the performance of the algorithm. Moreover, we examine the performance of the proposed splitting algorithm in the third test.

The metric used in the performance examination was introduced in [16]. It compares the signs of the estimated matrix Ae with A. The accuracy is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M52">View MathML</a>

(29)

where r11, r22, and r33 are the number of correctly identified positives, zeros and negatives, representing promotions, repressions, and no interaction, respectively.

The algorithm runs on a computer with Pentium (R) dual-core CPU E5200 2.50GHz, and RAM 2.0GB. The parameters of the algorithm are chosen as follows. In Test 1, γ in problem (8) is chosen from {0.05, 0.02, 0.008} to find possible GRNs and γτ = 0.02. In Test 2, γ is chosen from {0.02, 0.005, 0.001} and γτ = 0.005. In the first two tests, μ is chosen as 10γ for problem (8), λ1 as 0.2λ2 in the Block PCA model, and p0 as <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S9/S3/mathml/M53">View MathML</a> for Eq. (9), where d is the diameter of the corresponding network. The algorithm terminates in 3 iterations.

Test 1. A small gene network with 14 genes

Figure 1 shows the network and its two communities. The diameter of this gene network is 6. We choose different initial gene expression levels randomly for 30 times. The corresponding 30 accuracy rates of the calculated GRN are shown in Figure 2(A). "NCI" and "SGN" denote the NCI algorithm and sparse gene regulatory network method [16], respectively. Compared with the SGN algorithm, the NCI significantly improved prediction accuracy. In the noise case, 10% elements of the gene expression matrix X are incorporated with Gauss noise with zero mean and unit variance. The accuracy rates of two methods are shown in Figure 2(B). In both of the noise-free (Figure 2(A)) and noise cases (Figure 2(B)), the NCI method has much better performance in most of the 30 runs. In the noise-free case, the NCI algorithm increases the average accuracy from 78.8% to 83.5%. In the noise case, the NCI algorithm increases the average accuracy from 87.3% to 88.9%.

thumbnailFigure 1. Synthetic small gene network. This synthetic gene network consists of 14 genes and 27 interactions. There are two communities: gene 1-5 form a community, gene 6-10 form the other. "→" indicates promotion interaction, while "⊣" indicates repression.

thumbnailFigure 2. Accuracy of NCI with the 30 runs. The accuracies of the NCI method and SGN method are shown in (A), (B) with 30 random experiments. (C), (D) shows the efficiency of searching multiple possible GRNs at N-step of NCI method. "NCI" indicates normal NCI method searching multiple possible GRNs, "SGN" denotes the sparse gene regulatory network method [16], and "NCI-single" indicates the NCI method searching a single GRN at N-step. In (A) and (C), the expression levels are accurate, while in (B) and (D) 10% elements of the expression levels are incorporated with Gaussian noise with zero mean and unit variance.

To show the effectiveness of the NCI method at N-step in searching multiple possible GRNs, we compare the accuracy rates with the results of one iteration at N-step (γ = γτ at N-step). As shown in Figure 2(C), the average accuracy is improved from 78.5% to 83.5% in the 30 runs. In noise case (Figure 2(D)), the average accuracy is improved from 87.6% to 88.9%. Thus, a number of iterations at N-step is necessary for finding accurate GRNs with the NCI algorithm.

Test 2. A gene network with 50 genes

The network in the second test consists of 50 genes and 100 interactions (See Figure 3(A)). Network statistics are listed in Table 1. The nodes in red in Figure 3(A)) form a unique community. The inferred network by NCI algorithm contains 41 genes and 87 interactions. As shown in Figure 3(B), the community identified by the NCI algorithm has a very large overlap with the true community. Among 34 genes in the true community, 23 important ones (with large in-degree and out-degree) are successfully identified.

thumbnailFigure 3. Inference of the web100-023 gene regulatory network. The web100-023 gene network consists of 50 genes and 100 interactions which is shown in (A), where "→" indicates promotion interaction, while "┤" indicates repression. The 34 nodes in red in (A) form a unique community. With the 34 genes in the true community, the community identified by the NCI algorithm has an overlap of 23 important genes (with large in-degree and out-degree). The overlap is indicated by nodes in red in (B).

Table 1. Characteristics of web100-023 network

Test 3. The performance of the Block PCA and splitting algorithm

The following experiments are specially designed to test the efficiency of the Block PCA method and the performance of the splitting algorithm as well. We randomly generate three clusters with 30 points (See Figure 4(A)). Three clusters calculated by K-means are shown in Figure 4(B). Based on the distances between these points, matrix W1 is calculated by Eq. (13) with p0 = 0.684. The results of the splitting algorithm are shown in Figure 4(C). The corresponding three clusters calculated by SSVD are displayed with different colors in Figure 4(D). Among 30 data points, two points (the point 25 and 30) are outliers, not included in any cluster. The clustering result of the remaining 28 points is identical with that of K-means.

thumbnailFigure 4. Clusters produced by Block PCA and K-means. In (A), 30 points are randomly generated on a plane. Three clusters calculated by K-means are shown in (B). The low rank matrix calculated by the splitting algorithm is depicted in (C). The corresponding three clusters calculated by SSVD are displayed with different colors depicted in (D). Among 30 data points, two points (the point 25 and 30) are outliers, not included in any cluster. The clustering result of the remaining 28 points is identical with that of K-means.

To verify the effectiveness of the argument λ1, we choose different values and the calculated low rank matrices L are shown in Figure 5. It is shown that the number of nonzero elements of L (white pixels of the images) decreases, as λ1 increases. The numbers of nonzero elements of L are 804, 485, 284 and 100, with λ1 = 0.2λ2, 0.4λ2, 0.6λ2, and 0.7λ2 in Figure 5(A), (B), (C) and 5(D), respectively.

thumbnailFigure 5. Low rank matrices of the Block PCA model with various λ1. The figure depicts the low rank matrices L of the Block PCA model with different values of λ1. In (A), (B), (C) and (D), the value of λ1 of the Block PCA model is chosen as 0.2λ2, 0.4λ2, 0.6λ2, and 0.7λ2, respectively.

We compare the performance of the splitting algorithm with CVX and SDPNAL [26] by which the Block PCA model is solved via the SDP formulation (14). The results are listed in Table 2, where "Points30" indicates calculating on the data of 30 points on a plane, "funVal" indicates the calculated objective function value for the Block PCA model, "split", "cvx" and "sdpnal" indicate splitting method, CVX and SDPNAL, respectively. It is shown in Table 2 that splitting algorithm outperforms others in all the tests.

Table 2. Results of splitting algorithm, CVX and SDPNAL for solving the Block PCA model

Conclusion

We have developed the NCI method for gene regulatory network reconstruction from gene expression data. Based on the convex programming technology, the NCI method has shown the capability to identify multiple subnetworks within a large-scale gene regulatory network. The NCI method includes two main steps. At the first step, the algorithm infers a gene regulatory network. At the second step, the algorithm estimates potential community structures. These two steps repeat until the algorithm terminates. Furthermore, we have proposed an efficient Block PCA method for exploring communities within a GRN and the splitting algorithm for the Block PCA model. Numerical experiments have validated the effectiveness of the NCI method in identifying GRNs and inferring the communities.

Abbreviations

GRN: gene regulatory network; NCI: network and community identification; Block PCA: block principal component analysis.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

XL and ZX designed the NCI method and wrote the manuscript. XL and LZ designed the Block PCA method and splitting algorithm. XL and FW designed the ODE GRN model and experiments. All authors read and approved the final manuscript.

Acknowledgements

This article has been published as part of BMC Bioinformatics Volume 13 Supplement 9, 2012: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2011: Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/13/S9.

This research is partially supported by the Natural Science Foundation of China, Grant 11071029 and 11171049.

References

  1. Akutsu T, Miyano S, Kuhara S: Identification of genetic networks from a small number of gene expression patterns under the Boolean network model.

    Pac Symp Biocomput 1999, 17-28. PubMed Abstract OpenURL

  2. Bernard A, Hartemink A: Informative structure priors: joint learning of dynamic regulatory networks from multiple types of data.

    Pac Symp Biocomput 2005, 459-470. PubMed Abstract OpenURL

  3. Chen T, He H, Church G: Modeling gene expression with differential equations.

    Pac Symp Biocomput 1999, 29-40. PubMed Abstract OpenURL

  4. De Hoon M, Imoto S, Kobayashi K, Ogasawara N, Miyano S: Inferring gene regulatory networks from time-ordered gene expression data of Bacillus subtilis using differential equations.

    Pac Symp Biocomput 2003, 17-28. PubMed Abstract OpenURL

  5. Hu X, Ng M, Wu F, Sokhansanj B: Mining, modeling, and evaluation of subnetworks from large biomolecular networks and its comparison study.

    IEEE Trans Inf Technol Biomed 2009, 13(2):184-194. PubMed Abstract | Publisher Full Text OpenURL

  6. Huang Y, Tienda-Luna I, Wang Y: A survey of statistical models for reverse engineering gene regulatory networks.

    IEEE Signal Process Mag 2009, 26:76-97. PubMed Abstract | PubMed Central Full Text OpenURL

  7. Wu F: Inference of gene regulatory networks and its validation.

    Current Bioinformatics 2007, 2(2):139-144. Publisher Full Text OpenURL

  8. Yuan Y, Li C, Windram O: Directed partial correlation: inferring large-scale gene regulatory network through induced topology disruptions.

    PLoS One 2011, 6(4):e16835. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Newman M: Fast algorithm for detecting community structure in networks.

    Phys Rev E Stat Nonlin Soft Matter Phys 2004, 69(6):066133. PubMed Abstract | Publisher Full Text OpenURL

  10. Pothen A, Simon H, Liou K, et al.: Partitioning sparse matrices with eigenvectors of graphs.

    SIAM J Matrix Anal Applic 1990, 11(3):430-452. Publisher Full Text OpenURL

  11. Kernighan B, Lin S: An efficient heuristic procedure for partitioning graphs.

    Bell System Technical Journal 1970, 49(2):291-307. OpenURL

  12. Girvan M, Newman M: Community structure in social and biological networks.

    Proc Natl Acad Sci USA 2002, 99(12):7821-7826. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Radicchi F, Castellano C, Cecconi F, Loreto V, Parisi D: Defining and identifying communities in networks.

    Proc Natl Acad Sci USA 2004, 101(9):2658-2663. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Palla G, Derényi I, Farkas I, Vicsek T: Uncovering the overlapping community structure of complex networks in nature and society.

    Nature 2005, 435:814-818. PubMed Abstract | Publisher Full Text OpenURL

  15. Newman M: Detecting community structure in networks.

    The European Physical Journal B-Condensed Matter and Complex Systems 2004, 38(2):321-330. Publisher Full Text OpenURL

  16. Wu F, Liu L, Xia Z: Identification of gene regulatory networks from time course gene expression data.

    Conf Proc IEEE Eng Med Biol Soc 2010, 795-798. PubMed Abstract | Publisher Full Text OpenURL

  17. Lee M, Shen H, Huang J, Marron J: Biclustering via sparse singular value decomposition.

    Biometrics 2010, 66:1087-1095. PubMed Abstract | Publisher Full Text OpenURL

  18. Grant M, Boyd S: CVX: Matlab software for disciplined convex programming, version 1.21 (2010).

  19. Candés E, Li X, Ma Y, Wright J: Robust principal component analysis?

    ArXiv:0912.3599, in press. OpenURL

  20. Chandrasekaran V, Sanghavi S, Parrilo P, Willsky A: Rank-sparsity incoherence for matrix decomposition.

    SIAM J Optim 2011, 21(2):572-596. Publisher Full Text OpenURL

  21. He B, Tao M, Yuan X: A splitting method for separate convex programming with linking linear constraints.

    Tech rep 2010. OpenURL

  22. Cai J, Candés E, Shen Z: A singular value thresholding algorithm for matrix completion.

    SIAM J Optim 2010, 20(4):1956-1982. Publisher Full Text OpenURL

  23. Lin Z, Chen M, Wu L, Ma Y: The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices.

    ArXiv:1009.5055, in press. OpenURL

  24. Boyd S, Parikh N, Chu E, Peleato B, Eckstein J: Distributed optimization and statistical learning via the alternating direction method of multipliers.

    Foundations and Trends in Machine Learning 2010, 3:1-122. Publisher Full Text OpenURL

  25. AGN [http://www.comp-sys-bio.org/AGN/index.html] webcite

  26. Zhao X, Sun D, Toh K: A Newton-CG augmented Lagrangian method for semidefinite programming.

    SIAM J Optim 2010, 20(4):1737-1765. Publisher Full Text OpenURL

  27. web100-023 [http://www.comp-sys-bio.org/AGN/Web/Web100-023.html] webcite

  28. SDPNAL [http://www.math.nus.edu.sg/~matsundf/] webcite

  29. CVX [http://cvxr.com/cvx/] webcite