{'fraction_non_alphanumeric': 0.03427202653318183, 'fraction_numerical': 0.013512268525627246, 'mean_word_length': 4.794306049822064, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
NExUS: Bayesian simultaneous network estimation across unequal sample sizes

Network-based analyses of high-throughput genomics data provide a holistic, systems-level understanding of various biological mechanisms for a common population. However, when estimating multiple networks across heterogeneous sub-populations, varying sample sizes pose a challenge in the estimation and inference, as network differences may be driven by differences in power. We are particularly interested in addressing this challenge in the context of proteomic networks for related cancers, as the number of subjects available for rare cancer (sub-)types is often limited. We develop NExUS (Network Estimation across Unequal Sample sizes), a Bayesian method that enables joint learning of multiple networks while avoiding artefactual relationship between sample size and network sparsity. We demonstrate through simulations that NExUS outperforms existing network estimation methods in this context, and apply it to learn network similarity and shared pathway activity for groups of cancers with related origins represented in The Cancer Genome Atlas (TCGA) proteomic data.


Introduction
The last decade has seen a proliferation of large, complex datasets that quantify molecular variables such as gene, protein, microbiome, and population-wide genetic variation. The Cancer Genome Atlas (TCGA) is a prime example of recent large-scale consortium-level efforts, which has generated multi-platform 'omics measurements from >10K patients across 32 common and rare cancer types. This allows for systematic investigations into the molecular mechanisms behind various oncogenic processes. Multiple studies have established that cancer initiation and progression are not outcomes of a single mutation within a gene or protein, but rather the result of a perturbation to co-ordinated networks and pathways that correspond to basic oncogenic processes such as cell proliferation and apoptosis (Hanahan and Weinberg, 2000;. Therefore, it is important to understand and characterize the underlying network dependency structures within and across cancers (Cho et al., 2016;Hristov and Singh, 2017). From a discovery standpoint, this is not only crucial to identify new cancer biomarkers and mechanisms, but also to distinguish the key molecular regulators of networks in different cancers (Sonabend et al., 2014;Carro et al., 2010).
Most existing studies focus on identifying the characteristics of molecular networks in individual cancer populations one at a time to understand tumor-specific molecular interactions (Gill et al., 2014;Creixell et al., 2015). On the other hand, researchers are beginning to recognize the critical importance of simultaneous analysis of similar tumor types (e.g., in terms of cell origin, organ location, biological evolution) to understand fundamental commonalities and differences (Weinstein et al., 2013;Tamborero et al., 2013). This has led to multiple pan-cancer studies that encompass the examination of various genomic and immunologic features in related cancers such as different squamous carcinomas (Campbell et al., 2018), gynecologic and breast cancers (Berger et al., 2018), gastrointestinal adenocarcinomas (Liu et al., 2018) and urologic cancers (Chen et al., 2017). In this article, we propose a network modeling approach for simultaneous analysis of genomic data from related tumor types.
Probabilistic graphical models (Lauritzen, 1996) are well-established statistical tools to conduct estimation and inference of network structures. In particular, Gaussian graphical models (Whittaker, 1990) have gained immense popularity because of their ability to capture global dependency structures. In high-dimensional settings, such as genomics, sparse graphical models are widely used to identify important nodes and interactions in a network consisting of a large number of genes/proteins. Various approaches for Gaussian graphical model estimation (Yuan and Lin, 2007;Friedman et al., 2008;Wang, 2012Baladandayuthapani et al., 2014) have been proposed over the years. These methods recently have been generalized to link network estimation across multiple populations using penalization-based (Danaher et al., 2014) or Bayesian approaches (Peterson et al., 2015). The former approach is not optimal for groups with differing levels of similarity, while the latter has scaling limitations to the prior specification. Scalability in higher dimensions also is an issue for other existing methods for joint estimation of sparse graphs (Guo et al., 2011;Chun et al., 2015). Along with learning the sparse network for each cancer type while borrowing strength across groups, another objective of our proposed method is to avoid artefactual differences in network sparsity due to differing sample sizes. In an application to pan-cancer network analysis, Kling et al. (2015) proposed an algorithmic method for joint network estimation with a group-specific penalty correction based on sample size. Their method, however, lacks a formal statistical justification of the sample-size penalty term and its effect on the sparsity of the estimated networks.
In this article, we propose a Bayesian approach for simultaneous Network Estimation across Unequal Sample sizes (NExUS). The main advantages of the NExUS method over existing methods include (i) explicit incorporation of sample size correction in network estimation which not only allows control of the sparsity but also improves borrowing of strength between cancer-specific networks to enable a balancing of statistical power across groups; (ii) the ability to quantify a network similarity index (NSI) which can be used as a global measure of network similarity among different cancer-specific networks; (iii) automatic selection of penalty parameters within a fully Bayesian framework, avoiding the need for the cross-validation step required by most frequentist graphical modeling approaches; and (iv) an efficient strategy for updating the precision matrices, making the method more scalable both in terms of the number of variables and the number of groups than existing Bayesian methods (Peterson et al., 2015;Kundu et al., 2018).
Our methods are motivated by and applied to The Cancer Proteomic Atlas (TCPA, Li et al., 2013) that has collected high-quality Reverse Phase Protein Array (RPPA) data across 32 cancer types. Many of these cancers can be divided into groups based on their histological origin, location, or similarity of biological and oncogenic processes. In this paper, we consider 4 different group of cancers, namely pan-gynecological, pan-kidney, pan-squamous, and pan-gastrointestinal. The simultaneous analysis of the proteomic networks of related cancers can help in robustly identifying the common network features. In addition, it can improve power for the discovery of features for rare cancers, which are shared with other cancers. We analyze the shared network structure of these cancers and identify shared pathway activity between all pairs of protein networks within the same group of cancers. We are able to establish both global trends in network similarity (e.g., we find that uterine carcinoma, a rare cancer, is more similar to cervical squamous cell carcinoma and endocervical adenocarcinoma, than to other gynecological cancers), as well as pathway-specific activity sharing (e.g., we find that the hormone receptor pathway has common activity across the pan-gynecological group, while the RAS/MAPK and RTK pathways have shared activation across the various types of kidney cancers).

The NExUS method
Data structure and notation Let Xc represent the nc × p matrix of observed protein or gene expression for the c th cancer type of interest, where c = 1, 2, . . . , C. We assume that the same set of p proteins are observed for all subjects, but allow the sample sizes nc for each cancer type to differ. We assume that the data for each subject i follows a multivariate normal distribution with mean vector 0 ∈ R p and cancer-specific precision matrix Θc. The multivariate normal distribution is parametrized using the precision matrix Θc (rather than the covariance matrix Σc = Θ −1 c ) since there is a direct correspondence between the precision matrix and the conditional independence graph among variables. Specifically, in a Gaussian graphical model, entry θ c ij in the precision matrix Θc is exactly 0 if and only if the corresponding proteins i and j are conditionally independent for that cancer type i.e. they are not connected by an edge within the conditional dependence network (Dempster, 1972). Our statistical objectives are to learn a sparse network for each cancer type with an approach that both allows for borrowing of information across cancers, and avoids artefactual differences in network sparsity due to differing sample sizes. We now describe the prior formulation that allows us to achieve these goals.

×1
We put independent exponential priors with mean γ −1 on the diagonal elements (i.e., on each element of θ D ), and a gamma hyperprior on the parameter γ, with shape parameter αγ and rate parameter βγ.
To enable borrowing of strength while estimating sparse versions of the precision matrices for the set of cancers of interest, we introduce the shrinkage parameters λ 1 and λ 2 , where λ 1 controls the shrinkage within each cancer type and λ 2 induces similarity across the cancer types. Instead of taking the same λ 1 for each cancer type and the same λ 2 for each pair of cancers, we consider a more flexible scenario where the withincancer and cross-cancer shrinkage parameters are dependent on cancer type. Let λ c 1 denote the individual shrinkage parameter for cancer type c, 1 ≤ c ≤ C, and λ cc ′ 2 denote the cross-penalty across cancer types c and c ′ , 1 ≤ c < c ′ ≤ C. Following Kyung et al. (2010), the conditional prior of θ ND is given by (1). Here the first term aims to achieve sparsity within each cancer and the second term to achieve similarity across cancers.
Network similarity index: Taking a closer look into the role of the crossgroup penalty parameter in the model, it is evident that higher values of (λ cc ′ 2 ) 2 encourage more similarity between the networks of cancer c and c ′ i.e., the higher the value of (λ cc ′ 2 ) 2 the closer the network structure (edges) will be between two cancers (we will illustrate this in our simulations and real data analyses). Therefore within a group of related cancers, the similarity between cancer types c and c ′ can be estimated by the penalty (λ cc ′ 2 ) 2 , which we term the NSI. Within each group of related cancers, we transform the NSI values of all pairs of cancers to be within the unit interval using a linear monotonic mapping such that the pair of cancers with the lowest and the highest NSI values are mapped to 0 and 1 respectively. We term these transformed NSI values as normalized network similarity index (NNSI).

Incorporating sample size adjustment in the priors
When estimating networks for multiple cancer types using existing methods, network sparsity will vary depending upon sample size, with larger sample sizes generally resulting in denser graphs (more edges). From a pan-cancer context, however, sparsity should not be dependent on the sample size -since large cancers will overwhelm under-sampled or rare cancers. To counteract this dependency, we design priors using shrinkage parameters to mitigate the sample size effect. Since larger values of λ c 1 encourage more shrinkage of the elements of the precision matrix, and hence more sparsity, our objective is to have a larger prior mean of (λ c 1 ) 2 for the cancer types with larger sample sizes, and a smaller prior mean for rare cancer types with smaller sample sizes. Similarly, the shrinkage parameters λ cc ′ 2 encourage similarity between cancer types c and c ′ . To achieve these goals, we put gamma priors on the squared terms (λ c 1 ) 2 with shape parameter α 1 and rate parameter β c 1 , for 1 ≤ c ≤ C, and gamma priors on (λ cc ′ 2 ) 2 with shape parameter α 2 and rate parameter β cc ′ 2 , for 1 ≤ c < c ′ ≤ C. Following Kling et al. (2015), we define β c 1 and β cc ′ 2 as follows: where β 1 , β 2 > 0. Letn represent the average sample size across the C cancer types. We define the effective sample size of the c th cancer type as n e c =n δ n (1−δ) c , where 0 < δ < 1. Hence, we get the prior mean of (λ c 1 ) 2 to be α 1 β 1 (n e c ) 2 , and the prior mean of (λ cc ′ 2 ) 2 to be . At δ = 0, n e c = nc, so the prior mean depends only on the sample size for cancer type c, while at δ = 1, n e c =n, so the prior mean is equal across all cancer types. As δ approaches 0, relatively more sample size correction is induced by the prior.
An illustrative example: Consider three cancer types A, B, and C with sample sizes n A = 50, n B = 100, and n C = 200. Then for δ < 1, our objective is to have prior mean of (λ A 1 ) 2 to be the smallest and the prior mean of (λ C 1 ) 2 to be the largest, while prior mean of (λ B 1 ) 2 should be in between those values. Also, for δ < 1, for similarity shrinkage parameters, our objective is to have the lowest prior mean among the two cancers with the smallest sample sizes (here A and B) and the highest prior mean among the two cancers with the highest sample sizes (here B and C). The variation of the prior means of the shrinkage parameters for this scenario is plotted as a function of δ in Figure 1. Note that under the proposed way of prior construction, the desired order between the prior means of the shrinkage parameters are maintained for δ < 1. Also note that as δ increases to 1, the penalty terms converge in both scenarios.

MCMC and edge selection
After introducing latent variables T 2 and Ω 2 , with an appropriate choice of prior, it can be shown that the conditional prior of equation (1) can be re-expressed as the normal distribution where Σ Θ N D can be formulated in terms of T 2 and Ω 2 (see supplementary material Section S1). This property is important as it allows us to construct a computationally efficient Gibbs sampler.
Since the posterior distribution is not tractable, we rely on Markov Chain Monte Carlo (MCMC) sampling to obtain a sample from the posterior distribution. Specifically, we are able to construct a Gibbs sampler using the posterior conditional distributions for each parameter or set of parameters. Here is a high-level outline of the updates that take place in each iteration of the algorithm. For additional details on the joint posterior, posterior full conditional distributions, and sampling steps, see Sections S2-S4 in the supplementary material.
• Update the ith column of Θc following the approach in Wang (2012) by sampling a scalar from a gamma distribution and a vector from a multivariate normal, and then applying an appropriate transformation for c = 1, . . . , C and i = 1,…, p.
• Update the latent variable T 2 by sampling the element (τ c ij ) −2 from an inverse-gamma distribution for c = 1, . . . , C and 1 ≤ i < j ≤ p.
• Update the latent variable Ω 2 by sampling the element (ω cc ′ ij ) −2 from an inverse-gamma distribution for 1 ≤ c < c ′ ≤ C and 1 ≤ i < j ≤ p.
• Update γ by sampling from a gamma distribution.
Since the shrinkage prior results in sampled values for the precision matrices that have entries, which are very close to 0, but not exactly 0, thresholding is required to represent the sampled precision matrices as sparse. Therefore, to identify selected edges in the estimated networks, we choose a cut-off κ = 0.05 on the elements of the partial correlation matrices (derived from corresponding precision matrices). Instead of simply thresholding the elements of the posterior mean partial correlation matrices, we consider a slightly different probabilistic approach. From the posterior sample, we specifically calculate the posterior probability of each non-diagonal entry of the partial correlation matrices being greater than the cut-off κ. The edges corresponding to the entries with posterior probability of inclusion greater than 0.5 across all MCMC iterations are then selected (see supplementary material Section S5 for additional details on the edge selection procedure).

Simulation studies
In this section, we compare the performance of the proposed method to several existing methods on a simulated data set. To construct the simulated data, we consider a version of the set-up on p = 20 variables and C = 4 groups from section 5.1 of Peterson et al. (2015), modified so that all of the graphs have at least some shared structure and there are unequal sample sizes across groups. Specifically, the first precision matrix Θ 1 is ( 1 ) 2 ( 2 ) 2 ( 2 AB ) 2 (n A =50, n B =100) constructed so that all its diagonal entries are 1, For each cancer type c = 1, . . . , 4, we then generate nc observations by sampling from the multivariate normal distribution with mean zero and precision matrix Ωc.
For the methods comparison, we consider a variety of Bayesian and frequentist methods. In the Bayesian framework, in addition to the proposed NExUS method, we consider the Bayesian graphical lasso (BGL) and Bayesian adaptive graphical lasso (BayesAdapGL, Wang, 2012). The values of the tuning parameters for the BGL and BayesAdapGL are the same as those considered in Wang (2012). For edge selection, however, we use the same thresholding approach as NExUS. For NExUS, we use the following hyperparameter setting: α 1 = 1, α 2 = .1, β 1 = .1n 2 , β 2 = 1n 2 , αγ = 1, and βγ = 1. For all three Bayesian approaches, we discard the first 5,000 iterations as burn-in, and use the following 15,000 iterations as the basis for inference.
In the frequentist framework, we consider both single graph and joint estimation methods. Specifically, the single graph methods we compare to are the graphical lasso (FreqGL) (Friedman et al., 2008), adaptive graphical lasso (FreqAdapGL), and the SCAD penalized lasso (Fan et al., 2009). The joint estimation methods we consider are the group graphical lasso (GGL) and fused graphical lasso (FGL) (Danaher et al., 2014). Please see Section S6 of the supplementary material for further details regarding the application of these methods.  To compare the performance of these methods, we primarily rely on the area under the ROC curve (AUC), as it provides a single summary of network learning accuracy. In Table 1, we summarize the AUC values for the various methods for each of the four simulated graph structures. The reported value is the mean of AUC values estimated from 100 simulated datasets. The standard error values are reported inside parentheses. In these results, NExUS is consistently one of the top two performing methods. FGL also performs well for groups C 1 and C 2 , which have relatively smaller sample sizes and benefit most from the enforcement of similarity across graphs, but its performance is relatively worse for C 3 and C 4 . Some of the single graph methods (BGL and BayesAdapGL) become competitive when applied to the group with the largest sample size, but lack power in the classes with smaller sample sizes. In Table 2, we show the accuracy of learning shared edges, as summarized by the AUC for each pair of cancer types based on 100 simulated datasets. Again, NExUS demonstrates the best overall performance, with FGL and GGL also competitive. FGL and GGL, while achieving good performance overall, are somewhat less flexible than NExUS, as they rely a single value for λ 1 for each class and a single value of λ 2 for each pair of classes. The procedure of finding the optimal values of those penalty parameters is computationally challenging and not scalable for networks with a large number of nodes.
To gain insight into how NExUS improves on separate Bayesian inference, we plot the true positive rate (TPR), false positive rate (FPR), and Mathew's Correlation Coefficient (MCC) across a range of values for the threshold κ on the partial correlation matrix elements in the supplementary material Figure S1. This figure demonstrates that the NExUS method Table 2. Mean AUC values of shared edge learning across 100 simulated data sets, with standard errors in parentheses. The results for the top two performing methods for each cancer type are highlighted in bold.

NExUS analyses of pan-cancer proteomic data
We demonstrate the utility of NExUS to conduct pan-cancer analyses of proteomics data collected across various cancer types. We focus on the four sub-groups of related cancers: pan-gynecological (pan-gynae), pankidney, pan-squamous, and pan-gastrointestinal (pan-GI). Gynecologic cancers have similar embryonic origins, since female hormones influence their development (Berger et al., 2018). The pan-gynae group consists of breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), ovarian serous cystadenocarcinoma (OV), uterine corpus endometrial carcinoma (UCEC), and uterine carcinosarcoma (UCS). Since kidney cancers originate from the cells of the outer layer of the kidney (the renal cortex), we consider the following cancers in the pan-kidney group: kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), and kidney renal papillary cell carcinoma (KIRP). Squamous cell carcinomas (SCCs) arise from the epithelia of the aerodigestive and genitourinary tracts and share some histological characteristics that can be used for predicting the site of origin, clinical behavior, cause, prognosis, or optimal therapies (Campbell et al., 2018).
In the pan-squamous group, we consider esophageal squamous-carcinoma (ESCA(sq)), head and neck squamous cell carcinoma (HNSC), and lung squamous cell carcinoma (LUSC). Adenocarcinomas of the GI tract share similar endodermal developmental origins along with exposure to common insults that promote the tumor formation (Liu et al., 2018); hence the pan-GI group consists of colon/rectum adenocarcinoma (CORE), esophageal adeno-carcinoma (ESCA(ad)), and stomach adenocarcinoma (STAD).
Our proteomics data arises from The Cancer Proteome Atlas (TCPA, Li et al., 2013), which provides protein abundance data for TCGA samples, and consists of RPPA-based quantifications using antibodies that cover functions including proliferation, DNA damage, polarity, vesicle function, EMT, invasiveness, hormone signaling, apoptosis, metabolism, immunological and stromal function, as well as other critical cellular signaling pathways (Akbani et al., 2014). The proteins profiled allow for a focused exploration of the various functional mechanisms underlying oncogenic processes across tumor types.
A challenge in analyzing this data is that sample sizes for different cancers vary considerably, from 48 for UCS, to 879 for BRCA (see supplementary materials Table S4 for the number of samples available for each cancer type). Instead of estimating the proteomic networks for each cancer separately, joint estimation and borrowing of strength across the networks gives us a broader picture of the similarities and dissimilarities among the cancers belonging to the same group. Specifically, the cross-group penalty terms in NExUS help in identifying some of the weak signals (edges) for the cancer networks with smaller sample sizes when those edges are shared by other cancers of the same group, and the sample size adjustment encourages networks with more similar levels of sparsity across cancers belonging to the same pan-cancer group. In addition to learning the network structures per cancer type, we also are able to estimate the global proteomic network-based similarities and the pathway-specific similarities among the cancers belonging to the same group.
For the analysis, we use the same values of the prior parameters as in the simulation studies. We perform 20,000 iterations, discarding the first 5,000 iterations as burn-in for each pan-cancer group. We adopt the same posterior edge selection strategy as in the simulation studies. The posterior selected networks for each cancer type are shown in Figures S8-S21 of the supplementary materials. Here we focus on high-level take-aways of these pan-cancer analyses, in particular, on the global and pathway-specific similarity measures.

