In a gene network model, genes are represented by “nodes” and associations between them are represented by “edges.” Given measurable molecular units each corresponding to one gene (e.g. the encoded protein or mRNA), a network, or graph, can be constructed from their observed values.

Consider n observed values of the multivariate random vector x=(X1,…,Xp)T of node attributes, with each entry corresponding to one of p nodes. If we assume multivariate Gaussian node attributes, with an n×p observation matrix X with i.i.d. rows x1,…,xn∼N(0,Σ), a “partial correlation network” can be determined by estimating the inverse covariance matrix, or precision matrix, Θ=Σ−1. Specifically, the partial correlation between nodes i and j, conditional upon the rest of the graph, is given by ρij|V\{i,j}=−θijθiiθjj, where the θij’s are the entries of Θ and V the set of all node pairs (Lauritzen 1996). The partial correlations coincide with the conditional correlations in the Gaussian setting. Because correlation (resp. partial correlation) equal to zero is equivalent to independence (resp. conditional independence) for Gaussian variables, a conditional independence graph can thus be constructed by determining the non-zero entries of the precision matrix. To ensure invertibility, the precision matrix also required to be positive definite, Θ≻0.

In high-dimensional settings, the sample covariance matrix S=1n−1XTX is rarely of full rank and thus its inverse cannot be estimated directly. It is common to assume sparse network, meaning the number of edges in the edge set E is small relative to the number of potential edges in the graph (i.e. the sparsity measure 2|E|/(p2−p) is small). Penalized methods, such as the Glasso (Friedman et al. 2008) are well established for sparse Gaussian graphical model estimation. In the case of there being multiple (related) datasets available, such as from different tissue types, rather than estimating each network separately much statistical power could be gained by sharing information across networks through a joint approach.

Assume a network inference problem with K groups. We let {Θ}=(Θ(1),…,Θ(K)) be the set of their (unknown) precision matrices, and assume that the set of ∑k=1Knk observations are independent. We aim to solve the penalized log-likelihood problem (Danaher et al. 2014) (1){Θ^}=argmax{Θ≻0}{∑k=1Knk[log(det(Θ(k)))−tr(S(k)Θ(k))]−P({Θ})}, where S(k) is the sample covariance matrix of group k and P(⋅) is a penalty function. In (1), det(⋅) denotes the determinant and tr(⋅) denotes the trace. The JGL employs the fused penalty function (2)P({Θ})=λ1∑k=1K∑i≠jabs(θij(k))+λ2∑k<k′||Θ(k)−Θ(k′)||1, where λ1 and λ2 are positive penalty parameters, abs(⋅) denotes the absolute value function and ||⋅||1 denotes the L1 penalty. This penalty applies L1 penalties to each off-diagonal element of the K precision matrices as well as to the differences between corresponding elements of each pair of precision matrices. The parameter λ1 controls the sparsity, and the similarity parameter λ2 controls the degree to which the K precision matrices are forced toward each other, encouraging not only similar network structures but also similar precision matrix entries. This way, the dependency between observed datasets is modeled through their underlying graphical structures. It is important to note that a necessary assumption of the resulting model is that given their respective graphical structures, observations are independent across and within each dataset. The current penalty parameter selection approach for λ1 and λ2 is based on the AIC (Danaher et al. 2014), and while suitable for determining network similarities, likelihood-based selection criteria can lead to severe under- or over-selection and thus poor performance in high-dimensional settings (Liu et al. 2010).

To improve the performance of the JGL with the fused penalty for omics applications and other high-dimensional problems, we propose the stabJGL algorithm for stable sparsity and similarity selection in multiple-network reconstruction. StabJGL jointly estimates multiple networks by leveraging their common information and gives a basis for deeper exploration of their differences, as shown in Fig. 1. Below we outline the algorithm, which comprised two steps: (i) selecting the sparsity parameter λ1 in the fused penalty (2) based on the notion of model stability, and (ii) selecting the similarity parameter λ2 based on model likelihood. The full StabJGL algorithm is given in Supplementary Algorithm S1.

