Robust Bayesian Graphical Regression Models for Assessing Tumor Heterogeneity
in Proteomic Networks–Robust Bayesian Graphical Regression Models for
Assessing Tumor Heterogeneity in Proteomic Networks December
10.1111/j.1541-0420.2005.00454.x

# Robust Bayesian Graphical Regression Models for Assessing Tumor
Heterogeneity in Proteomic Networks

Tsung-Hung Yao1,∗<EMAIL_ADDRESS>Yang Ni2 and Anindya Bhadra3 and Jian Kang1
and Veerabhadran Baladandayuthapani1,∗∗
1Department of Biostatistics<EMAIL_ADDRESS>University of Michigan at Ann
Arbor  Ann Arbor  Michigan  U.S.A.
2Department of Statistics  Texas A&M University  College Station  Texas
U.S.A.
3Department of Statistics  Purdue University  West Lafayette  Indiana  U.S.A

(Received October 2023\. Revised February 2024\. Accepted March 2024.; 2024)

###### Abstract

Graphical models are powerful tools to investigate complex dependency
structures in high-throughput datasets. However, most existing graphical
models make one of the two canonical assumptions: (i) a homogeneous graph with
a common network for all subjects; or (ii) an assumption of normality
especially in the context of Gaussian graphical models. Both assumptions are
restrictive and can fail to hold in certain applications such as proteomic
networks in cancer. To this end, we propose an approach termed robust Bayesian
graphical regression (rBGR) to estimate heterogeneous graphs for non-normally
distributed data. rBGR is a flexible framework that accommodates non-normality
through random marginal transformations and constructs covariate-dependent
graphs to accommodate heterogeneity through graphical regression techniques.
We formulate a new characterization of edge dependencies in such models called
conditional sign independence with covariates along with an efficient
posterior sampling algorithm. In simulation studies, we demonstrate that rBGR
outperforms existing graphical regression models for data generated under
various levels of non-normality in both edge and covariate selection. We use
rBGR to assess proteomic networks across two cancers: lung and ovarian, to
systematically investigate the effects of immunogenic heterogeneity within
tumors. Our analyses reveal several important protein-protein interactions
that are differentially impacted by the immune cell abundance; some
corroborate existing biological knowledge whereas others are novel findings.

###### keywords:

Bayesian graphical models; Cancer; Conditional sign independence; Covariate-
dependent graphs; Protein-protein interactions.

††volume: 64

##  1 Introduction

Graphical models are ubiquitous and powerful tools to investigate complex
dependency structures in high-throughput biomedical datasets such as genomics
and proteomics (Airoldi, 2007). They allow for holistic exploration of
biologically relevant patterns that can be used for deciphering biological
processes and formulating new testable hypotheses. However, most existing
graphical models make one of the two canonical assumptions: (i) a homogeneous
graph that is common to all subjects; or (ii) an assumption of normality as in
the context of Gaussian graphical models (Ni et al., 2022). However, in some
biomedical applications such as inference of proteomic networks in cancer,
both assumptions are violated, as we illustrate next.

Proteomic networks and tumor heterogeneity. Proteins control many fundamental
cellular processes through a complex but organized system of interactions,
termed protein-protein interactions (PPIs) (Cheng et al., 2020). Moreover,
aberrant PPIs are associated with various diseases including cancer, and
investigating PPI can lead to effective strategies and treatments, including
immunotherapies, tailored to different individuals (Cheng et al., 2020).
Consequently, it is highly desirable to elucidate PPIs in cancer and construct
flexible graphical models that can identify multiple types and ranges of
dependencies. Modern data collection methods have allowed systematic
assessment of multiple proteins simultaneously on the same tumor samples,
often referred to as high-throughput proteomics (Baladandayuthapani et al.,
2014). However, the resulting data are typically not normally distributed even
after extensive preprocessing and data transformations (e.g. logarithmic). As
an illustration, Figure 1 shows the level of non-normality in protein
expression data after logarithmic transformation for two cancers: lung
adenocarcinoma (LUAD) and ovarian cancer (OV) samples from The Cancer Genome
Atlas (TCGA) (Weinstein et al., 2013) that are used as case-studies in this
paper. Panels (A) and (B) display the empirical density and the quantile-
quantile (q-q) plots of four exemplar proteins: Akt and PTEN for LUAD, and
E-Cadherin and Rb for OV. Both the empirical distributions and q-q plots
demonstrate deviations from normal distribution with heavier tails as shown in
Panels (A) and (B). The level of non-normality is further quantified using the
H-score, defined as $H(\bm{Y})=2\Phi[\log\\{1-\text{pval}(\bm{Y})\\}]$, where
$\Phi$ is the cumulative distribution function of the standard normal
distribution, and $\text{pval}(\bm{Y})$ is the p-value of the Kolmogorov-
Smirnov test for the normality of $\bm{Y}$ (Chakraborty et al., 2021). The
H-score is bounded between zero and one, and a higher H-score implies
increased departure from normality. The H-scores for all four proteins are
$>0.999$, consistent with the conclusions from the empirical and q-q plots.
Panel (C) shows the H-score across all the proteins in our datasets,
indicating a high degree of non-normality across both cancers.

Another axes of complexity that arises in cancer research is tumor
heterogeneity. It is now well-established that tumors are heterogeneous with
distinct proteomic aberrations even for the same type of cancer across
different patients (Janku, 2014). Accumulating evidence suggests that
considering tumor heterogeneity, in general, and specifically at the level of
PPIs can enhance our understanding of tumorigenesis and the development of
anti-cancer treatments (Cheng et al., 2020). Specifically, tumor heterogeneity
differentially impacts the PPIs across different patients and results in
varied treatment responses (Cheng et al., 2020). Hence, incorporating patient-
specific information, i.e., accounting for tumor heterogeneity, could provide
valuable clues to identify PPIs disrupted during carcinogenesis. In summary,
constructing PPI networks poses two main statistical challenges
simultaneously: (i) coherently accounting for non-normality in proteomic
networks, and (ii) incorporating heterogeneous patient-specific information in
graphical modeling.

Figure 1: Non-normality levels of protein expression data from lung
adenocarcinoma (LUAD) and ovarian cancer (OV). The empirical density plots
from actual data (black) and the normal distribution (blue) for the expression
of four proteins with the H-score are shown in the upper-right of each protein
in Panel (A). Panel (B) illustrates the expression of four proteins in LUAD
(Akt and PTEN) and OV (E-Cadherin and Rb) with qq-plots. Panel (C) shows the
violin plot of H-scores across all proteins in LUAD and OV cancers. The
H-score is bounded between zero and one, and a higher H-score implies a higher
level of non-normality (see manuscript for description).

Existing methods and modeling background. Most existing methods address the
aforementioned challenges separately, i.e., either accommodating non-normality
without accounting for the sample-specific information (e.g. Pitt et al.,
2006; Dobra and Lenkoski, 2011) or requiring normality when incorporating
patient-specific information (Ni et al., 2022). To accommodate non-normality,
existing approaches transform the original variables into normal variables
either via deterministic functions (e.g. Dobra and Lenkoski, 2011; Liu et al.,
2012; Chung et al., 2022) or via random transformations (e.g. Finegold and
Drton, 2011, 2014). For instance, Bhadra et al. (2018) use a Gaussian scale
mixture technique that generalizes the $t$-distribution and introduce a new
graph characterization for undirected graphs. Chakraborty et al. (2021)
further generalize this concept to characterize chain graphs with both
directed and undirected edges. However, all existing models mentioned above
assume a common graph across all patients and fail to incorporate the subject-
specific information.

