There are two principally different ways to construct co-expression networks from microarray data. One approach follows a hard- and the other a soft-thresholding of correlation coefficients (Zhou et al., 2002; Horvath and Dong, 2008). In the following r_ij = |cor(i,j)| is the absolute value of the Pearson correlation coefficient. In the case of hard-thresholding, the correlation coefficient r_ij between gene i and j is assessed with respect to a threshold τ and both genes are connected by an edge, A_ij = A_ji = 1, if r_ij ≥ τ, otherwise the genes remain unconnected, A_ij = A_ji = 0. The threshold or cut-off parameter τ can be obtained from randomization of the data allowing the assessment of statistical significance (Carter et al., 2004). This results in an undirected, unweighted network.

For the soft-thresholding two types of adjacency functions are frequently used (Zhang and Horvath, 2005). The sigmoid function,

and the power adjacency function,

Both types of adjacency functions lead to undirected but weighted networks. In order to choose the above parameters appropriately Zhang and Horvath (2005) suggested to adjust them in a way that the resulting network has approximately a scale-free degree distribution.

We would like to remark that the purpose for the construction of co-expression networks is different to all other methods discussed in this paper. Co-expression networks serve as means to explore the functionality of genes on a systems level (Zhang and Horvath, 2005) and do not aim to be causal representations of regulatory networks. Nevertheless, we included them in this review to point out that also networks that are not causal can be very useful in gaining biological understanding and insights.

Asymmetric-N is an algorithm that takes the fact into account that biological networks contain hubs. It is a modified version of Symmetric-N (Agrawal, 2002) which was designed for the construction of co-expression networks (Chen et al., 2008a). The two-step algorithm of Agrawal (2002), Symmetric-N, utilizes the N-nearest-neighbor concept. In a first step, for each node in the network its correlation to other nodes is calculated and sorted in descending order. In a second step, each pair of nodes is evaluated one by one and if for a pair of nodes both are in the corresponding N-nearest-neighbors list, a connection between them is included. Otherwise, they are not connected. Here N is a pre-defined cut-off value for the number of neighbor nodes to be considered (Agrawal, 2002). It is demonstrated in Agrawal (2002) that the inferred network has a scale-free topology.

Instead of using N neighbors for all nodes as potential connections, a higher number N_C is assigned to some core nodes (e.g., TFs) and a smaller number N_P is assigned to periphery nodes (e.g., non-TFs). Since the numbers of potential neighbors are different for core and peripheral nodes, the algorithm was called Asymmetric-N (Chen et al., 2008a). Asymmetric-N employs the intuition that transcription factors (TFs) are more frequently connected to other genes than regulated genes. An interesting result found in Chen et al. (2008a) is that Asymmetric-N performed better than Bayesian networks for small sample sizes.

Xiong et al. (2004) use a structural equation model (SEM; Jöreskog, 1973; Bollen, 1989) to reconstruct gene networks. Structural equation models combine path analysis (Wright, 1921, 1934) and confirmatory factor analysis (Jöreskog, 1969; Brown, 2006).

Xiong et al.’s (2004) approach starts by assuming that the expression level of genes, X, can be described linearly by

Here X is a vector of p endogenous and ξ a vector of q exogenous variables, and B and Γ are p × p respectively p × q coupling matrices. The last term in equation (4), ε, represents noise. The structure of the regulatory gene network is coded by the entries of the coupling matrices B and Γ. The structural equation in equation (4) is learned in two steps. First, assuming the structure of the network is known the coupling parameters are estimated via maximum likelihood optimizing the distance between the (p + q) × (p + q) covariance matrix Σ(Θ) and the sample covariance matrix S estimated from the data. Here, Θ = f (B, Γ, Φ, Ψ), is a function of the parameters of the SEM and the covariances of ξ (Φ) and ε (Ψ) that can be expressed analytically. The second step consists in finding the network structure (of B and Γ) by an optimization method whereas the different models are assessed by Akaike’s information criterion (AIC).