Global proteomic similarity between cancers
To understand the extent of shared network structure among the cancers belonging to the same group, we estimate the normalized network similarity indexes (NNSI) for each pair of cancers. We plot the NNSI in Figure 2 for all pan-cancer groups. In the pan-gynae group we observe 3 distinct clusters among the pairs of cancers based on NNSI: CESC-UCS are the closest; OV-UCS, CESC-OV, UCEC-UCS, CESC-UCEC, and OV-UCEC belong to the second cluster with intermediate NNSI values; and the third cluster, which represents the most distant pairs of cancers, contains BRCA-OV, BRCA-UCS, BRCA-CESC and BRCA-UCEC. These results show that the protein network for the rare cancer UCS is more similar to that of CESC and OV cancers as compared to other gynae cancers. BRCA has the least network similarity with other gynae cancers, suggesting BRCA has distinctive network features. In the pan-kidney group, the KICH and KIRP networks are the closest, while KIRC and KIRP networks are least similar. For the pan-squamous and pan-GI groups, ESCA(sq)-LUSC and ESCA(ad)-STAD are most similar while HNSC-LUSC and CORE-STAD are the least similar. To provide further insight into the meaning of the network similarity index, we examined the relationship between it and the L 1 distance between pairs of precision matrices, and found the relation to be approximately linear, where pairs of cancer with relatively high values of the network similarity index have relatively lower L 1 distances between them (see supplementary material Figure S5).
Finally, to visualize the extent of shared structure across each group of cancers, we construct heatmaps of the posterior edge inclusion probabilities ( Figure 3). The ordering of the cancer types, which was determined by hierarchical clustering of the edge probabilities, supports the NNSI results, with pairs of cancer scoring higher on the NNSI grouped more closely within the estimated hierarchy.

