Mathematical modeling of the structural a priori

We first introduce our notations before detailing our structural models and variational formulation. Let G represent the total number of genes for which expression data is collected. Expression data is gathered in a symmetric weighted adjacency matrix \(\mathbf {W} \in \mathbb {R}^{G\times G}\). Its (i,j) element corresponds to a statistical measure reflecting the strength of the link, or information shared, between the expression profiles of gene i∈{1,…,G} and gene j∈{1,…,G}. Our approach uses non-negative weights. A convenient choice for ω _i,j is the normalized mutual information (ω _i,j ∈ [ 0,1]) computed between the expression profiles of genes i and j.

Let us comment the first term in the above functional. In order to select edges of strong weights ω _i,j , the first term reflects a biological data fidelity term. It represents a gene-to-gene edge deletion cost. Thus, if ω _i,j is large (respectively, small), its edge deletion cost is high (respectively, low), disfavoring (respectively, favoring) its deletion. We now explore the two last penalty terms of (2) corresponding to our biologically-related structural a priori regularization.

The second term counterbalances the first one. Independently from the fact that actual TF genes are less numerous than \({\overline {\text {TF}}}\) genes, regulatory relationships between couples of TFs are expected to be less frequent than between one TF and one \({\overline {\text {TF}}}\). This expectation may promote biological graphs with a modular structure [19, 20]. An illustration is presented in Fig. 1. As we are looking for gene regulatory knowledge, we infer edges linked to at least one TF. In addition, we want to favor the preservation of \({\text {TF}}\textrm {-}{\overline {\text {TF}}}\) edges over TF-TF links. This edge selection capability is driven by positive weights λ _i,j . Their values depend on the three types of pairs of nodes i and j. We define these case-dependent weights as follows:

$$ \lambda_{i,j} = \left\{ \begin{array}{ll} 2\,\eta & \text{if}~ i \notin \mathcal{T} \text{and}~ j \notin \mathcal{T}, \\ 2\, \lambda_{\text{TF}} & \text{if}~ i \in \mathcal{T} \text{and}~ j \in \mathcal{T}, \\ \lambda_{\text{TF}} + \lambda_{\overline{\text{TF}}} & \text{otherwise.} \end{array} \right. $$

(3)

Hence, \({\overline {\text {TF}}}\textrm {-}{\overline {\text {TF}}}\) edges have weights assigned to 2η. The parameter λ _TF (respectively, \({\lambda _{\overline {\text {TF}}}}\)) acts in the neighborhood of TF genes (respectively, \({\overline {\text {TF}}}\) genes). They may be interpreted as two threshold parameters. This double threshold promotes grouping between strong and weaker edges among functionally-related genes. A similar approach is used in image segmentation [21] to enhances object detection with reduced sensitivity to irrelevant features [22]. To promote \({\text {TF}}\textrm {-}{\overline {\text {TF}}}\) interactions, the λ _TF parameter should be greater than \({\lambda _{\overline {\text {TF}}}}\). To ensure that any TF involved interaction is selected first, we should verify that \(\eta \geq {\lambda _{\text {TF}}} \geq {\lambda _{\overline {\text {TF}}}}\). Additionally, removing all \({\overline {\text {TF}}}\textrm {-}{\overline {\text {TF}}}\) edges amounts to setting their corresponding x _i,j to zero. Consequently, η should exceed the maximum value of the weights ω. Since we address different data types and input weight distributions, we can easily renormalize them all to ω _i,j ∈ [ 0,1], and choose η=1. When \({\lambda _{\text {TF}}}={\lambda _{\overline {\text {TF}}}}\), no distinction is made on the type of edges. This is equivalent to using a unique threshold value, as in classical gene network thresholding. This can be interpreted as if, without further a priori, all genes were indistinguishable from putative TFs. However, different λ _TF and \({\lambda _{\overline {\text {TF}}}}\) may be beneficial. We indeed show in the Additional file 1 that for any fixed value of λ _TF, smaller values for \({\lambda _{\overline {\text {TF}}}}\) improve graph inference results. A simple linear dependence \({\lambda _{\text {TF}}} = \beta {\lambda _{\overline {\text {TF}}}}\), with β≥1 suffices to define a generalized inference formulation encompassing the classical formulation. We fixed here β as a parameter based on the gene/TF cardinal ratio: \(\beta = \frac {|\mathcal {V}|}{|\mathcal {T}|}\). This choice is consistent when no a priori is formulated on the TFs (i.e. all genes are considered as putative TFs). Hence, β=1 and \({\lambda _{\text {TF}}}={\lambda _{\overline {\text {TF}}}}\). As mentioned above, without knowledge on TFs, we recover classical gene network thresholding. The λ _i,j parameter now only depends on a single free parameter \({\lambda _{\overline {\text {TF}}}}\), similarly to the large majority of inference methods requiring a final thresholding step on their weights.

Finally, the third term of the proposed functional aims to enforce a regulator coupling property (see Fig. 2). If two transcription factors are co-expressed, and co-regulate at least one gene, we consider plausible that any gene regulated (respectively non regulated) by one of these TFs is regulated (respectively, non regulated) by the other TF. We quantitatively translate the co-expression of two TFs j and j ^′ by \(\phantom {\dot {i}\!}\omega _{j,j^{\prime }}>\gamma \), where \(\gamma \in \mathbb {R}^{+}\) is a threshold reflecting the strength of the co-expression between j and j ^′. Similarly, the regulation of a \({\overline {\text {TF}}}\) k by a TF j (respectively, j ^′) is numerically expressed by ω _j,k >γ (respectively, \(\phantom {\dot {i}\!}\omega _{j^{\prime },k}>\gamma \)). We define γ from robust statistics [23] as the (G−1)^th quantile of the weights. We thus choose the coupling parameter as:

$$ \rho_{i,j,j'} = \mu \frac{\sum\limits_{k \in {\mathcal{V}\setminus (\mathcal{T} \cup\{i\})}} 1(\min \{\omega_{j,j'}, \omega_{j,k}, \omega_{j',k}\} > \gamma)}{|\mathcal{V} \setminus \mathcal{T}|-1}, $$

where 𝟙 is the characteristic function (equals to 1 when the condition in argument is satisfied and 0 otherwise) and μ≥0 is a regularization parameter controlling the impact of the third term on the global cost. The proposed numerator counts the number of \({\overline {\text {TF}}}\) genes co-regulated by j and j ^′. As we exclude the gene i, the maximal number of \({\overline {\text {TF}}}\) genes co-regulated by j and j ^′ equals \(|\mathcal {V} \backslash \mathcal {T}|-1\). Hence, using the latter quantity as the denominator casts the \(\phantom {\dot {i}\!}\rho _{i,j,j^{\prime }} / \mu \) parameter as a co-regulation probability relative to couples of TFs (j,j ^′). The greater the co-regulation probability, the stronger the influence of the third term. This penalty requires that at least two \({\overline {\text {TF}}}\) target genes exist (hence, the denominator does not vanish). Otherwise, when \(|\mathcal {T}| = |\mathcal {V}|\), we set μ=0.

We now turn our attention to the strategy for computing an optimal labeling vector x ^∗ solution to Problem (2).

Optimization strategy

By using elements from graph theory, we now explain how a maximum flow algorithm can solve Problem (2). It relies on the maximum flow/minimum cut duality [24]: the computation of an optimal edge labeling minimizing (2) can be performed by maximizing a flow in a network G _f .

A flow (or transportation) network \(\mathcal {G}_{f}\) is a directed, weighted graph including two specific nodes, called source (a node with 0-in degree) and sink (a node with 0-out degree), respectively denoted by s and t. We recall that the degree of a node is defined as the number of edges incident to that node.

We now introduce the concept of flow in the transportation network \(\mathcal {G}_{f}\). A flow function f assigns a real value to each edge under two main constraints: the capacity limit and the flow conservation. The capacity limit constraint entails that the flow in each edge has to be less than or equal to the capacity (i.e. the weight) of this edge. If the flow equals the capacity, the edge is said saturated. The flow conservation constraint signifies that, at each node, the entering flow equals the exiting flow. Subject to these two constraints, the aim is to find the maximal flow from s to t in the flow network \(\mathcal {G}_{f}\). According to the graph construction rules provided by [25], the flow network for solving Problem (2) is composed of:

Figure 3 illustrates this graph construction on a small-size example. Computing a maximum flow from the source to the sink in this flow network saturates some edges, thus splitting the nodes a _i,j into two different groups: nodes that are reachable through a non saturated path from the source, and those that are not. Assuming that the source node s is labeled with 1, and the sink node t is labeled with 0, binary values are thus attributed to the edge labels x _i,j (secondarily, binary values are also assigned to the y labels of nodes v in the flow network), and this final labeling returns the set of selected edges \(\mathcal {E}^{*}\) which minimizes (2). We use the C++ code implementing a max-flow algorithm from [26].

Problem dimension reduction

This formula also provides a better insight into the role played by thresholding parameters λ _i,j . We now have a fast optimization strategy to generate a solution to the proposed variational formulation. One of the advantages of employing the BRANE Cut algorithm is the optimality guaranty of the resulting inferred network with respect to the proposed criterion. We next describe quantitative gains that can be achieved using BRANE Cut.

Validation datasets

Evaluation measures

where TP is the number of true positive, FP is the number of false positive and FN is the number of false negative. The Precision value indicates the proportion of correctly inferred edges compared to the total number of inferred edges. The Recall value reveals the proportion of correctly inferred edges compared to the total number of expected edges given by the gold standard. In order to evaluate and to rank the different tested methods, Precision-Recall (PR) curves are commonly used [8]. As the best results correspond to both high precision and high recall values, the Area Under the Precision-Recall Curve (AUPR) is an appropriate quantitative criterion to measure the quality of an inference method (higher is better).

Results on DREAM4

Computed AUPR in Table 1(a) highlight in bold that, globally, first and second best performances are always produced with BRANE Cut. Furthermore, each method tested (CLR, GENIE3, ND-CLR or ND-GENIE3) used as initialization exhibits an improved AUPR with BRANE Cut post-processing. Indeed, the average improvement reaches 10.6 % based on the CLR weights, 8.4 % for the GENIE3 weights, 5.9 % with ND-CLR weights and 7.2 % compared to the ND-GENIE3 weights, see Table 1(b).

We finally compare ND and BRANE Cut as post-processing methods on original weights. As shown in Table 1(c), BRANE Cut outperforms network deconvolution except for a practically unnoticeable degradation on the fifth network for GENIE3 weights. This relative improvement is essentially due to the fact that network deconvolution degrades results on the first two networks.

In the associated Precision-Recall curves, reported in the Additional file 1, we notice that the improvements of our results are mostly obtained in the first part of the curves, corresponding to a Precision greater than 50 % in the inference. Thus, such inferred graphs are expected to be more reliable for a biological interpretation. From this observation, looking at the AUPR for different Precision ranges, from the whole scale to precisions above 50 %, provides a finer assessment of the predictive power of inference methods. Thus, Fig. 4 highlights relative AUPR improvements for given Precision ranges. It illustrates that BRANE Cut improvement ratios over GENIE3 AUPR are clearly visible at higher Precision ranges, typically over 65 %.

Based on the AUPR criterion, we conclude that BRANE Cut outperforms state-of-the-art methods. Specifically, single-threshold results are sensibly refined by our approach, regardless of initial weights.

Results on the E. coli dataset

We now present the results of the BRANE Cut method on the real E. coli dataset. Our approach uses the normalized weights ω _i,j defined by CLR (with μ = 1000) and GENIE3 (with μ = 10). A discussion about the choice of the μ parameter is given in the Additional file 1. The parameter γ is set as in the previous section. The different Precision-Recall curves are reported in Fig. 5, to compare BRANE Cut to GENIE3 and CLR.

Best performance is expected toward the upper right side of Precision-Recall curves, with both high recall and precision. However, GRN Precision-Recall curves traditionally exhibit low Precision values over the whole curve on real datasets due to the difficulty in inferring accurate regulation relationships among large amounts of genes. For instance with the E. coli dataset, we observe that a recall below 0.05 corresponds to small inferred graphs, with less than 300 edges and a high precision (more than 75 %). Due to their higher predictive power and their readability, such small networks are often preferred by biologists. Hence, we focus on the upper-left part of the Precision-Recall curves in Fig. 5, emphasized in a close-up, corresponding only to high precision and small graphs. Here, BRANE Cut initialized with GENIE3 weights proves to be the best performer on smaller graphs (less than 100 edges corresponding to a recall below 0.02). However, graphs of larger size (up to a recall of 0.08) are more accurately reconstructed with CLR and BRANE Cut initialized with CLR weights. Again, the BRANE Cut approach improves the prediction results of both CLR and GENIE3.

Inferred network example on E. coli

An example of regulatory network on the E. coli dataset obtained with BRANE Cut, initialized with GENIE3 weights, is displayed in Fig. 6. The inferred network obtains a Precision score of 85 %, with a better predictive power than the network produced by the GENIE3 method alone. The binary network for GENIE3 is obtained by selecting edges having a weight higher than 0.707. This threshold renders a network with 85 % of Precision. In comparison to the reference, we discover 20 additional plausible regulatory interactions. Among these 20 predictions, ten were also predicted by the GENIE3 method, leading to ten predictions specific to BRANE Cut. By analyzing the predictions using STRING [29] and EcoCyc [30] databases, we observe that the predicted groups of genes were already identified as co-expressed and are known to belong to the same functional mechanism.

Results on DREAM5

We have evaluated BRANE Cut on three DREAM5 networks (1, 3 and 4) for which a validation exists. BRANE Cut parameters are initialized with the proposed heuristics and results are obtained using the validation procedure previously detailed. BRANE Cut outperforms CLR and GENIE3 by 7 % and 5 % respectively on Network 1. The improvement is 2.8 % and 2.1 % for Network 3. For the fourth network, the maximum Precision only reaches about 35 %. The AUPR computed with every method is exceptionally low. As such, the relative AUPR differences are insignificant, within the numerical precision. The detailed AUPR are given in the Additional file 1. Regarding the results in these additional datasets, the proposed heuristics lead to improvements over state-of-the-art.