More recently, several studies incorporate the subject-specific information
under explicit Gaussian assumptions. Multiple Gaussian graphical models were
first proposed to estimate graphs that vary across heterogeneous sub-
populations (e.g. Peng et al., 2009; Danaher et al., 2014; Peterson et al.,
2015). Ni et al. (2019) introduce a more general framework called “graphical
regression” that constructs covariate-dependent graphs through a regression
model and incorporates both continuous and discrete covariates, in directed as
well as undirected settings (Ni et al., 2022). Similarly, Zhang and Li (2022)
provide a penalized procedure to estimate undirected graphs through Gaussian
graphical regression and introduce continuous covariates in both the mean and
the covariance structures. However, all these models are developed under the
normality assumption for inferential and computational reasons. To the best of
our knowledge, no existing method incorporates subject-specific information
under non-Gaussian settings and motivates development of new methodology. We
summarize six important and relevant models mentioned above in Table 1 and
compare these models in four different aspects.

Table 1: Comparison of existing methods and our proposed rBGR method across four different properties. Method | Uncertainty | Undirected | Sample- | Non-Normality
---|---|---|---|---
Quantification | specific
GGMx (Ni et al., 2022) | ✓ | ✓ | ✓ | ✗
RegGMM (Zhang and Li, 2022) | ✗ | ✓ | ✓ | ✗
GSM (Bhadra et al., 2018) | ✓ | ✓ | ✗ | ✓
BGR (Ni et al., 2019) | ✓ | ✗ | ✓ | ✗
RCGM (Chakraborty et al., 2021) | ✓ | ✓ | ✗ | ✓
rBGR (the proposed) | ✓ | ✓ | ✓ | ✓

To address these challenges simultaneously, we develop a unified and flexible
modeling strategy called robust Bayesian graphical regression (rBGR), which
allows construction of subject-specific graphical models for non-normally
distributed continuous data. rBGR makes three main contributions:

  * (a)

Robust framework to build subject-specific graphs for non-normal data. rBGR
robustifies the normal assumption via random transformation and incorporates
covariates employing graphical regression strategies. By accommodating non-
normality via random scale transformations, we obtain a Gaussian scale
mixture, which presumes an underlying latent Gaussian variable, allows
explicit incorporation of covariates in the precision matrix (Section 2.2),
and admits efficient posterior sampling procedures (Section 4).

  * (b)

New characterization of dependency structures for non-normal graphical models.
The introduction of the random marginal transformations engenders a new type
of edge characterization of the conditional dependence for non-normal data,
called conditional sign independence with covariates (CSIx, Proposition 2).
CSIx is a generalization of the notion of conditional sign independence (CSI),
introduced by Bhadra et al. (2018), which explicitly characterizes the sign
dependence between two variables that holds for a much broader class of models
than Gaussian graphical models. We demonstrate via multiple simulations that
rBGR can accurately recover dependency structures under different levels of
non-normality and against competing graphical regression approaches that
assume normality (Section 5).

  * (c)

Deciphering impact of immunogenic heterogeneity in proteomic networks. We use
rBGR to assess proteomic networks across two cancers, lung and ovarian, to
systematically investigate the effects of the inherent immunogenic
heterogeneity within tumors. Specifically, we quantify immune cell abundance
across tumors and build PPI networks that vary across different immune cell
abundance. Our analyses reveal several important hub proteins and PPIs that
are impacted by the immune cell abundance; some corroborate existing
biological knowledge whereas others are novel associations (Section 6).

The rest of the paper is organized as follows: we introduce rBGR models and
characterization in Section 2. Section 3 focuses on priors and estimation, and
Section 4 delineates the posterior inference via Gibbs sampling. In Section 5,
we conduct a series of simulations to evaluate the operating characteristics
of rBGR against competing approaches. Section 6 provides a detailed analysis
of the TCGA dataset, results, biological interpretations, and implications.
The paper concludes by discussing implications of the findings, limitations,
and future directions in Section 7. A general purpose R package and datasets
used in this paper for constructing PPI networks is also provided at
https://github.com/bayesrx/rBGR.

##  2 Robust Bayesian Graphical Regression (rBGR)

We start by reviewing the Gaussian graphical regression (Section 2.1), which
is a special case of rBGR under the normality assumption, and then generalize
it to the robust case by random transformations (Section 2.2). Subsequently,
the introduction of the random scale changes the interpretation of graph and
motivates a new edge characterization (Section 2.3).

###  2.1 Gaussian Graphical Regression

Consider $p$-dimensional random vectors
$\bm{Y}_{i}=(Y_{i1},\ldots,Y_{ip})^{\sf T}\in\mathbb{R}^{p}$ as (continuous)
responses with $q$-dimensional random vectors of
$\bm{X}_{i}=(X_{i1},\ldots,X_{iq})^{\sf T}\in\mathbb{R}^{q}$ as covariates for
subject $i=1,\ldots,n$. A subject-specific PPI network from proteomics data
$\bm{Y}_{i}$ is constructed to vary based on the immune cell abundance
$\bm{X}_{i}$ (Section 6). Let $G_{i}=(V,E_{i})$ be an undirected graph over
$p$ nodes, where $V=\\{1,\ldots,p\\}$ is the set of nodes representing
$\bm{Y}_{i}$ and $E_{i}\subset V\times V$ is the set of undirected edges in
the graph for subject $i$. An undirected edge exists between nodes $j$ and $k$
if $\\{j,k\\}\in E_{i}$. Under Gaussian assumption, given the covariates
$\bm{X}_{i}$, suppose $\bm{Y}_{i}$ follows a multivariate normal distribution,

$\displaystyle\bm{Y}_{i}\mid\bm{X}_{i}\sim\bm{N}_{p}({\bf
0},\widetilde{\bm{\Omega}}^{-1}(\bm{X}_{i})),\textrm{ for }i=1,\ldots,n,$ (1)

where
$\widetilde{\bm{\Omega}}(\bm{X}_{i})=\\{\widetilde{\omega}^{j,k}(\bm{X}_{i})\\}_{p\times
p},j,k\in V$ is a functional precision matrix (of covariates) with each
element $\widetilde{\omega}^{j,k}(\bm{X}_{i})$ as a function that depends on
$\bm{X}_{i}$. The functional precision matrix characterizes the graph $G_{i}$
through zero precision elements. Specifically, a zero element of the precision
matrix represents a missing edge in the graph, e.g., for the case of scalar
precision matrix,
$\widetilde{\omega}^{j,k}(\bm{X}_{i})=\widetilde{\omega}^{j,k}$, zero
precision implies conditional independence (CI) and an missing edge in the
graph of CI under Gaussianity (Lauritzen, 1996). For the functional precision
matrix, Ni et al. (2022) introduced covariate-dependent graphs in $G$ and
generalized the concept of CI to CI with covariates (CIx, henceforth). In
essence, given a covariate $\bm{X}_{i}$, the zero precision of
$\widetilde{\omega}^{j,k}(\bm{X}_{i})=0$ implies an missing edge of CIx
between nodes $j$ and $k$. Contrarily, when the functional precision is non-
zero $\widetilde{\omega}^{j,k}(\bm{X}_{i})\neq 0$, $Y_{j}$ and $Y_{k}$ are
conditional dependent with covariates (CDx, henceforth) and an edge exists
between nodes $j$ and $k$ given the covariate $\bm{X}_{i}$. By modeling the
functional precision matrix, CIx defines covariate-specific graphs that vary
based on different covariates.

###  2.2 Robust Graphical Regression via Random Transformation

In practice, normal assumption does not always hold (as shown in Figure 1).
Violation of the normal assumption results in the failure of modeling graphs
through normal precision matrices and motivates new modeling strategies
(Finegold and Drton, 2011; Bhadra et al., 2018). In this paper, we adapt the
random transformation approach (Bhadra et al., 2018) that allows for various
non-normal distributions with different tail behaviors. We focus on continuous
distributions with heavy tails as observed in our motivating data. To this
end, let $0<d_{j}<\infty$ for $j=1,\ldots,p$ be independent positive random
scales and have distribution as $d_{j}\sim p_{j}$ with $\int dp(d_{j})<\infty$
almost surely. Let $\mathbf{D}_{i}=\text{diag}(1/d_{i1},\ldots,1/d_{ip})$ be a
diagonal matrix for subject $i$. Given random scales $d_{ij},j=1,\ldots,p$ and
the covariates $\bm{X}_{i}$, we assume the distribution of
$\mathbf{D}_{i}\bm{Y}_{i}$ conditional on $\mathbf{D}_{i}$ and $\bm{X}_{i}$
follows a multivariate distribution,