Pathway-specific similarity between cancers
To understand important aspects of the network similarity between cancer types, we conducted a deeper investigation of the proteomic networks using available pathway information. We focus our exploration on 12 wellestablished curated pathways with translational relevance (Akbani et al., 2014;Ha et al., 2018): apoptosis, cell cycle, DNA damage response, EMT, hormone receptor, hormone signaling breast, PI3K/AKT, RAS/MAPK, RTK, TSC/mTOR, breast reactive, and core reactive. A few of the proteins profiled are shared by more than one pathway (see supplementary material Table S5). One of our main motivations for deconvolving proteomic activity into pathways is to investigate which pathways are activated in each of the cancer types, and thereby gain a better understanding of the mechanistic and regulatory sharing of information between proteomic pathways.
To quantify pathway-specific activation, we estimate the proportion of shared edges within each pathway and across each pair of pathways for all cancers belonging to the same group. The heatmap for the proportion of shared edges within and across pathways for each pair of cancers is given in Figures S6-S7 of the supplementary materials. In general, the proportion of shared edges are higher within pathways than across different pathways, which is along expected lines. This makes biological sense, as core pathway activities are likely to be preserved across different cancer lineages. To identify the strongest shared pathway activity for each pair of cancers within each group, we computed the percentage of shared edges within each pathway for each pair of cancer types. In Figure 4, we include edges based on the ranking of these percentages of edges shared: specifically, we include a colored link between cancer types in the figure for the 40% highest sharing percentages across all pathways and cancer pairs for all cancer types, except pan-gynae. As the pan-gynae group has a larger number of cancer types included, we focus on the top 20% there to improve legibility. For the pan-gynae group, the hormone receptor pathway and the core reactive pathway are the most active shared pathways across the pairs of cancers. The hormone receptor pathway is important for OV cancer, where hormonal-based systemic therapies are used to treat ovarian stromal tumors (Dacheux et al., 2013). As UCS is a rare cancer, it is challenging to pinpoint its active pathways and to draw network-based inference. We are able to identify the cell cycle pathway as one of the top three actively shared pathways between UCS and UCEC. This is clinically relevant as a cell cycle pathway inhibitor has been identified as one of the therapeutic options for UCS (Cherniack et al., 2017). Finally, PI3K/AKT pathway activity is shared between the BRCA and UCEC as well as BRCA and OV cancer types. TCGA Research Network (2012b) mentions many components of PI3K pathway were amplified in basal-like breast cancers. Figure 4: Shared pathway activity: Proportion of shared edges within pathways across all pairs of cancers, with edges denoting top 20% (for gynae) and 40% (for other cancer groups) of all pathways for each pair of cancers within each group are plotted for (a) pan-gynae, (b) pan-kidney, (c) pan-squamous, (d) pan-GI cancer groups.
For the pan-kidney group, our analysis identifies the TSC/mTOR, RAS/MAPK, and RTK pathways as the top 3 actively-shared pathways. We also observe shared pathway activity for the PI3/AKT, EMT, core reactive, cell cycle, and hormone receptor pathways. Chen et al. (2016) and TCGA Research Network (2013) recognized the PI3K/AKT and mTOR pathways to be active for renal cell carcinoma (RCC) and clear cell renal cell carcinoma (ccRCC). In addition, Gibbons and Creighton (2018) found EMT as one of the important activated pathways in kidney cancers.
For the pan-squamous group, the cell cycle, RTK, and PI3K/AKT pathways are activated between all pairs of cancers. Other activated pathways are TSC/mTOR, EMT, core reactive, and hormone receptor. TCGA Research Network (2015) and TCGA Research Network (2012a) described the PI3K, RAS, and AKT pathways as major pathways influencing HNSC and LUSC, respectively. TCGA Research Network (2015) and TCGA Research Network (2012a) also reported that the cell cycle and RTK pathways play an important role in HNSC and LUSC. Finally, for the pan-GI group, we find the hormone receptor, RTK, EMT, and core reactive pathways as the top 3 activated shared pathways.