Partial correlations of low-order have been employed in de la Fuente et al. (2004), Magwene and Kim (2004), Wille and Bühlmann (2006). This approach is capable to infer undirected networks only, because a correlation, of any order, is symmetric in its arguments. The principle working mechanism of this approach is as follows. Starting with a fully connected network, interactions between genes are iteratively excluded by moving toward an increasing order of the correlation coefficient. More specifically, starting from a fully connected adjacency matrix A, one calculates in the first step all pair-wise correlation coefficients, ρ_ij = ρ(X_i, X_j). If ρ_ij = 0 the connection between gene X_i and X_j is deleted, A_ij = A_ji = 0. Inthe second step all partial correlation coefficients of first-order, ρij|k=ρXi,Xj|Xk, are calculated for all triplets of genes with ρ_ij ≠ 0. Now, if there is at least one gene X_k for which ρ_ij|k = 0 holds the connection between gene X_i and X_j is deleted, A_ij = A_ji = 0. Continuation to higher orders is analogously. The exclusion of interactions is based on testing the null hypothesis ρ_ij = 0 against the alternative hypothesis ρ_ij ≠ 0 for all orders of the (partial) correlation coefficients. In (de la Fuente et al., 2004) this approach was applied up to order three. The resulting network from this procedure is frequently called an undirected dependency graph (UDG; Shipley, 2000). An implementation of this approach is called ParCorA (de la Fuente et al., 2004).

GGM, also known as covariance selection model, concentration graph, or Markov random field (Dempster, 1972; Whittaker, 1990; Koller and Friedman, 2009), is a graphical model which assumes that all variables are distributed according to a multivariate normal distribution with a specific structure of the inverse of the covariance matrix, Ω = Σ⁻¹, also called precision or concentration matrix. Network inference methods based on GGM make use of the relation,

connecting the partial correlation of full-order with the elements of Ω, w_ij ∈ Ω. The partial correlation is of full-order (with respect to the number of genes) because V\{ij} is the set of all genes excluding i and j, i.e., the largest possible set of genes not considering i and j. In case of a vanishing partial correlation in equation (5) one can write

From this conditional independence relation, the principle way to infer a network structure from GGM becomes apparent estimating a network in the following way. If ρ_ij|V\{ij} ≠ 0, according to a hypothesis test, we include an edge, A_ij = A_ji = 1, otherwise there is no edge between i and j.

Several approaches have been made to infer gene regulatory networks based on GGM (Wille et al., 2004; Schäfer and Strimmer, 2005; Li and Gui, 2006). These methods differ in the way the inverse of the covariance matrix, Σ⁻¹, is estimated and in the statistical tests employed to define significance. The reason for these technical variants comes from a variety of problems. First, if the number of samples is smaller than the number of genes, which is typically the case for genomics data, the sample covariance matrix is not positive definite and, hence, not invertible. Another problem is caused by the small samples size problem (Schäfer and Strimmer, 2005).

Bayesian networks allow identifying a DAG structure of a network. BN were among the first methods that have been applied to expression data to infer GRN (Friedman et al., 2000; Hartemink et al., 2001). In order to overcome the problem that only acyclic networks can be inferred by BN, which would be a severe limitation considering the fact that real biological networks contain many feedback loops and are, hence, cyclic, a modified method called dynamic Bayesian network (DBN) is used instead. Briefly, DBN unfold cyclic processes by mapping them onto a sequence of acyclic events, which can then be analyzed with a BN (Dean and Kanazawa, 1990; Koller and Friedman, 2009). For the inference of GRN, this approach has been applied in Husmeier (2003), Perrin et al. (2003), Zou and Conzen (2005). The variations in these approaches provide different methods for the structure learning problem, which is the integral part of a BN learned from data. Principally, a major practical problem of a BN is that the computational complexity of algorithms to learn their structure has been shown to be NP-hard for score-based approaches (Chickering et al., 1995; Chickering, 1996). Hence, BN or DBN can be either only applied to quite small networks or heuristic approximation methods have to be used that separate the score into tractable units that can be optimized by locally constraint search techniques making the computational complexity manageable (Friedman et al., 1999, 2000).

In Figure 6, we categorize BN as correlation-based methods because the simplest statistical realization of a BN is obtained by estimating (X ╨ Y|S)_P, see section 1, by means of partial correlation coefficients (Spirtes et al., 1993; Pearl, 2000). However, also non-linear approaches are possible.