$\displaystyle\mathbf{D}_{i}\bm{Y}_{i}=\left[\frac{Y_{i1}}{d_{i1}},\ldots,\frac{Y_{ip}}{d_{ip}}\right]^{\sf
T}\sim\bm{N}_{p}(\mathbf{0},\bm{\Omega}^{-1}(\bm{X}_{i})),\textrm{ for
}i=1,\ldots,n,$ (2)

where $\bm{\Omega}(\bm{X}_{i})=\\{\omega^{j,k}(\bm{X}_{i})\\}_{p\times
p},j,k\in V$ is the functional precision matrix that characterizes the graph
with the covariates $\bm{X}_{i}$.
The model in (2) generalizes several existing approaches: (i) Equation (1) is
a special case of Equation (2) with $d_{ij}$ as a degenerated distribution of
a point mass at one; (ii) when $d_{1}=\ldots=d_{p}=\tau$ with $\tau^{2}$
following an inverse gamma distribution, Equation (2) reduced to a
multivariate t-distribution on $\bm{Y}$ as used by Finegold and Drton (2014),
and (iii) for general $d_{ij}$, (2) establishes a rich family of Gaussian
scale mixtures for the marginal distribution of $Y_{j}$ with the density
$p(Y_{j})=\int(2\pi
d_{j})^{-1/2}\exp\\{-y_{j}^{2}/(2d_{j})\\}\textrm{d}p(d_{j})$.

The introduction of random scales in Equation (2) allows us to construct
various marginal distribution of $Y_{j}$ with different tail behaviors.
Specifically, by matching tail behaviors of random scales and the target
distribution, random scales allow us to model different marginal distributions
that the data might exhibit. For example, letting $Y_{j}$ decay polynomially,
$Y_{j}/d_{j}$ follows a normal distribution if the random scale $d_{j}$ also
has a polynomial tail. In Figure 2, Panel (A) shows that the target
distribution of $Y$ with a polynomial tail deviates from the normal
distribution but with the introduction of random scales the distribution of
$Y/d$ is normally distributed. Similar idea can be used for target
distribution with other tail behaviours e.g. exponential tails. While the
random scales robustify the model to accommodate non-normality, the resulting
functional precision matrix $\bm{\Omega}(\bm{X}_{i})$, however, requires
careful characterization and interpretation.

Figure 2: Robustification of non-normal distribution with random scales and
the visualization of CSIx and CSDx. Panel (A) shows the qq-plot to illustrate
that random scale $d$ accommodates the non-normal distribution $Y$ with $Y/d$
following the normal distribution. Panel (B) demonstrates CSIx (Case (i) and
(ii)) and CSDx (Case (iii) and (iv)) of $Y_{1}$ and $Y_{2}$ with the partial
correlation $\omega^{1,2}(X_{i})=X_{i}$ conditioning on $Y_{3}$. Cases (i) and
(ii) represent two examples of CSIx with zero precision of $X_{i}=0$ given
$Y_{3}=1$ and $0$. Cases (iii) and (iv) demonstrate the cases of CSDx with
non-zero precision of $X_{i}=0.7$ given $Y_{3}=1$ and $0$. Panel (B) is
centered on the values between $[-10,10]$. Panel (C) shows the nested
relationship between CSIx and CIx (top) and CSDx and CDx (bottom). See more
details in Section 2.3.

###  2.3 Characterization of Functional Precision Matrix

The functional precision matrix in (2) determines the graphical dependence as
a function of covariates, but the random (marginal) scale changes the standard
conditional independence interpretations in the resulting precision matrix,
which requires a new characterization. Bhadra et al. (2018) introduced the
concept of conditional sign independence (CSI) in non-normal graphs that is
defined as follows. Consider two protein expression of interest as random
variables $Y_{1}$ and $Y_{2}$ with the expression data from the rest of the
proteins denoted by a random vector $\bm{Y}_{3}$. Given $\bm{Y}_{3}$, $Y_{1}$
and $Y_{2}$ are conditional sign independence (CSI) if
$\mathcal{P}(Y_{2}>0\mid Y_{1},\bm{Y}_{3})=\mathcal{P}(Y_{2}>0\mid\bm{Y}_{3})$
and $\mathcal{P}(Y_{1}>0\mid
Y_{2},\bm{Y}_{3})=\mathcal{P}(Y_{1}>0\mid\bm{Y}_{3})$. Otherwise, $Y_{1}$ and
$Y_{2}$ are conditional sign dependent (CSD) given $\bm{Y}_{3}$. The CSI of
$Y_{1}$ and $Y_{2}$ implies that the information of $Y_{1}$ does not affect
the sign of $Y_{2}$ given $\bm{Y}_{3}$. That is, conditioning on the rest of
the protein expression data $\bm{Y}_{3}$, the probability of over- or under-
expression for $Y_{2}$ is independent of the expression level of $Y_{1}$.
Under the multivariate distribution of (2) with a constant precision matrix
$\bm{\Omega}(\bm{X}_{i})=\bm{\Omega}$, zero precision of $\omega^{j,k}=0$ and
the CSI of $Y_{j}$ and $Y_{k}$ given the rest are equivalent, which can be
represented by a missing edge between nodes $j$ and $k$ in an undirected graph
(Bhadra et al., 2018; Chakraborty et al., 2021).

In this paper, we generalize the concept of CSI to incorporate covariates and
consider subject-specific CSI of two random variables given all the other
random variables and a realization of covariates $\bm{X}_{i}$, as formalized
in the following proposition:

######  Proposition 1 (Conditional Sign Independence with Covariate (CSIx))

Given random scales $\mathbf{D}_{i}=\text{diag}(1/d_{i1},\ldots,1/d_{ip})$ and
the covariates $\bm{X}_{i}$, consider the conditional distribution of
$\mathbf{D}_{i}\bm{Y}_{i}$ as Equation (2) with functional precision matrix
$\bm{\Omega}(\bm{X}_{i})$. If $\omega^{j,k}(\bm{X}_{i})=0$, then $Y_{ij}$ and
$Y_{ik}$ are CSI. Otherwise, when $\omega^{j,k}(\bm{X}_{i})\neq 0$, then
$Y_{ij}$ and $Y_{ik}$ are CSD.

The proof of Proposition 1 follows the fact that $\omega^{j,k}(\bm{X}_{i})=0$
implies the CSI of $Y_{ij}$ and $Y_{ik}$ given $\bm{X}_{i}$, and we call
$Y_{ij}$ and $Y_{ik}$ are conditional sign independence with covariates
$\bm{X}_{i}$ to highlight the role of the covariates in the graph. Otherwise,
$Y_{ij}$ and $Y_{ik}$ are called CSDx. See Supplementary Material Section S1
for more details.

An illustrative example. We use a simple low-dimensional example to visually
demonstrate and interpret CSIx and CSDx. Following Proposition 1, we show two
examples with a general functional precision matrix $\bm{\Omega}(\bm{X}_{i})$.
Consider a trivariate distribution of (2) with unit diagonal elements and
$\omega^{1,2}(X_{i})=X_{i}$. We illustrate two scenarios shown in Panel (B) of
Figure 2:

  * •

When $X_{i}=0$, we obtain the CSIx of $Y_{1}$ and $Y_{2}$ given two different
values of $Y_{3}=0$ (Case (i)) and $1$ (Case (ii)).

  * •

When $X_{i}=0.7$, $Y_{1}$ and $Y_{2}$ are CSDx and we observe that the
distribution of the sign of $Y_{2}$ varies based on the value of $Y_{1}$ (see
Case (iii) and (iv)). Specifically, as $Y_{1}$ increases, $Y_{2}$ tends to be
negative.