We first assess the performance of stabJGL on simulated data. We compare the network reconstruction ability of stabJGL to that of state-of-the-art methods, including the JGL with the fused penalty (FGL) and group penalty (GGL) with penalty parameters selected with the default AIC-based criterion (Danaher et al. 2014). To assess the performance of another selection criterion specifically designed for high-dimensional graph selection, we also consider FGL with penalty parameters tuned by the extended BIC for multiple graphs (7). We further include the Bayesian spike-and-slab JGL (SSJGL) of Li et al. (2019b), as well as the Glasso of Friedman et al. (2008) tuned by StARS (Liu et al. 2010). The latter estimates each network separately. We generate data that closely resembles our omics application of interest, featuring partial correlations between 0.1 and 0.2 in absolute value, while also exhibiting the “scale-free” property—a typical assumption for omics data where the “degree distribution” (i.e. the distribution of the number of edges that are connected to the nodes) adheres to a power-law distribution (Chen and Sharp 2004). We consider a wide range of settings, simulating K∈{2,3,4} networks with p∈{100,200,300,1000} nodes. We manipulate the degree of similarity in their “true” graphical structures to assess the performance of the method over a wide range of scenarios. We then apply different network reconstruction techniques to determine the networks from the data and report the results averaged over N=100 replicates. The parameter specifications for the different methods are given in Supplementary Section S2. We also investigate the effect of the variability threshold β1 in stabJGL on the results in a setting with p=100 nodes and K=2 networks. Finally, to compare the scalability of the respective methods, we consider the time needed to infer networks for various p and K. Further details and code for the simulation study can be found at https://github.com/Camiling/stabJGL_simulations.

Estimation accuracy is assessed with the “precision” (positive predictive value), and the “recall” (sensitivity). The precision gives the fraction of predicted edges that were correct, while the recall is the fraction of edges in the true graph that were identified by the inference. Because the sparsity of estimated networks will vary between methods, the precision–recall trade-off should be taken into consideration. In general, the recall will increase with the number of selected edges while the precision will decrease. Since sparsity selection is a main feature of our proposed method, we do not consider threshold-free comparison metrics, such as the AUC. We therefore put emphasis on the following characteristics in our comparative simulation study; (i) suitable sparsity level selection, (ii) utilization of common information at any level of network similarity, i.e. inference improves with increased network similarity, and (iii) a suitable precision–recall trade-off that overly favors either measure.

The results for K=2 and K=4 networks with p=100 nodes are summarized in Fig. 2. The full tables of results for these settings and the additional settings with K∈{2,3,4} networks and p∈{100,200,300,1000} nodes are shown in Supplementary Section S2, where details on selected sparsity and penalty parameters are given for all methods. Additional simulations investigating the effect of random noise on the graph reconstruction performance of the different methods are also given in Supplementary Section S6. In Fig. 2, FGL tuned with eBIC did not select any edges in any of the settings and is therefore not shown. In all settings considered, the FGL and GGL with the default AIC-based penalty parameter selection strongly over-select edges. This leads to high recall, but very low precision. Second, they do not appear to sufficiently utilize network similarities; the performance of the two methods, particularly GGL, differs little between completely unrelated and identical networks. Notably, in all cases, the selected value of λ2 is smaller for FGL and GGL tuned by AIC than it is for stabJGL. Consequently, similarity is not sufficiently encouraged even in settings where the networks are identical. The AIC criterion does not seem to provide sufficient penalization to encourage suitable sparsity and similarity. On the other hand, the alternative eBIC criterion gives extremely sparse FGL estimates, selecting an empty graph in most settings, i.e. no edges. Although the extended BIC is developed specifically for graphical model selection, likelihood-based criteria for sparsity selection tend to perform poorly in high-dimensional settings and risk both severe under- and over-selection (Foygel and Drton 2010). This issue is avoided in the stabJGL algorithm as the eBIC only is used to select similarity and not sparsity.