RN (relevance networks). The principle idea of RN (Butte and Kohane, 2000) is to compute all mutual information (MI) values for all pairs of genes, for a given data set, and declare mutual information values as relevant if their corresponding value is larger than a given threshold I₀. The resulting network is constructed based on this threshold by including an edge between two genes in the respective adjacency matrix of the network, A_ij = A_ji = 1, if I_ij > I₀, otherwise no edge is included between i and j. In Butte and Kohane (2000) the threshold I₀ was found by randomization of the expression data set. From this randomization, mutual information values were re-calculated from which a reference distribution of mutual information values, resembling a null-distribution, was obtained. Based on this reference distribution the threshold I₀ was obtained by heuristic arguments. For this reason, mutual information values that are larger than I₀ are called relevant but cannot necessarily be called statistically significant.

ARACNE (algorithm for the reconstruction of accurate cellular networks). ARACNE (Basso et al., 2005; Margolin et al., 2006) is similar to RN because it also uses mutual information and a randomization of the data to obtain a threshold I₀ allowing to declare mutual information values significant if I_ij > I₀¹. If I_ij is found to be significant, than an edge is included in the corresponding adjacency matrix between gene i and j, A_ij = A_ji = 1, otherwise no edge is included. However, in contrast to RN, ARACNE performs a second step testing all gene-triplets (three genes with mutual information values larger than I₀) such that, for each triplet (ijk), the edge corresponding to the lowest mutual information value I1=Ii′j′, with (i′j′) = argmin{I_ij, I_jk, I_ik}, is eliminated from the adjacency matrix, if it is smaller than the second smallest MI value I₂ multiplied by a factor, i.e.,

Here 0 ≤ ∈ ≤ 1. The introduction of this step has been motivated by the so called data processing inequality (DPI; Cover and Thomas, 1991). The DPI is a relation between mutual information values which means loosely that a post-processing of data cannot increase its information content. Specifically, one can show (Cover and Thomas, 1991) that the DPI for the following relation between the three random variables,

implies that I(X,Z) ≤ I(X,Y). Due to the fact that the criteria in equation (8) is for ∈ > 0 less stringent than the DPI [equation (9)], ∈ is called tolerance parameter.

To ensure that the application of equation (8) results in an unique solution, independent of the order the triples have been selected, the procedure starts by listing all possible triplets in the network G. Then all of these triplets are tested. Hence, the results of these tests have no influence on subsequent tests of triplets. At the end of the second step, the resulting network represents the final result.

ARACNE employs two parameters, I₀ and ∈. The cut-off parameter I₀ is determined by a resampling method estimating the distribution of the null hypothesis corresponding to a vanishing mutual information. This allows to assign p-values to mutual information values. In contrast, the DPI is not directly connected to statistical inference, but serves as a filtering step. Optimal values of ∈ are found from simulation studies that allow a comparison with the underlying true network. Hence, I₀ is found in an unsupervised and ∈ in a supervised way of learning.

CLR (context likelihood of relatedness). CLR (Faith et al., 2007) is also an extension of RN. It starts by estimating the pair-wise mutual information values for all genes. Then, CLR estimates the statistical likelihood of each MI value I_ij, for a particular pair of genes (ij), by comparing this MI value to a background. Specifically, for each gene pair (ij) two z-scores are obtained, one for gene i and one for gene j, by comparing the mutual information value I_ij with gene specific distributions, p_i and p_j. Here the two distributions p_i and p_j correspond to the distributions of mutual information values related to gene i ({I_ik|k ∈ V}) and gene j ({I_jk|k ∈ V}). By making a normality assumption about these distributions, corresponding z-scores, z_i and z_j, can be obtained from which the joint likelihood measure zij¯=zi2+zj2 is calculated. In contrast to RN and ARACNE, which employ a global threshold I₀ for each mutual information value correspondingly pair of genes, CLR estimates individual thresholds by considering an individual background for each pair of genes. This procedure can be seen as a first attempt to take the network context of gene pairs into account.