By modeling the functional precision matrix, we can build covariate-specific
precision matrix that depends on the different realization of the covariates
$\bm{X}_{i}$. Consequently, we can construct a graph of CSI corresponding to
the precision matrix and the covariates.

We can now conceptually compare models (1) and (2). Both models incorporate
the covariates in the functional precision matrix, which characterize the
covariate-specific graph. However, the interpretation of the graph differs.
The graph from model (2) encodes CSIx whereas the graph from model (1) encodes
CIx. We further visualize the relationship between CSIx and CIx in Panel (C)
of Figure 2 and summarize as follows:

  * •

For $\omega=0$, CSIx is a weaker condition than CIx since CSIx only considers
the sign while CIx depends on both the sign and the magnitude.

  * •

When $\omega\neq 0$, CSDx is a stronger condition than CDx since CDx allows
either magnitude or the sign to be dependent while CSDx only focuses on the
sign.

##  3 Priors and Estimation

The functional precision matrix $\bm{\Omega}(\bm{X}_{i})$ lives in a high-
dimensional space. For example, the PPI network for ovarian cancer from our
application considers $\mathcal{O}(n\times p\times(p-1)/2)=197,620$ possible
edges – which makes joint estimation difficult if not untenable, especially
since variability across each subject $i$ is allowed. Hence, we employ a
neighborhood selection procedure (Meinshausen and Bühlmann, 2006) to estimate
the functional precision matrix that has been used in several graphical
modeling approaches (e.g. Ni et al., 2019; Zhang and Li, 2022). This procedure
offers three main benefits: (i) tractable estimation, (ii) reduced computation
burden, and (iii) flexible prior elicitation. Specifically, we regresses one
node $Y_{j}$ on the rest nodes $Y_{k},\;k\neq j$ and build the graph based on
zero coefficients (Section 3.1 and 3.2). This use of neighborhood selection,
which employs conditional estimation as opposed to joint estimation in (2)
effectively reduce the number of parameters (edges) to $\mathcal{O}(q\times
p\times(p-1)/2)=3,280$ – according a 100-fold reduction. Additionally, the
effective number of edges can be further reduced by different model
specification such as a thresholding mechanism (Section 3.3) and different
priors such as spike-and-slab (Section 3.4).

###  3.1 Regression-based Approach for Functional Precision Matrix Estimation

The rBGR model leverages a regression-based framework on model (2) to relate
the regression coefficients and precision matrix. Given random scales
$\bm{D}_{i}$, we regress one variable on all other variables and relates the
partial correlation with regression coefficients. Zero coefficients is then
equivalent to zero partial correlations (Meinshausen and Bühlmann, 2006).
Specifically, we construct node-specific regressions as:

$\displaystyle\frac{Y_{ij}}{d_{ij}}=\sum_{k\neq j}^{p}\beta_{j,k}({\bf
X_{i}})\frac{Y_{ik}}{d_{ik}}+\epsilon_{ij},$ (3)

where $\epsilon_{ij}\sim N(0,1/\omega^{j,j}(\bm{X}_{i}))$ and the functional
coefficient
$\beta_{j,k}({\bm{X}_{i}})=-\frac{\omega^{j,k}(\bm{X}_{i})}{\omega^{j,j}(\bm{X}_{i})}$.
Under this specification, $\beta_{j,k}(\bm{X}_{i})=0$ if and only if
$\omega^{j,k}(\bm{X}_{i})=0$, which enables the functional coefficients to
characterize the covariate-specific graphs. However, the interpretation of the
coefficients changes from the standard Gaussian graphical models due to the
introduction of the random scales, which is detailed in the next subsection.

###  3.2 Graph Construction through Regression Coefficients

We build graphs with a missing edge between nodes $j$ and $k$ when $Y_{j}$ and
$Y_{k}$ are CSIx given the remaining variables and the covariates
$\bm{X}_{i}$. Consider $\bm{Y}_{i}$ and $\bm{X}_{i}$ with the regression (3).
We call $\beta_{j,k}(\bm{X}_{i})$ the conditional sign independence function
(CSIF) because zero CSIF $\beta_{j,k}({\bf X}_{i})=0$ implies that $Y_{j}$ and
$Y_{k}$ are CSIx given all the other nodes $\bm{Y}_{-\\{j,k\\}}$ and
covariates $\bm{X}_{i}$, as formally characterized in the following
proposition.

###### Proposition 2

Consider model (3). If $\beta_{j,k}({\bf X}_{i})=0$, then
$\mathcal{P}(Y_{j}>0\mid Y_{k},\bm{Y}_{-\\{j,k\\}},{\bf
X}_{i})=\mathcal{P}(Y_{j}>0\mid\bm{Y}_{-\\{j,k\\}},{\bf X}_{i})$ and
$\mathcal{P}(Y_{k}>0\mid Y_{j},\bm{Y}_{-\\{j,k\\}},{\bf
X}_{i})=\mathcal{P}(Y_{k}>0\mid\bm{Y}_{-\\{j,k\\}},{\bf X}_{i})$.

We sketch the proof and leave the details in Supplementary Section S1. The
proof follows from the fact that the CSIF
$\beta_{j,k}(\bm{X}_{i})=-\frac{\omega^{j,k}(\bm{X}_{i})}{\omega^{j,j}(\bm{X}_{i})}$
is related to the partial correlation, and a zero partial correlation is
equivalent to a zero precision of $\omega^{j,k}(\bm{X})=0$, which ensures the
CSIx between $Y_{j}$ and $Y_{k}$ (see the example in Section 2.3). Therefore,
zero CSIF indicates the CSIx between $Y_{j}$ and $Y_{k}$ given the remaining
response variables $\bm{Y}_{-\\{j,k\\}}$ and covariates $\bm{X}_{i}$. In this
paper, we further assume $\omega^{j,j}$ in CSIF to be scalar,
$\beta_{j,k}(\bm{X}_{i})=-\frac{\omega^{j,k}(\bm{X}_{i})}{\omega^{j,j}}$, for
ease of computation. Note that our main interest is edge selection and the
CSIF is zero if and only $\omega^{j,k}(\bm{X}_{i})=0$, which is unrelated to
$\omega^{j,j}$.

###  3.3 Modeling the Conditional Sign Independence Function

Proposition 2 transforms the problem of robust graph construction to a more
tractable regression coefficient selection (i.e., selecting which part of CSIF
is exactly zero). Therefore, modeling the CSIF is crucial to the graph
estimation. To this end, we parameterize the CSIF as a product of two
components:

$\displaystyle\beta_{j,k}(\bm{X}_{i})=\underbrace{\theta_{j,k}(\bm{X}_{i})}_{\textrm{Covariate
function}}\underbrace{\mathbb{I}(\mid\theta_{j,k}(\bm{X}_{i})\mid>t_{j})}_{\textrm{Thresholding
function}}.$ (4)

We elaborate the role and justification of each component below.

Covariate functions [$\theta_{\bullet}(\bullet)$]. For simplicity, we consider
only the linear effects of covariates $\bm{X}_{i}$,
$\theta_{j,k}(\bm{X}_{i})=\sum_{h=1}^{q}\alpha_{j,k,h}X_{ih}$, where
$\alpha_{j,k,h}$ represents the coefficients for the $h$-th covariate. The
covariate function allows similar edge sets for individuals with a similar
level of $\bm{X}_{i}$ and varies the graph thus borrowing strength. If
desired, it is relatively straightforward to extend it to nonlinear effects
with e.g., using basis expansion techniques such as splines.