The Bayesian SSJGL tends to select very few edges, leading to high precision but low recall. Its performance deteriorates drastically as the network differences increase, leading to extremely low recall. This implies a lack of flexibility to adapt to varying network similarity levels, as has previously been observed (Lingjærde et al. 2022). Out of all the joint methods, stabJGL gives the most accurate sparsity estimate. This ensures that we neither get very low precision like FGL and GGL tuned by AIC, nor very low recall like SSJGL and FGL tuned by eBIC. Replacing the mean of the λ1’s over the N=100 replicates by alternative measures, such as the median yields the same reported penalty parameter values and hence results in all settings considered in Fig. 2. The mean and median could be expected to differ slightly if a finer grid of λ1 values were to be considered, but our findings suggest they would yield highly similar results. StabJGL also appears to adapt well to the similarity between networks, with the prediction accuracy increasing with the number of shared edges. As a result, the method either outperforms the Glasso tuned by StARS for highly similar networks or performs comparably to it for unrelated networks. The similar performance for unrelated networks can be explained by the fact that the sparsity controlling penalty parameters of both methods are tuned with a stability-based approach. The results suggest that stabJGL adapts well to the level of similarity and hence can be employed agnostically in settings where there is no prior knowledge about the degree to which the networks have shared structures.

A key question is whether stabJGL can achieve as high precision as the methods that give sparser networks (i.e. SSJGL) by using a lower variability threshold. Similarly, we want to see if stabJGL can achieve as high recall as the methods that infer more edges (i.e. FGL and GGL). To investigate this, we consider the same setting as in Fig. 2 with K=2 networks. Supplementary Figure S4 compares the performance of stabJGL for different values of the variability threshold β1 to the other methods, and we find that by decreasing (resp. increasing) β1 we can obtain at least as high precision (resp. recall) as the other methods at any level of similarity. This illustrates that the method can be adapted to reflect the priorities of the user (i.e. concern for false positives versus false negatives). For most applications, a middle-ground value, such as 0.1 yields a good balance between false positives and false negatives as demonstrated in the simulations.

Figure 3 shows the CPU time used to jointly infer networks for K∈{2,3,4} networks and various numbers of nodes p, with n∈{100,150} observations, for the JGL with the fused penalty (FGL) with penalty parameters tuned with the AIC and stabJGL with the same parameter specifications as in the previously described simulations. Due to an efficient parallelized implementation, stabJGL has an almost identical runtime to FGL when the same number of λ1 and λ2 values are considered. Thus, the increased estimation accuracy of stabJGL does not come at a computational cost. It is important to note that due to the generalized fused lasso problem having a closed-form expression for the updates in the special case of K=2 (Danaher et al. 2014), stabJGL is substantially faster for only two networks than for K>2. As stabJGL uses the fused penalty this comparison is the most relevant, but a runtime comparison of all methods considered in our simulation study can be found in Supplementary Fig. S1. In the Supplementary Material, we also demonstrate that stabJGL can be applied to problems with p≥2500 nodes or K=10 networks within reasonable time (Supplementary Fig. S3).

We perform a proteomic network analysis of Reverse Phase Protein Array (RPPA) data from The Cancer Genome Atlas (TCGA) across different pan-Cancer tumor types (Cancer Genome Atlas Network and others 2012). In a large proteomic pan-Cancer study of 11 TCGA tumor types, Akbani et al. (2014) identified a major tumor super cluster consisting of hormonally responsive “women’s cancers” [Luminal breast cancer (BRCA), ovarian cystadenocarcinoma (OVCA), and uterine corpus endometrial carcinoma (UCEC)]. Our objective is to map the proteomic network structure of the respective tumor types, so that we can get a better grasp of the common mechanisms at play in the hormonally responsive tumors. We are also interested in highlighting the differences.

We consider RPPA data from Luminal BRCA (n=273), high-grade serous OVCA (n=412), and UCEC (n=404). All data are downloaded from the UCSC Xena Browser (Goldman et al. 2020). The data are measured with p=131 high-quality antibodies that target (phospho)-proteins. To alleviate batch effects, the RPPA data are normalized with replicate-base normalization (Akbani et al. 2014). We use stabJGL to jointly estimate the proteomic networks of the respective tumor types and interpret the results and their implications. We compare the output with that obtained with the FGL of Danaher et al. (2014) with the default penalty parameter tuning with AIC as described in Section 3.1. Further details and code for the analysis are given at https://github.com/Camiling/stabJGL_analysis.