C3NET (conservative causal core). C3NET consists of two main steps (Altay and Emmert-Streib, 2010a). The first step is for the elimination of non-significant edges, whereas the second step selects for each gene the edge among the remaining ones with maximum mutual in formation value. The first step is similar to RN (Butte and Kohane, 2000), ARACNE (Margolin et al., 2006), or CLR (Faith et al., 2007) essential for eliminating non-significant links, according to a chosen significance level α, between gene pairs. In the second step, the most significant link for each gene is selected. This link corresponds also to the highest MI value among the neighbor edges for each gene. This implies that the highest possible number of edges that can be inferred by C3NET is equal to the number of genes under consideration. This number can decrease for several reasons. For example, when two genes have the same edge with maximum MI value. In this case, the same edge would be chosen by both genes to be included in the network. However, if an edge is already present another inclusion does not lead to an additional edge. Another case corresponds to the situation when a gene does not have significant edges at all. In this case, apparently, no edge can be included in the network. Since C3NET employs MI values as test statistics among genes, there is no directional information that can be inferred thereof. Hence, the resulting network is undirected and unweighted. For a detailed explanation of C3NET and technical details, the reader is referred to Altay and Emmert-Streib (2010a, 2011), de Matos Simoes and Emmert-Streib (2011).

SA-CLR (synergy augmented CLR). As indicated by its name, SA-CLR is based on CLR but includes synergistic effects between genes (Anastassiou, 2007; Watkinson et al., 2009). Synergy,

is defined as the difference between the generalized two-way mutual information,

and the mutual information values for two random variables. A visualization of synergy is shown in Figure 4. In this figure, the different intersections of the three variables X, Y, and Z are labeled by lower case letters, whereas “e” corresponds to the synergy, Syn(X,Y;Z). Using the definitions of the occurring mutual information values one obtains a simplified form of the synergy, given by

which is the negative of the three-dimensional interaction information (Watkinson et al., 2009). We want to remark that I(X; Y; Z) can assume positive as well as negative values (McGill, 1954). For this reason the intersection “e” in Figure 4 is shown in blue and red to emphasize this. Due to the fact that a Venn diagram visualizes logical relations between variables Figure 4 is no Venn diagram in the strict sense but a generalization thereof, allowing also to express negative values by including colors in the representation. The idea of SA-CLR consists in utilizing the synergy given in equation (10) in the following way,

This term is called synergistic regulation index (SRI; Watkinson et al., 2009). The interpretation of SRI is that for two given genes, X and Y, we are searching a third one, Z, that maximizes the synergy Syn(X, Y; Z) obeying the above constraints. This gene, Z, can be be arbitrarily chosen among all available genes, besides X and Y. The two constraints that involve mutual information, are chosen assuming that gene X and Z regulate Y. In this case, the mutual information between the two regulators, X and Z, should be smallest, motivating the constraints. Interestingly, these constraints introduce an asymmetry in X and Y² making it possible to speak about regulator and modulator genes. In the case of a positive synergy, a significant high value of SRI may indicate that gene X regulates genes Y (Watkinson et al., 2009).

In order to detect also non-cooperative effects, the measure actually used by SA-CLR is I(X, Y) + S(X, Y), the mutual information between two genes plus their synergistic regulation index. The significance of I(X, Y) + S(X, Y) is detected in a similar way as described above for CLR.

It is interesting to note that

holds formally but is not identical to the (general) conditional mutual information I(X; Y|Z) because the synergy for X, Y, and Z has to be positive. To make this clear, we included this constraint visibly in equation (20). This allows SA-CLR, in contrast to CLR, to obtain a directed network by assigning an edge from gene X to Y if I(X; Y) + S(X, Y) is significant.

MRNET (maximum relevance, minimum redundance). MRNET (Meyer et al., 2007) is an iterative algorithm that identifies potential interaction partners of a target gene Y that maximize a scoring function. Its working mechanism is as follows,

Whenever a gene, X_j, is found with a score that maximizes equation (21) and s_j is above a threshold, s₀, then this gene is added to the set S. The algorithm is iterated until no further gene can be found that would pass the threshold test. The basic idea of MRNET is to find interaction partners for Y that are of maximal relevance [first term in equation (22)] for Y, but introduce a minimum redundancy [second term in equation (22)] with respect to the already found interaction partners in the set S. Starting with a fully connected, undirected network among all genes, MRNET reduces successively edges between Y and VS, which have not maximized the score in equation 22.