Thresholding functions
[$\mathbb{I}(\mid\theta_{\bullet}(\bm{X})\mid>t_{\bullet})$]. The edge
thresholding mechanism is desired to achieve sparse graphs in rBGR due to the
large number of parameters and according multiplicity adjustments. For
example, the ovarian PPI network in our application requires $qp(p-1)/2=3,280$
parameters and results in a dense graph with inefficient inference. To solve
the problem, we truncate edges with small magnitudes with an indicator
function $\mathbb{I}(\mid\theta_{j,k}(\bm{X})\mid>t_{j})$, where $t_{j}$ is
the threshold parameter specific to the node $j$. An edge is shrunk to zero
and removed when the magnitude is smaller than the threshold parameter,
resulting in a sparse graph. One might consider threshold parameter as
$t_{j,k}$. However, $t_{j,k}$ is not fully identifiable when
$\alpha_{j,k,h}=0$ for all $h=1,\ldots,q$ since when
$\theta_{j,k}(\bm{X}_{i})=0$, the value of $t_{j,k}$ can be arbitrary. To
alleviate the problem, we assume $t_{j,k}=t_{j}$ to improve the
identifiability as long as one of of $\theta_{j,k}\neq 0$.

###  3.4 Prior Specification

To complete the model specification, rBGR contains three parameters: (a)
random scales $d_{j}$, (b) threshold parameter $t_{j}$, and (c) covariate
coefficients $\alpha_{j,k,h}$. Specifically, we assign priors as follows:

$\displaystyle
d_{j}\sim(1-\pi_{j})\delta_{1}(d_{j})+\pi_{j}p_{j}(d_{j});\thickspace
t_{j}\sim\textrm{Unif}(0,t_{\max})$
$\displaystyle\alpha_{j,k,h}\sim\gamma_{j,k,h}N(0,\nu_{j,k,h})+(1-\gamma_{j,k,h})N(0,v_{0}\nu_{j,k,h}),$

where $v_{0}$ and $t_{\max}$ are pre-specified hyperparameters, $\pi_{j}$
models the degree of non-normality with beta prior as $\pi_{j}\sim
Beta(a_{\pi},b_{\pi})$, $p_{j}$ is a function to accommodate the non-
normality, $\gamma_{j,k,h}$ is a binary variable with Bernoulli prior as
$\gamma_{j,k,h}\sim Ber(\rho_{j})$, and $\nu_{j,k,h}$ decides the variance of
$\alpha_{j,k,h}$ with a inverse Gamma prior of $\nu_{j,k,h}\sim
InvGa(a_{\nu},b_{\nu})$. Specifically, when $d_{j}=1$, $Y_{j}$ is normally
distributed. When $d_{j}\sim p_{j}$, $Y_{j}$ follows a non-normal
distribution. We match tail behavior of $p_{j}$ and the marginal distribution
of $Y_{j}$, and allow each marginal distribution $Y_{j}$ to have different
level of non-normality by specific $\pi_{j}$. For the current model, we focus
on the $Y_{j}$ with polynomial decay as illustrated by the motivating data in
Figure 1 and assign a inverse Gamma distribution on $p_{j}(d_{j}^{2})\sim
InvGa(a_{d},b_{d})$. For threshold parameter $t_{j}$, we assign a uniform
prior on $t_{j}$ to model the thresholding mechanism and control the graph
sparsity. Intuitively, when $t_{j}\rightarrow 0$, no edge is truncated with a
fully connected graph. When $t_{j}\rightarrow\infty$, all edges are shrunk to
zero with all nodes disconnected. For covariate coefficients $\alpha_{j,k,h}$,
we assign a spike-and-slab prior to achieve the covariate sparsity with a
small $v_{0}$ because not all covariates necessarily contribute to the varying
structure of our graph.

##  4 Posterior Inference

Gibbs sampler. In this Section, we introduce an efficient Gibbs sampler for
the proposed rBGR model. Instead of the Metropolis-Hastings algorithm, we
implement a Gibbs sampler except for the random scales, which largely improves
the computation and convergence compared to Ni et al. (2019). Recently, Li et
al. (2023) derived a closed-form of the conditional distributions for Gibbs
sampler by formulating the thresholded coefficients as mixture distributions.
Specifically, if we view the distribution with one component as a special case
of mixture distribution, the mixture distribution from the thresholded
coefficient then can achieve conjugacy. We derive the full conditional
distributions for covariate coefficients $\alpha_{j,k,h}$ and the threshold
parameter $t_{j}$, which are the mixture of truncated normal and the mixture
of uniform distribution, respectively. By assigning normal priors on covariate
coefficients and a uniform prior on threshold parameter, we obtain conjugacy
on all thresholded parameters. We further use the parameter expansion
technique (Ni et al., 2019; Scheipl et al., 2012) on covariate coefficients to
improve the mixing of MCMC. We implement the Metropolis-Hasting algorithm for
the random scales. The full details of the full conditional distributions and
sampling algorithm are provided in Supplementary Materials (Section S2.2).

Covariate and edge selection. The estimated coefficients from rBGR of (3) do
not guarantee symmetry required for undirected graphs. Also, due to the
introduction of random scales with CSIx characterization, we focus on the sign
of the edge. We describe algorithms to symmetrize the estimated covariate
coefficients $\widehat{\alpha}_{j,k,h}$ and the sign of graph edges of
$\widehat{\beta}_{j,k}(\bm{X}_{i})$. For covariate coefficients, we compare
the posterior inclusion probability (PIP) of asymmetric coefficients from two
directions ($\widehat{\alpha}_{j,k,h}$ and $\widehat{\alpha}_{k,j,h}$) and
assign the symmetrized coefficients as the directed coefficient with a smaller
PIP. Given a cutoff $c_{0}$, the rule above requires both directions of
coefficients to have PIPs bigger than $c_{0}$ encouraging a parsimonious
network. For the edge, we symmetrize the edge based on the edge posterior
probability (ePP). Specifically, we symmetrize ePP by taking the maximum of
two asymmetric ePP. Given a cutoff $c_{1}$, we call an undirected edge if at
least one of the directed ePPs is bigger than $c_{1}$. We then decide the sign
of the edge by comparing the posterior probability of positive and negative
for the chosen direction. The full details are provided in Supplementary
Materials (Section S2.3).

##  5 Simulation Studies

We empirically demonstrate the performance of rBGR under a variety of non-
normal contaminations and against other competing models in terms of edge and
covariate selection. To the best of our knowledge, no other existing method
estimates covariate-specific graphs for non-normal data. Therefore, we compare
rBGR to two models that estimate the covariate-specific graph without
addressing the violation of normality assumption. Specifically, we consider
Bayesian graphical regression (BGR) (Ni et al., 2019) and the Gaussian
graphical model regression (RegGMM) (Zhang and Li, 2022) representative of a
fully Bayesian and a frequentist penalization-based models for the covariate-
specific graph under normal assumption, respectively. For RegGMM, we run the
algorithm with various tuning parameters to obtain the probability of the
signs of edges and covariate coefficients and select the optimal tuning
parameter by cross-validation using their default algorithm. For rBGR and BGR,
we symmetrize the graph mentioned in Section 4 and set $c_{0}=c_{1}=0.5$. We
run $10,000$ and $30,000$ iterations and discard the first $90\%$ iterations
for rBGR and BGR, respectively.

Data generating mechanism. We generate the observed non-normal data by
multiplying the random scale to the latent normal data
$\bm{Y}_{i}^{*}=\left[Y_{i1}^{*},\ldots,Y_{ip}^{*}\right]^{\sf T}$ that
follows an multivariate normal distribution with a functional precision matrix
that represents the undirected graph. Specifically, we generate the covariates
$\bm{X}_{i}\overset{iid}{\sim}U(-1,1)$ and latent data
$\bm{Y}_{i}^{*}\overset{iid}{\sim}\bm{N}_{p}(\mathbf{0},\bm{\Omega}^{-1}(\bm{X}_{i}))$
with the true precision matrix $\bm{\Omega}(\bm{X}_{i})$. For
$\bm{\Omega}(\bm{X}_{i})$, we assign unit diagonal elements and randomly pick
$2\%$ of the off-diagonal to be non-zero. We let the non-zero precision depend
on the covariates linearly and truncate the precision with a magnitude smaller
than $0.15$. We obtain the random scales from a mixture distribution of the
point mass at one and a inverse gamma distribution and assign three different
levels non-normal contamination: $\pi\in\\{0,0.5,0.8\\}$. We multiply the
random scales to generate the observed data of
$[Y_{i1},\ldots,Y_{ip}]=[Y_{i1}^{*}d_{i1},\ldots,Y_{ip}^{*}d_{ip}]$. For all
simulations, we set the sample size and the dimensions of $\bm{Y}_{i}$ and
$\bm{X}_{i}$ as $(n,p,q)=(250,50,3)$ based on our real data case-studies. We
show the results for $50$ independent replicates.