Conclusion and Remarks
In this article, we propose NExUS, a fully Bayesian method for estimating a group of related networks using Gaussian graphical modeling. The incorporation of a penalty on dissimilarity across the precision matrices in the prior specification enables borrowing of strength across the networks. This aspect of the model is especially helpful for identifying edges in the networks with smaller sample sizes. In addition, the resulting NSI helps order the relative closeness of each pair of networks within the set analyzed. Another novel feature of NExUS is the inclusion of a sample size correction which encourages similar sparsity levels across networks with different sample sizes. Our simulation studies show that NExUS outperforms other existing methods for individual and joint estimation of networks. NExUS is motivated by the TCGA-based RPPA dataset, wherein we estimate the pan-cancer proteomic networks for 4 groups of related cancers lineages. The existence of rare cancers (e.g. UCS in pan-gynae, and ESCA(ad) in pan-GI) makes the proposed model particularly appropriate for this application.
We developed the NExUS model under Gaussian assumptions and applied it for datasets containing a moderate number of proteins obtained using targeted profiling. In the future, it can be modified for skewed data and discrete variables, potentially following the methods proposed in Bhadra et al. (2018). The method would then be applicable to non-normal data sets such as mutation or copy number variation. Finally, more efficient algorithms could be developed to allow NExUS to scale to a larger number of variables, enabling an application to data sets such as the Clinical Proteomic Tumor Analysis Consortium (CPTAC), which was acquired using untargeted proteomics, and therefore has a much higher dimensionality.