We want to remark that the score s_j can be approximated by the difference between two mutual information values,

where the (auxiliary) random variable, A, represents the influence of the variables in the set S. Equation (23) is visualized in the left Figure 5. Here, I(X_j; Y) corresponds to “d” and I(X_j; A) to “b.” Due to the fact that the mutual information is always positive, I(X_j; Y) adds a positive and I(X_j; A) a negative value to s_j, indicated by blue and red colors. It has been argued in (Meyer et al., 2007) that due to the maximization of equation (22), MRNET approximates the conditional mutual information I(X_j;Y|S). This is also motivated by Figure 5.

MI3 (mutual information 3). The MI3 algorithm (Luo et al., 2008) uses three-way mutual information for the inference, hypothesizing that gene regulation commonly involves more than one regulator genes. The value of the measure MI₃ is defined as,

which equals the difference of mutual informations between the target gene and the two regulators and the target gene with one of the regulators. Alternatively, MI₃ can also be written as

Equation (25) is the sum of two conditional mutual informations between the target gene and a regulator given the other regulator, which emphasizes the three-way nature of this measure. Interestingly, MI₃ has also a simple relation to synergy, as shown in equation (26). This follows directly from equations (10 and 24).

In the right Figure 5, we show a visualization of MI₃. Here, the shown intersections result from the correspondence of I(T; R₁, R₂) to “d + e + f,” I(T; R₁) to “d + e,” and I(T; R₂) to “e + f.” Subtraction of these terms according to equation (24) results in “d + f.” Note that I(T; R₁, R₂) is counted twice.

The MI3 algorithm learns gene regulatory networks in two steps: first, local regulatory networks consisting of only three genes (T, R₁, and R₂) are learned. Starting from a given gene, T, two regulator genes are searched which maximize MI3 (T, R₁, R₂). As a result, directed edges between R₁ and T and R₂ and T are added constituting a small regulatory network. Second, the local regulatory networks learned in step one are assembled. To ensure that the resulting network is acyclic, there may be a need to remove edges forming cycles in the assembled network. This is solved heuristically, by identifying all local three-gene networks that contribute to a cycle and elimination of all edges of them from the overall network, except from the one with the highest MI₃ value. Overall, the MI3 algorithm aims to learn the optimal two-parent causal model for each target variable in the form R1 → T and R2 → T.

CMI (conditional mutual information). Soranzo et al. (2007) use only the conditional mutual information (CMI) to infer regulatory networks. They estimate all conditional mutual information values I(X_i; X_j | X_k) between triplets of genes. Starting from a fully connected network, they successively remove edges A_ij = A_ji = 0 if I(X_i; X_j | X_k) = 0 for at least one gene X_k, according to a threshold condition.

MI-CMI. The MI-CMI algorithm (Liang and Wang, 2008) uses both mutual information and conditional mutual information to estimate the network. Starting from a fully unconnected network G, the algorithm consists of three steps. First, all pairs of mutual information values I(X_i; X_j) are estimated. If I(X_i; X_j) ≥ I₀, then an edge is included, A_ij = A_ji = 1. Second, for all triplets (X_i, X_j, X_k) with I(X_i; X_j) ≥ I₀, I(X_i; X_k) ≥ I₀, and I(X_j; X_k) ≥ I₀ the conditional mutual information for all combinations of I(X_i, X_j, X_k) is estimated. Then based on heuristic rules, is it decided which edges in G are likely to result from indirect regulations (interactions). These edges will be removed from G. Third, for all triplets (X_i, X_j, X_k) with I(X_i; X_j) < I₀, I(X_i; X_k) < I₀, and I(X_j; X_k) < I₀ the conditional mutual information for all combinations of (X_i; X_j; X_k) is estimated. Again, following some heuristic rules, edges that are likely the result from interactive regulations, which could not be detected by mutual information, are included in G. Special emphasize was given to the used statistical estimators, for the mutual information and conditional mutual information values.

Mutual information in combination with conditional mutual information was also used in Zhao et al. (2008), Zhang et al. (2011). However, in contrast to (Liang and Wang, 2008), these approaches are conceptually closer to the pc algorithm (Spirtes et al., 1993; Shipley, 2000).