Performance metrics. We evaluate the graph recovery through the edge and
covariate selection. For covariate selection, we report the true positive rate
(TPR), true false rate (TFR), and Matthew’s correlation coefficient (MCC) with
the cut-off for PIP at $c_{0}=0.5$. We also report the area under the ROC
curve (AUC) and partial area under ROC curve (pAUC) between specificity
ranging from 0.8 to 1. For edge selection, we use AUC and three metrics of
TPR, TNR and MCC with the cut-off for ePP at $c_{1}=0.5$. We further
investigate the sign consistency by examining the agreement between the
posterior probability for the signs of CSIF
$\textrm{sgn}(\widehat{\beta}_{j,k}(\bm{X}_{i}))$ and the true signs of
$\textrm{sgn}(\beta_{j,k}(\bm{X}_{i}))$. Specifically, we exclude the zero
CISF and focus on the subset of the data with both true and estimated non-zero
CSIF to restrict the problem as two-class classification (positive versus
negative). We assess the sign consistency by MCC (referred to as sign-MCC).

Simulation results. Panel (A) of Figure 3 shows the simulation results for
covariate selection. We observe that rBGR outperforms BGR and RegGMM across
all non-normality levels, as indicated by higher MCC and AUC. The difference
of MCC and AUC between rBGR and the other competing methods increases when the
non-normality contamination level increases, which is expected. For TNR, rBGR
performs slightly worse than BGR but better than RegGMM across all non-
normality levels. However, all three methods select correct covariates
($>93\%$) with small difference ($<5\%$) in terms of TNR. For TPR, rBGR
outperforms BGR under all levels of non-normality and the advantage of rBGR
becomes more prominent as the non-normality increases. Compared to RegGMM,
rBGR’s performance is comparable under normal distribution in TPR, but rBGR is
preferred when the level of non-normality increases. Overall, modeling the
non-normality from random scales in rBGR is favored compared to models without
random scales in terms of covariate selection.

We show the graph recovery for the edge selection in Panel (B) of Figure 3.
For edge selection, rBGR outperforms BGR and RegGMM in AUC, and the advantage
of rBGR increases with a larger discrepancy between rBGR and the competing
methods when the non-normality level increases. For MCC, rBGR outperforms
RegGMM under all levels of non-normality, but is slightly inferior than BGR
under the normal distribution. However, rBGR is favored when the non-normality
level increases. For TPR, rBGR is better than BGR under all levels of non-
normality, and slightly worse than RegGMM under normal assumption. However,
when non-normality increases, rBGR starts to surpass the RegGMM. Both TNR and
sign-MCC show excellent selection performance ($>95\%$) for all three methods,
with minimal differences ($<5\%$) across the three non-normality levels.

In summary, modeling the non-normality through random scales in rBGR result in
equivalent (under normal distribution) or better performances in all metric
for edge selection compared to the other methods.

Figure 3: Simulation results: graph recovery for BGR (red), rBGR (green) and
RegGMM (blue) under different levels of non-normality in terms of (A)
covariates selection (top row) and (B) edge selection (bottom two rows). Panel
(A) shows the covariate selection through four metrics (from left to right:
TPR, TNR, MCC and AUC) that are measured under three different levels of non-
normality. Panel (B) demonstrates the edge selection by four criteria (from
upper left to lower right: TPR, TNR, MCC, AUC) and the sign consistency by
sign-MCC (lower left) for non-zero edges. All values for TPR, TNR and MCC are
measured at a cut-off at $c_{0}=c_{1}=0.5$.

Additional simulations and model evaluations. We provide additional simulation
of data generating mechanism and model evaluation results for (i) convergence
of the algorithm, and (ii) different cut-off of $c_{0}$ and $c_{1}$
controlling for false discovery rates – which are summarized in Supplementary
Material Section S3. Overall, rBGR generates equivalent or better performances
compared to other methods given cut-offs based on different criteria.

##  6 Analyses of Proteomic Networks under Immunogenic Heterogeneity

Key scientific questions and dataset overview. Aberrant protein-protein
interactions (PPIs) are associated with various diseases including cancer
(Cheng et al., 2020), and immune cells around the tumor can modulate
malfunctioning PPIs to influence tumor growth and progression (Joyce and
Fearon, 2015). In cancer, cells around the tumor form the tumor
microenvironment (TME) that closely interacts with the tumor (Whiteside,
2008). For example, the dysregulated PPIs in tumor suppress multiple immune
cells in TME to escape the detection from immune system (Whiteside, 2008)
while immune cells in TME can alter the aberrant PPIs to eliminate cancerous
cells (Joyce and Fearon, 2015). This demonstrates the connection between the
dysregulated PPIs and the TME and shows the importance of immunogenic
heterogeneity in tumor behavior. A better understanding of the impact of the
immune cells on aberrant PPIs offers a foundational paradigm for potential
targeted therapies in cancer (Cheng et al., 2020). To this end, our key
scientific questions were as follows: (i) identify important PPIs across
different cancer types and (ii) discover the effect of immunogenic
heterogeneity on aberrant PPIs as potential targets for future investigation.

We exemplify the utility of rBGR, using data from The Cancer Genome Atlas
(TCGA) to build patient-specific PPI networks and investigate the impact of
immunogenic heterogeneity across two different cancers. Specifically, we used
reverse-phase protein array for proteomic data ($\bm{Y}$) to build the PPI
network of CSIx graph and incorporated the immune cell transcriptome
signatures as covariates ($\bm{X}$) as markers of immunogenic heterogeneity.
Our analysis focuses on ovarian cancer (OV) and lung adenocarcinoma (LUAD) as
representative examples of two different types of cancers that elicit distinct
immune responses. OV represents a immunologically “cold” tumor with a weaker
immune response, while LUAD is considered a immunologically “hot” tumor with a
stronger immune response (Galon and Bruni, 2019). We focus on proteins in $12$
important cancer-related pathways (Ha et al., 2018) and obtained $p=41$
proteins with $n=241$ and $n=360$ patients for OV and LUAD, respectively. For
covariates, we included mRNA-derived immune cell gene signatures and
quantified the immune cell abundance corresponding to T cells and two crucial
members of myeloid-derived suppressor cells, monocytes and neutrophils, for
both OV and LUAD. Both T cells and myeloid-derived suppressor cells are
essential in both OV and LUAD since T cells is the main immune component that
kills cancer cells while myeloid-derived suppressor cells regulates T cells
(Whiteside, 2008). We ran rBGR on OV and LUAD with $20,000$ iterations and
discarded the first $19,000$ iterations. The convergence diagnostics and the
details of data preprocessing procedures are provided in Supplementary
Material Section S4.1.

###  6.1 Population-Level Proteomic Networks

We first focus on the covariate dependent population-level networks for OV and
LUAD that are estimated by $\widehat{\alpha}_{j,k,h}$. The corresponding
networks are shown in in Figure 4 (Panels (A) for LUAD and (B) for OV). We
observed that the number of edges is much less in OV compared to LUAD for all
immune components (T cells: (7, 15), monocytes: (5, 82) and neutrophils: (7,
260) for (OV, LUAD)). This is further evidenced in Panel (C) that shows the
distribution of PIPs for OV and LUAD. Interestingly, we observe that the PIPs
for LUAD are higher than those for OV for all immune components (median of
(OV, LUAD) for T cells: (0.123, 0.271), monocytes: (0.131, 0.307), and
neutrophils: (0.127, 0.380)). The higher PIPs in LUAD imply that immune
components have a greater impact on PPIs in LUAD compared to OV. This finding
is consistent with the existing biology, as LUAD belongs to the immune hot
tumors with a stronger immune response (Galon and Bruni, 2019). Furthermore,
we identify HER2, Rb and Bax as the top three hub proteins with the highest
degree in LUAD. In LUAD, HER2 mutation is associated inferior survival (Pillai
et al., 2017), Bcl-2 family protein including Bax is a prognostic biomarker
(Sun et al., 2017), and Rb mutation predicts poor clinical outcomes (Bhateja
et al., 2019). For OV, AR protein is identified as a hub protein with the
highest degree (AR: 13 with the rest protein $\leq 10$). Recent evidence
supports the critical role of AR for the progression of OV (Zhu et al., 2017).

Population graphs also confer specific information about the interaction
between proteins. For example, we observe an edge between Akt and PTEN with
the highest PIP regulated by T cell for LUAD (Panel (A)) suggesting the impact
from T cell on the PPI between Akt and PTEN. It is well-known that PTEN down-
regulates Akt and the loss of tumor suppressor PTEN often leads to
dysregulated PI3K pathway including Akt and the following tumor growth for
LUAD (Conciatori et al., 2020). For OV, despite the smaller number of PPIs, we
still identify important PPIs. For example, rBGR suggests a PPI regulated by T
cells between Caveolin-1 and PR. In OV, Caveolin-1 is regulated by
progesterone, which is mediated by PR, and suggests a consitent result with
the estimated PPI between Caveolin-1 and PR (Syed et al., 2005). Overall, our
analyses capture important hub proteins and characterize the cancer PPIs, and
the results are highly concordant with the existing cancer literature.

Figure 4: Population-Level Proteomic Networks. The posterior inclusion
probability (PIP) (Panel (C)) and the population networks of PPIs with the
cut-off at $c_{0}=0.5$ for LUAD (Panel (A)) and OV (Panel (B)). For each panel
of LUAD and OV, PPI networks of specific immune component are shown from the
left to right for T cells, monocytes, and neutrophils. The degree of each
protein is shown by the node size – with a bigger node representing a higher
degree of connectivity.

###  6.2 Patient-Specific Networks: $\beta_{jk}(\bm{X}_{i})$

We next focus on patient-specific PPI networks to examine the effect of immune
component abundance ($\bm{X}_{i}$) on PPIs. Specifically, we vary one immune
component with the rest components fixed at their mean and generate networks
of CSIx for different individuals at five percentiles (5th, 25th, 50th, 75th
and 95th percentiles) of the varying immune component. We set the cut-off for
the ePP at $c_{1}=0.5$ and show the networks for LUAD in Figure 5 with the
networks for OV in Supplementary Material S4.

We focus on PPIs of CSDx that show dependency on the abundance of specific
immune component. We present PPIs that change the signs in the 5th and 95th
percentiles indicating specific PPIs that are impacted by the immune
components, such as Akt-PTEN for T cells, Bid-PCNA for monocytes, and Bax-
GATA3 for neutrophils. Interestingly, we discovered that the sign of Akt-PTEN
is positively correlated to the T cell abundance. Specifically, when T-cell
abundance is higher, Akt-PTEN is positive; vice-versa, Akt-PTEN is negative
when T cell is scarce. It is well-established that PTEN suppresses Akt
signaling and the loss of PTEN results in the hyper-activation of Akt in
cancer cells and the low T cell abundance in lung cancer (Conciatori et al.,
2020). In addition, we find Bid-PCNA edge is positively correlated with
monocytes abundance. It has been shown that PCNA promotes Bid through caspase
proteins and is crucial to immune evasion in cancers (Wang et al., 2021).
Finally, we discover that Bax-GATA3 edge is positively correlated with
neutrophil abundance. Recently, GATA3 has been found to down-regulate BCL-2
(Cohen et al., 2014), which inhibits the Bax protein (Antonsson et al., 1997),
and neutrophils promotes the Bax to induce the apoptosis (Li et al., 2020).
These findings highlight specific PPIs that are influenced by the abundance of
immune components and suggest potential targets for further investigation of
immunotherapy.

Figure 5: Covariate-specific networks of LUAD under five different percentiles
immune component of (A) T cells, (B) monocytes and (C) neutrophils with the
rest two components fixed at mean zero. The estimated network for varying
immune components are shown from the left to right for 5, 25, 50 ,75, and
95-th percentiles. Edges are identified with signs (green: positive and red:
negative) when the ePPs are bigger than $c_{1}=0.5$.

##  7 Discussion

In this paper, we develop a flexible Bayesian framework called robust Bayesian
graphical regression (rBGR) to construct heterogeneous networks that
explicitly account for covariate-specific information for non-normally
distributed data. By accommodating the non-normal marginal tail behaviors
through random scales, we construct covariate-specific graph through graphical
regression-based approaches. This framework allows us to explicitly
characterize, infer and interpret covariate-specific edge dependencies through
conditional sign independence functions. We also propose an efficient Gibbs
sampler for posterior inference. Our simulations demonstrate that rBRG
outperforms other existing Gaussian-based methods that construct the
covariate-specific graphs under a variety of settings, which display non-
normal marginal behavior such as heavy tails or skewness.

We employ rBGR on proteogenomic datasets in two cancers to build patient-
specific PPI networks and identify PPIs that are impacted by tumor
heterogeneity. Specifically, we quantify the immune cell abundance to
elucidate the effects of immunogeneic heterogeneity on aberrant PPI for lung
and ovarian cancers that are triggered by different levels of immune
responses. Our analyses align with existing biology along three major axes:
(i) immune responses, (ii) hub proteins, and (iii) PPIs. For example, higher
connections in LUAD are consistent with existing biology since LUAD belongs to
the class of the immunlogically “hot” tumors. We identify a hub protein of
HER2, which is associated with a poor survival in LUAD. Another example is a
PPI of Akt-PTEN, which is consistent with the knowledge that PTEN down-
regulates Akt. Our study further suggests PPIs that vary with specific immune
component. For example, we discover that PPIs of Akt-PTEN, Bid-PCNA, and Bax-
GATA3 vary positively with T cells, monocytes and neutrophils, respectively.
These findings suggest potential future targets for immunotherapy in lung
cancer.

In the current implementation, we consider only linear effect of covariates to
reduce the inferential and computation burden. It is possible to include the
non-linear functionals through basis expansion techniques such as splines (Ni
et al., 2019) – however this will increase the computational burden. Another
possible extension is other types of graphical dependencies. For example, a
chain graph considers ordered multi-level structure via directed and
undirected edges (e.g. Chakraborty et al., 2021). By introducing the random
scales and generalizing the regression coefficients as functional
coefficients, the model can include the covariates in the precision matrix to
build the subject-specific chain graphs. Another direction could be include
discrete nodes and the concept of CSIx can be extended for discrete data
(Bhadra et al., 2018). All these directions are left for future
investigations.

Code and Data Availability. We also provide a general purpose code in R that
accompanies this manuscript along with all the necessary documentation and
datasets required to replicate our results (available at
https://github.com/bayesrx/rBGR with necessary documentation in Supplementary
Material Section S5).

## References

  * Airoldi (2007) Airoldi, E. M. (2007).  Getting started in probabilistic graphical models.  PLoS Comput Biol 3, e252.
  * Antonsson et al. (1997) Antonsson, B., Conti, F., Ciavatta, A., Montessuit, S., Lewis, S., Martinou, I., et al. (1997).  Inhibition of Bax channel-forming activity by Bcl-2.  Science 277, 370–372.
  * Baladandayuthapani et al. (2014) Baladandayuthapani, V., Talluri, R., Ji, Y., Coombes, K. R., Lu, Y., Hennessy, B. T., Davies, M. A., and Mallick, B. K. (2014).  Bayesian sparse graphical models for classification with application to protein expression data.  Ann Appl Stat 8, 1443–1468.
  * Bhadra et al. (2018) Bhadra, A., Rao, A., and Baladandayuthapani, V. (2018).  Inferring network structure in non-normal and mixed discrete-continuous genomic data.  Biometrics 74, 185–195.
  * Bhateja et al. (2019) Bhateja, P., Chiu, M., Wildey, G., Lipka, M. B., Fu, P., Yang, M. C. L., et al. (2019).  Retinoblastoma mutation predicts poor outcomes in advanced non small cell lung cancer.  Cancer Med 8, 1459–1466.
  * Chakraborty et al. (2021) Chakraborty, M., Baladandayuthapani, V., Bhadra, A., and Ha, M. J. (2021).  Bayesian robust learning in chain graph models for integrative pharmacogenomics.
  * Cheng et al. (2020) Cheng, S. S., Yang, G. J., Wang, W., Leung, C. H., and Ma, D. L. (2020).  The design and development of covalent protein-protein interaction inhibitors for cancer treatment.  J Hematol Oncol 13, 26\.
  * Chung et al. (2022) Chung, H. C., Gaynanova, I., and Ni, Y. (2022).  Phylogenetically informed Bayesian truncated copula graphical models for microbial association networks.  Ann Appl Stat 16, 2437–2457.
  * Cohen et al. (2014) Cohen, H., Ben-Hamo, R., Gidoni, M., Yitzhaki, I., Kozol, R., Zilberberg, A., and Efroni, S. (2014).  Shift in GATA3 functions, and GATA3 mutations, control progression and clinical presentation in breast cancer.  Breast Cancer Res 16, 464\.
  * Conciatori et al. (2020) Conciatori, F., Bazzichetto, C., Falcone, I., Ciuffreda, L., Ferretti, G., Vari, S., et al. (2020).  PTEN function at the interface between cancer and tumor microenvironment: implications for response to immunotherapy.  Int J Mol Sci 21,.
  * Danaher et al. (2014) Danaher, P., Wang, P., and Witten, D. M. (2014).  The joint graphical lasso for inverse covariance estimation across multiple classes.  J R Stat Soc Series B Stat Methodol 76, 373–397.
  * Dobra and Lenkoski (2011) Dobra, A. and Lenkoski, A. (2011).  Copula Gaussian graphical models and their application to modeling functional disability data.  Ann Appl Stat 5, 969 – 993.
  * Finegold and Drton (2011) Finegold, M. and Drton, M. (2011).  Robust graphical modeling of gene networks using classical and alternative t-distributions.  Ann Appl Stat 5, 1057 – 1080.
  * Finegold and Drton (2014) Finegold, M. and Drton, M. (2014).  Robust Bayesian graphical modeling using Dirichlet $t$-distributions.  Bayesian Analysis 9, 521 – 550.
  * Galon and Bruni (2019) Galon, J. and Bruni, D. (2019).  Approaches to treat immune hot, altered and cold tumours with combination immunotherapies.  Nat Rev Drug Discov 18, 197–218.
  * Ha et al. (2018) Ha, M. J., Banerjee, S., Akbani, R., Liang, H., Mills, G. B., Do, K.-A., and Baladandayuthapani, V. (2018).  Personalized integrated network modeling of the cancer proteome atlas.  Scientific Reports 8, 14924\.
  * Janku (2014) Janku, F. (2014).  Tumor heterogeneity in the clinic: is it a real problem?  Ther Adv Med Oncol 6, 43–51.
  * Joyce and Fearon (2015) Joyce, J. A. and Fearon, D. T. (2015).  T cell exclusion, immune privilege, and the tumor microenvironment.  Science 348, 74–80.
  * Lauritzen (1996) Lauritzen, S. L. (1996).  Graphical Models.  New York : Oxford University Press.
  * Li et al. (2023) Li, M., Li, L., and Kang, J. (2023+).  Bayesian inference of spatially varying correlations via thresholded correlation Gaussian processes.
  * Li et al. (2020) Li, R., Zou, X., Zhu, T., Xu, H., Li, X., and Zhu, L. (2020).  Destruction of neutrophil extracellular traps promotes the apoptosis and inhibits the invasion of gastric cancer cells by regulating the expression of bcl-2, bax and nf-κb.  OncoTargets and Therapy 13, 5271–5281.
  * Liu et al. (2012) Liu, H., Han, F., Yuan, M., Lafferty, J., and Wasserman, L. (2012).  High-dimensional semiparametric Gaussian copula graphical models.  Ann Stat 40, 2293 – 2326.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the Lasso.  Ann Stat 34, 1436 – 1462.
  * Ni et al. (2022) Ni, Y., Baladandayuthapani, V., Vannucci, M., and Stingo, F. C. (2022).  Bayesian graphical models for modern biological applications.  Stat Methods & Appl 31, 197–225.
  * Ni et al. (2019) Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2019).  Bayesian graphical regression.  J Am Stat Assoc 114, 184–197.
  * Ni et al. (2022) Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2022).  Bayesian covariate-dependent Gaussian graphical models with varying structure.  Journal of Machine Learning Research 23, 1–29.
  * Peng et al. (2009) Peng, J., Wang, P., Zhou, N., and Zhu, J. (2009).  Partial correlation estimation by joint sparse regression models.  J Am Stat Assoc 104, 735–746.
  * Peterson et al. (2015) Peterson, C. B., Stingo, F. C., and Vannucci, M. (2015).  Bayesian inference of multiple Gaussian graphical models.  J Am Stat Assoc 110, 159–174.
  * Pillai et al. (2017) Pillai, R. N., Behera, M., Berry, L. D., Rossi, M. R., Kris, M. G., Johnson, B. E., Bunn, P. A., Ramalingam, S. S., and Khuri, F. R. (2017).  HER2 mutations in lung adenocarcinomas: A report from the Lung Cancer Mutation Consortium.  Cancer 123, 4099–4105.
  * Pitt et al. (2006) Pitt, M., Chan, D., and Kohn, R. (2006).  Efficient Bayesian inference for Gaussian copula regression models.  Biometrika 93, 537–554.
  * Scheipl et al. (2012) Scheipl, F., Fahrmeir, L., and Kneib, T. (2012).  Spike-and-slab priors for function selection in structured additive regression models.  J Am Stat Assoc 107, 1518–1532.
  * Sun et al. (2017) Sun, P. L., Sasano, H., and Gao, H. (2017).  Bcl-2 family in non-small cell lung cancer: its prognostic and therapeutic implications.  Pathol Int 67, 121–130.
  * Syed et al. (2005) Syed, V., Mukherjee, K., Lyons-Weiler, J., Lau, K. M., Mashima, T., Tsuruo, T., and Ho, S. M. (2005).  Identification of ATF-3, caveolin-1, DLC-1, and NM23-H2 as putative antitumorigenic, progesterone-regulated genes for ovarian cancer cells by gene profiling.  Oncogene 24, 1774–1787.
  * Wang et al. (2021) Wang, Y. L., Lee, C. C., Shen, Y. C., Lin, P. L., Wu, W. R., Lin, Y. Z., et al. (2021).  Evading immune surveillance via tyrosine phosphorylation of nuclear PCNA.  Cell Rep 36, 109537\.
  * Weinstein et al. (2013) Weinstein, J. N., Collisson, E. A., Mills, G. B., Shaw, K. R., Ozenberger, B. A., Ellrott, K., et al. (2013).  The Cancer Genome Atlas Pan-Cancer analysis project.  Nat Genet 45, 1113–1120.
  * Whiteside (2008) Whiteside, T. L. (2008).  The tumor microenvironment and its role in promoting tumor growth.  Oncogene 27, 5904–5912.
  * Zhang and Li (2022) Zhang, J. and Li, Y. (2022).  High-dimensional Gaussian graphical regression models with covariates.  J Am Stat Assoc 0, 1–13.
  * Zhu et al. (2017) Zhu, H., Zhu, X., Zheng, L., Hu, X., Sun, L., and Zhu, X. (2017).  The role of the androgen receptor in ovarian cancer carcinogenesis and its clinical implications.  Oncotarget 8, 29395–29405.