# Functional Integrative Bayesian Analysis of High-dimensional Multiplatform
Genomic Data

Rupam Bhattacharyyaa, Nicholas C. Hendersona
and Veerabhadran Baladandayuthapania,∗

aDepartment of Biostatistics, University of Michigan, Ann Arbor, MI 48105
∗Corresponding author. E-mail<EMAIL_ADDRESS>
###### Abstract

Rapid advancements in collection and dissemination of multi-platform molecular
and genomics data has resulted in enormous opportunities to aggregate such
data in order to understand, prevent, and treat human diseases. While
significant improvements have been made in multi-omic data integration methods
to discover biological markers and mechanisms underlying both prognosis and
treatment, the precise cellular functions governing these complex mechanisms
still need detailed and data-driven de-novo evaluations. We propose a
framework called Functional Integrative Bayesian Analysis of High-dimensional
Multiplatform Genomic Data (fiBAG), that allows simultaneous identification of
upstream functional evidence of proteogenomic biomarkers and the incorporation
of such knowledge in Bayesian variable selection models to improve signal
detection. fiBAG employs a conflation of Gaussian process models to quantify
(possibly non-linear) functional evidence via Bayes factors, which are then
mapped to a novel calibrated spike-and-slab prior, thus guiding selection and
providing functional relevance to the associations with patient outcomes.
Using simulations, we illustrate how integrative methods with functional
calibration have higher power to detect disease related markers than non-
integrative approaches. We demonstrate the profitability of fiBAG via a pan-
cancer analysis of 14 cancer types to identify and assess the cellular
mechanisms of proteogenomic markers associated with cancer stemness and
patient survival.

Keywords: Bayesian variable selection, cancer genomics, Gaussian processes,
multi-omic data integration, proteogenomic analyses.

##  1 Introduction

Rapid advancements in collection, processing, and dissemination of multi-
platform molecular and genomics (multi-omics, in short) data has resulted in
enormous opportunities to aggregate such data in order to understand, prevent,
and treat diseases. This has catalyzed development of integrative methods that
can collectively mine multiple types and scales of multi-omics data, in order
to provide a more holistic view of human disease evolution and progression
(Subramanian et al., 2020). Specifically, in the context of cancer, a disease
driven predominantly by agglomerations of several molecular changes (Sun et
al., 2021), the importance of synthesizing information from multi-platform
omics and clinical sources to understand the cellular basis of the disease is
even further underscored. Cellular oncological mechanisms, triggered at
different molecular levels of the DNA $\rightarrow$ RNA $\rightarrow$ Protein
path, can confer profound phenotypic advantages/disadvantages. While
significant improvements have been made in multi-omics data integration
methods to unveil such mechanisms, focused on both prognosis (Duan et al.,
2021) and treatment (Finotello et al., 2020), the precise functions governing
them need detailed and data-driven de-novo evaluations. Our work, in the same
vein, aims at two different but inter-related scientific axes: (i) selection
of biomarkers associated with cancer prognosis and clinical outcomes, and (ii)
learning the mechanism of these biomarkers’ effects upon such outcomes via
integrating upstream molecular information - we provide some additional
scientific context below.

#### Classes of Integrative Omics Models

First, we briefly discuss existing integrative omics approaches in order to
contextualize the need for our framework. Broadly, most of the existing
integrative statistical methods can be classified into two categories -
horizontal (meta-analysis type) and vertical (multi-omics) integration
procedures (Tseng et al., 2015). Horizontal meta-analysis methods focus on
integrating data on similar omics features from different sources such as
laboratories, cohorts, sites, etc; examples include works by Tu et al., (2015)
and Angel et al., (2020). Vertical integrative methods, on the other hand, are
focused on integrating data sets on the same cohort of samples obtained from
different omics experiments, wherein the data sets can be vertically aligned;
examples include works by Cheng et al., (2015) and Kaplan and Lock, (2017);
see Richardson et al., (2016) and Morris and Baladandayuthapani, (2017) for a
comprehensive review of integrative methods. Most, if not all such studies
perform the integration in an agnostic manner – they neither take into account
known biological structures nor utilize data illustrating functional roles of
the markers of interest. Incorporating such structures and molecular
regulatory information into integrative models can improve both the power to
detect true biomarkers of a disease and the understanding of their cellular
roles in the progression of it, as we discuss next.

#### Importance of Functional Information

Broadly, by functional information, we mean the knowledge of the specific
molecular functions of the cellular genomic, epigenomic, and transcriptomic
elements, leading to disease outcomes. Incorporation of such information in
biomarker association models is important due to several reasons. First,
different omic components of the molecular configuration of a disease, while
interconnected and hierarchical, can provide complementary information. A
recent review by Buccitelli and Selbach, (2020) indicates how in the DNA
$\rightarrow$ RNA $\rightarrow$ Protein path, termed the ‘gene expression
pathway’, lower or higher correlations of the expressions of proteins and
their coding genes may be observed due to changes in the functional regulatory
elements. Second, specifically for cancer, recent literature indicates that
recurrent regulatory structures drive tumor progression through aberrations at
different omics levels via common master regulatory mechanisms (Califano and
Alvarez, 2017). Finally, experimental validation and characterization of
functional information on a biomarker-by-biomarker basis is resource-intensive
- especially with many plausible candidates. Thus, computational models that
can identify and incorporate functional information about genes/proteins
(referred to as proteogenomic data henceforth) into the models rather than
post-hoc analyses in a natural, inherent way can facilitate this understanding
and can lead to de-novo data-driven prioritization of the relevant biomarkers,
especially for translational and clinical utility. To this end, we propose a
framework called Functional Integrative Bayesian Analysis of High-dimensional
Multiplatform Genomic Data (fiBAG, in short), that allows simultaneous
identification of upstream functional evidence of proteogenomic biomarkers and
the incorporation of such knowledge in Bayesian (biomarker) selection models
to improve signal detection.

#### Goals and Utility of fiBAG

Our scientific goals are multifold. We focus on integrating high-dimensional
multi-omics data and clinical responses in an approach similar to that of Wang
et al., (2013); Jennings et al., (2013), deciphering and delineating
functional roles of proteogenomic markers using mechanism-driven (mechanistic)
models, and incorporating this functional information into outcome models,
thus providing functional relevance to the findings. In particular, we want to
contextualize the available functional information such as the type of
proteogenomic activity following existing literature such as Gevaert et al.,
(2013) and Song et al., (2019). Using functional evidence from these
mechanistic models, we then guide selection and hence the degree of
penalization (thus prioritization) of covariates in the final outcome models.
Figure 1 provides a broad-scale summary of how the fiBAG procedure achieves
these goals. The first column describes the upstream data utilized to infer
functional information. For the purpose of this work, we only use copy number
and DNA methylation; other data platforms having potential functional
relevance (such as microRNA) may also be utilized. The middle column describes
the three axes of mechanistic information that we intend to infer on using
gene and protein expression data, namely, driver gene (dashed line), driver
protein (dashed-dotted line), and cascading protein (dotted line). We describe
the construction of these mechanistic models in more detail in the following
Section(s). The third and final column describes the “calibrated” outcome
model, where summary information from the mechanistic models are incorporated
alongside outcome data to improve selection of genes and proteins. Briefly,
our study offers both methodological novelty via proposing a calibrated
Bayesian variable selection procedure for the outcome model, and scientific
innovation via performing integration of both patient outcomes and tumor
features with multiomic data.

Figure 1: Conceptual schematic summarizing the fiBAG procedure. Panels (A),
(B), and (C) respectively describe the upstream platforms, the mechanistic
models, and the outcome model. The dashed, dotted-dashed, and dotted arrows
indicate the regulatory paths for the driver gene, driver protein and
cascading protein functional axes, respectively.

#### Methodological Novelty

We employ a mapping function to calibrate numerical evidences of significance
obtained from the mechanistic models to a prior inclusion probability scale,
that subsequently inform the outcome model of prior functional evidence in
favor of specific proteogenomic candidates. Our method is flexible – the
calibration function can be adapted according to the choice of the mechanistic
model and the resulting quantification of significance. This hierarchical
evidence sharing procedure allows our method to integrate data across any
number of genomic, epigenomic, and other relevant platforms of choice. Using
Gaussian processes to identify the mechanistic evidence, our model is better
equipped to detect nonlinear cellular associations than a standard linear
model, as used in Wang et al., (2013), and is computationally simpler,
eliminating the needs of choosing the number of knots and incorporating
penalization in a spline-type setting, such as the approach taken by Jennings
et al., (2013). We calibrate the evidence summarized from these models to the
Bayesian variable selection setting by proposing a generalized version of the
spike-and-slab prior originally proposed by George and McCulloch, (1997),
termed the calibrated spike-and-slab prior. This calibrated prior structure
improves the selection of the covariates by borrowing strength across multi-
platform data, choosing to continuously up-weight prior inclusion
probabilities of biomarkers in a data-driven manner. While we take the spike-
and-slab route to build an adaptive and flexible mechanism to incorporate
external knowledge in this work, other existing methods perform the
incorporation of such a priori information in an outcome model via adaptive
shrinkage, penalization, or some different prior structure. Using simulation
studies under multiple synthetic and real data-based scenarios, we compare
both the selection and estimation performances of our method against standard
penalized regression (Tibshirani, 1996), grouped penalized regression
(Boulesteix et al., 2017), and prior-informed selection (Velten and Huber,
2021; Zeng et al., 2021) methods. Our method exhibits comparable selection and
estimation performances with state-of-the-art methods for higher sample size
to number of covariates ratios but accords significantly better false
discovery rate controls, and exhibits substantial improvements in performance
for low-sample high-dimensional settings. We also offer computational
flexibility using both a Markov chain Monte Carlo (MCMC) implementation and a
computationally efficient expectation-maximization based variable selection
procedure (Ročková and George, 2014). Our calibrated Bayesian variable
selection (cBVS) procedure has broad applicability to general regression
problems even outside the multi-omic integration domain where some type of
covariate importance quantification is available, allowing the user to
incorporate prior evidence quantities via other calibration functions tuned to
the scale of such evidence.

#### Scientific Innovation

Multiple works from recent biostatistical literature have focused on
incorporating existing evidence or external information into final models of
interest via various approaches, such as data-adaptive shrinkage (Boonstra et
al., 2015), adaptive Bayesian updates (Boonstra and Barbaro, 2020), or
calibrated maximum-likelihood type procedures (Chatterjee et al., 2016). Our
method is different from such approaches in the sense that it offers a
framework to both learn de novo evidence within the pipeline and incorporate
the said evidence into the final outcome models. As discussed before, omics
elements at different hierarchical levels of the gene expression pathway may
provide partly independent and complementary information (Buccitelli and
Selbach, 2020). Our integrative approach allows learning such information
across interconnected axes of functional activity such as DNA, RNA, and
protein level quantifications. Additionally, our integrative analysis of pan-
cancer proteogenomic data from The Cancer Genome Atlas utilizes both
traditional prognostic outcomes (survival data) and recently developed
cellular descriptors of cancer growth (stemness indices) to identify the
proteogenomic signatures driving such features, along with the molecular basis
of these signatures. Cancer stem-like cells lead to sustained proliferation
via resisting apoptosis, evading growth suppression, and exhibiting increased
invasive and metastatic potential (Fulda, 2013; Adorno-Cruz et al., 2015). We
look at the challenge of identifying the cellular molecular basis of the
behavior of such stem-like cells, by using the mRNA-based stemness index (SI)
proposed by Malta et al., (2018) as an outcome variable. From the survival and
stemness outcome analyses across four common groups of cancers: Pan-
gynecological, Pan-gastrointestinal, Pan-squamous and Pan-kidney, we identify
both known and novel markers associated with the outcomes alongside insights
on their functional roles. In particular, the genes RPS6KA1 (protein p90RSK)
and YAP1 (proteins YAP and YAPPS127) are identified as top driver genes in the
pan-gyne cancers, along with significant associations with the stemness
outcome across multiple cancers in the group; both have been known to be
crucial agents impacting gynecological cancers. Similarly, our analyses found
the gene ERBB2 (protein HER2) to be positively associated with stemness for
gastrointestinal cancers, supported by existing literature.

The rest of the paper is organized as follows. Sections 2.1-2.3 describe fiBAG
both conceptually and with the mathematical details of the mechanistic and
outcome models and the evidence calibration function. Section 2.4 describes
the computational steps behind the fitting, estimation, and selection
procedures for the two models. Section 3 summarizes simulation studies in both
synthetic and real data-based settings comparing our method to existing
benchmarks. Section 4 summarizes results from our pan-cancer integrative
proteogenomic analyses. We conclude the paper with a discussion on the
methodological and biological aspects of our work, along with some potential
future directions in Section 5. All our results are available in an
interactive R-Shiny dashboard (Supplementary Materials Section 5, denoted by
SM henceforth). All the software codes used to perform the analyses, along
with the processed datasets are available as a zip file along with the
submission.

##  2 fiBAG Model

###  2.1 Conceptual Factorization of the Multi-omics Model

We begin with the data structure and some notations. Let $n$ denote the number
of samples, $q_{g}$ be the number of genes, and $q_{p}$ be the number of
proteins in the dataset of interest. The gene (mRNA) expression data matrix
$\mathbf{G}$ and the protein expression matrix $\mathbf{P}$ respectively have
dimensions ${n\times q_{g}}$ and ${n\times q_{p}}$. Further, let the data
matrices corresponding to the upstream covariates copy number alteration and
DNA methylation be denoted as $\mathbf{C}$ and $\mathbf{M}$, respectively,
each having $n$ rows and are matched to their respective genes and proteins
(i.e. cis-level summarizations). Let $\mathbf{Y}$ denote the $n\times 1$
outcome data vector. Thus, the proteogenomic, upstream, and outcome data
available for a cohort of samples can be aligned vertically (i.e. matched
across samples). We write the joint model of the outcome and the proteogenomic
data conditional on the upstream data as

$\mathcal{P}[\mathbf{Y},\mathbf{G},\mathbf{P}|\mathbf{C},\mathbf{M},\boldsymbol{\theta}]=\underbrace{\mathcal{P}[\mathbf{Y}|\mathbf{G},\mathbf{P},\boldsymbol{\theta}_{Y}]}_{Outcome}\underbrace{\mathcal{P}[\mathbf{G},\mathbf{P}|\mathbf{C},\mathbf{M},\boldsymbol{\theta}_{F}]}_{Mechanistic}.$
(1)

Here, $\boldsymbol{\theta}=(\boldsymbol{\theta}_{Y},\boldsymbol{\theta}_{F})$
denotes a conceptual parameter (possibly multi-dimensional) that connects the
two layers of models. The omics-only part
($\mathcal{P}[\mathbf{G},\mathbf{P}|\mathbf{C},\mathbf{M},\boldsymbol{\theta}_{F}]$)
represents the mechanistic model, which concerns the functional mechanisms of
proteogenomic expression as driven by DNA-level cellular activities. Via the
parameter $\boldsymbol{\theta}_{F}$, this mechanistic information is then
learned and incorporated into the outcome model
($\mathcal{P}[\mathbf{Y}|\mathbf{G},\mathbf{P},\boldsymbol{\theta}_{Y}]$). The
parameter vector $\boldsymbol{\theta}$ enables the mechanistic layer to inform
the outcome layer, in line with our scientific aims of integrating functional
information.

#### Biological Rationale for the Factorization

The connection between the two layers drives potentially improved
identification of proteogenomic features in the final outcome model. This
conceptual framework aligns with the idea that different tumors, although
potentially driven by changes in different agents or genomic locations, are
controlled by a recurrent regulatory architecture where genomic alterations
cluster upstream of functional proteins (master regulators) (Califano and
Alvarez, 2017). The interconnections between these regulators form the tumor
checkpoints that can potentially be useful as biomarkers and therapeutic
targets. Further, epigenetic and genetic changes determining a cancer cell
state can be intertwined, and the mutual dependencies between such traits can
contribute to tumor progression via sequential layers of cellular activity
(Alizadeh et al., 2015). Thus, it makes sense to model multi-platform omics
data in a way where functional contributions to the variability in expressions
of genes and proteins can be inferred separately and can be utilized to
calibrate the identification of the roles of those genes and proteins in tumor
progression. Over the next few subsections, we describe the specific
approaches undertaken in this work to formulate these two model layers to
utilize this biological framework.

###  2.2 Mechanistic Models

We focus on three axes of mechanistic information (one for each gene and two
for each protein) based on the available upstream data (Figure 1). We first
describe the mathematical settings and the interpretations of the three axes.
Let $j$ denote the index for the specific proteogenomic biomarker of interest
in a mechanistic model, with the understanding that it is a gene if
$j\in\\{1,\ldots,q_{g}\\}$, and a protein if
$j\in\\{q_{g}+1,\ldots,q_{g}+q_{p}\\}$. For biomarker $j$ and sample $i$, let
the corresponding sub-vectors of $\mathbf{C}$ and $\mathbf{M}$ be denoted by
$\mathbf{C}_{ij}$ and $\mathbf{M}_{ij}$, respectively. We also assume that all
proteogenomic expression data are mean-centered. The general forms of the
models corresponding to the three mechanistic axes are as follows:

$\displaystyle\textbf{Driver gene model:}\;\;\;G_{ij}$ $\displaystyle=$
$\displaystyle
f_{1j}\big{(}(\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T}\big{)}+e_{1ij},$
$\displaystyle\textbf{Driver protein model:}\;\;\;P_{ij}$ $\displaystyle=$
$\displaystyle
f_{2j}\big{(}(\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T}\big{)}+e_{2ij},$
$\displaystyle\textbf{Cascading Protein model:}\;\;\;P_{ij}$ $\displaystyle=$
$\displaystyle
f_{3j}\big{(}(G_{ij},\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T}\big{)}+e_{3ij},$

where $f_{\bullet}$ are nonlinear functionals and $e_{\bullet}$ denote the
(normal) error components. The upstream to gene correspondences are defined by
the physical co-location of the coding segment of the gene in the genome,
where-in a window of $\pm 500$ kb are taken into account. The gene to protein
correspondences are defined using which gene codes for which protein. The
upstream to protein correspondences are defined analogously. The three models
can be biologically interpreted in the following manner:

  1. 1.

Driver gene: Whether the regulations corresponding to the gene are unique at
transcript levels, and whether there is a significant relationship between the
genomic/epigenomic events and the resulting gene expression which, in turn,
drives cancer progression and outcomes (Gevaert et al., 2013).

  2. 2.

Driver protein: Whether the regulations corresponding to the protein are
unique at transcript levels, and whether there is a significant relationship
between the genomic/epigenomic events and the resulting protein expression
which, in turn, drives cancer progression and outcomes.

  3. 3.

Cascading protein: Whether the protein-specific regulations transit through
multiple prior omics levels, i.e., whether there is a cascade of effects via
the DNA $\rightarrow$ RNA $\rightarrow$ protein path (Song et al., 2019).

The middle column of Figure 1 describes the data structures as presented in
the above equations. For each model, we are interested in a null hypothesis of
the type $f_{\bullet}=0$, indicating that the corresponding proteogenomic
expression does not depend on the covariates i.e. $\mathbf{C}$ and
$\mathbf{M}$. Compelling evidence against such hypotheses would provide
evidence for the corresponding mechanistic activity being present for the
gene/protein in question. For example, a strong level of evidence for the
driver gene model would mean that the upstream DNA-level events impact the
expression of the gene of interest significantly; the other models can be
interpreted similarly. We now describe the characterization of these
relationships via suitable modeling choices for the $f_{\bullet}$s, and the
hypothesis testing setting to quantify the strength of evidence in the data.

#### Characterization of $f$ via Gaussian Processes

The flexibility of the mechanistic model setting lies in the freedom in
choosing the specific form of $f$. The simplest choice would be a linear
function - which would lead to multiple linear regression models. These
models, explored by Wang et al., (2013), are easy to handle computationally
and easy to interpret, since the regression parameters are explicitly
available for inference (estimation/testing). However, they could miss
potentially nonlinear associations prevalent across the multi-omics levels of
cellular activity, as has been shown previously (Solvang et al., 2011;
Litovkin et al., 2014). A more sophisticated choice of $f$, allowing such
nonlinear associations would be to use a set of basis functions to describe
$f$ (such as using a spline model, as explored by Jennings et al., (2013)).
Such models, however, require specifications of the knots based on a priori
knowledge and demand additional penalization to obtain a stable fit, rendering
the procedure to be more computationally intensive and less interpretable. To
allow computationally tractable and interpretable identification of nonlinear
associations while avoiding the need to specify knot locations over a
multivariate domain, we use Gaussian process (GP) models. GPs have been
utilized in the context of genomic data in past literature including modeling
gene expression dynamics (Rattray et al., 2019), transcriptional regulations
(Lawrence et al., 2007), and pathway analyses (Liu et al., 2007). Our
simulation studies indicate that in scenarios with a high degree of non-
linearity among the covariates in the generating model, GPs are better
equipped in capturing significant associations than linear models (Section 3).

We now describe the GP specifics for the driver gene mechanistic model here;
the other models can be expressed similarly. To recall, the $j^{\textrm{th}}$
driver gene model is written as
$G_{ij}=f_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T})+e_{1ij}$, where
we assume $e_{1ij}\overset{\textrm{iid}}{\sim}\textrm{N}(0,\tau_{1j}^{2})$.
Let us also denote
$f_{1j}^{(i)}=f_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T})$ for all
$i$. Then, the GP prior on $f_{1j}$ is placed as follows:

GP prior:
$\displaystyle(f_{1j}^{(1)},\ldots,f_{1j}^{(n)})^{T}\sim\mathbf{N}(\mathbf{0},\mathbf{K}_{1j}),$
Covariance matrix:
$\displaystyle\mathbf{K}_{1j(i,k)}=K_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T},(\mathbf{C}_{kj}^{T},\mathbf{M}_{kj}^{T})^{T}),$
Kernel function: $\displaystyle
K_{1j}(\mathbf{u},\mathbf{v})=g\tau_{1j}^{2}\textrm{exp}\big{(}-\frac{||\mathbf{u}-\mathbf{v}||^{2}}{\lambda_{1j}^{2}}\big{)}.$

The hyperpriors specify $\tau_{1j}^{2}\sim\textrm{Inverse-
Gamma}(\frac{\nu_{01j}}{2},\frac{\nu_{01j}\tau_{01j}^{2}}{2})$ and
$\lambda_{1j}\sim\textrm{exp}(\lambda_{01j})$. For all our analyses, we set
$g=n$. Although we use the standard squared exponential kernel, a common
default choice (Micchelli et al., 2006), other kernels can be adopted as well.
As stated before, the gene expression data is mean-centered.

#### Hypothesis Tests for Drivers and Cascades via Bayes Factors

We describe our Bayesian hypothesis testing procedure for the mechanistic
models using our driver gene models; the other models can be tested similarly.
For the $j^{\textrm{th}}$ driver gene model, we test $H_{0j}:f_{1j}=0\textrm{
vs }H_{1j}:H_{0j}\textrm{ is false (equivalent to
}H_{0j}:K_{1j}(\bullet,\bullet)\equiv 0)$. We compute a log of the Bayes
factor (lBF) corresponding to the comparison of the full model vs the null
(zero-mean, error-only) model to perform the test. Bayes factors (BFs) are
particularly useful in this setup from both a statistical and a scientific
point of view. From a methodological perspective, elegant almost sure
convergence results for Bayes factors are available for rather general
settings under standard assumptions (Chatterjee et al., 2020). Further, Bayes
factors have been successfully utilized to quantify significance and compare
model performances in omics models in past literature (Stephens and Balding,
2009). For our case, the final expression of the lBF for the driver gene model
$j$, ignoring the constants $a_{n},b_{n},c_{n}$, and $a$ (dependent on data
dimensions and hyperparameters), is given below (detailed derivations in SM
2.1). The integral in Equation 2 is computed numerically. For any matrix
$\mathbf{A}$, $\mathbf{A}_{\bullet j}$ denotes the $j^{\textrm{th}}$ column of
$\mathbf{A}$, and $\mathbf{A}_{i\bullet}$ denotes the $i^{\textrm{th}}$ row of
$\mathbf{A}$.

$\displaystyle\textrm{lBF}_{1j}=\big{[}a_{n}+b_{n}\ln\Big{(}a+\sum_{i=1}^{n}G_{ij}^{2}-\frac{(\sum_{i=1}^{n}G_{ij})^{2}}{c_{n}}\Big{)}$
(2) $\displaystyle+$
$\displaystyle\ln\int_{0}^{\infty}\exp(-\lambda_{01j}\lambda_{1j})\frac{2^{b_{n}}|\mathbf{K}_{1j}/\tau_{1j}^{2}+I|^{-\frac{1}{2}}}{\big{\\{}\mathbf{G}_{\bullet
j}^{T}(\mathbf{K}_{1j}/\tau_{1j}^{2}+I)\mathbf{G}_{\bullet
j}+a)\big{\\}}^{b_{n}}}d\lambda_{1j}\big{]}/\ln(10).$

We decide the strength of the evidence posed by the lBF from a mechanistic
model using the following standard significance ranges: $<0.5$ (no evidence),
$0.5-1$ (substantial), $1-2$ (strong), and $>2$ (decisive) (Kass and Raftery,
1995). For each gene, we have one lBF from the driver gene model, whereas for
each protein, we have a maximum of two such lBFs from the driver and cascading
protein models. We now describe our approach to calibrate these evidence
quantities into the variable selection models for the outcomes of interest.

###  2.3 Calibrated Bayesian Variable Selection

Functional information is captured by the lBFs from the mechanistic models -
each lBF summarizes the strength of evidence for the functional role of the
corresponding gene/protein. The lBF metric is particularly useful for two
reasons - first, it provides us with an interpretable, unidirectional and
continuous scale of evidence strength, and second, across a large
proteogenomic panel, it provides a highly parallelizable procedure to gather
scalar evidence of functional relevance which are useful numerical quantities
in their own terms and pertinent prior information for future models, such as
our outcome model.

Thus, following Equation 1, we want to incorporate this mechanistic
information into our final outcome model. We pose this problem in context of a
general regression framework where some quantitative summaries of covariate
importance are available beforehand, and such information is to be combined
with the mechanics of a typical variable selection setting. Generally, we
denote such prior information as $\boldsymbol{\mathcal{E}}$, with
$\boldsymbol{\mathcal{E}_{j}}$ denoting the possibly multi-dimensional
evidence summary for covariate $j$ (i.e., lBF for gene/protein $j$ in our
case). We intend to inform the final outcome model using these evidences in
terms of selection/non-selection of the predictors. Specifically, if there is
sufficiently strong evidence in favor of a covariate, we want to up-weight its
probability of inclusion. Otherwise, we want to put a uniform probability on
selection/non-selection for that particular covariate. To achieve this, we
utilize a hierarchical Bayesian framework with spike-and-slab priors for each
effect, with the spike probabilities calibrated using the evidence available.
The rest of this subsection describes the components of our fiBAG outcome
model, called calibrated Bayesian variable selection, or cBVS in short.

#### Notations

Following notations introduced at the beginning of this section, let $Y_{i}$
denote the outcome for individual $i$. For the purpose of exposition, we
assume that the outcome of interest is continuous, and each $Y_{i}$ is assumed
to be observed. Generalizations to censored or categorical outcomes are
straightforward. Let us denote the design matrix corresponding to any
additional covariates (other than genes and proteins, such as clinical
information or demographics) by $\mathbf{B}_{n\times q_{b}}$. The combined
design matrix will have a dimension $p=1+q_{b}+q_{g}+q_{p}$. We then propose a
hierarchical Bayesian outcome model, as described below.

$\displaystyle Y_{i}$ $\displaystyle=$
$\displaystyle\beta_{0}+\mathbf{B}_{i\bullet}^{T}\mbox{\boldmath$\beta$}_{B}+\mathbf{G}_{i\bullet}^{T}\mbox{\boldmath$\beta$}_{G}+\mathbf{P}_{i\bullet}^{T}\mbox{\boldmath$\beta$}_{P}+\epsilon_{i},\quad\forall
i\in\\{1,...,n\\},$ $\displaystyle\epsilon_{i}$
$\displaystyle\overset{\textrm{iid}}{\sim}$
$\displaystyle\textrm{Normal}(0,\sigma^{2}),\quad\forall i\in\\{1,...,n\\}.$

Let $\mbox{\boldmath$\beta$}_{p\times 1}$ denote the complete regression
coefficient vector. We set up a hierarchical calibrated spike-and-slab prior
on its last $q_{g}+q_{p}$ components, as described below. A standard conjugate
prior is put on the residual variance parameter as
$\sigma^{2}\sim\textrm{Inverse-Gamma}(\tfrac{\nu}{2},\tfrac{\nu\lambda}{2})$.

$\displaystyle(\beta_{0},\mbox{\boldmath$\beta$}_{B}^{T},\mbox{\boldmath$\beta$}_{G}^{T},\mbox{\boldmath$\beta$}_{P}^{T})^{T}=\mbox{\boldmath$\beta$}|\mbox{\boldmath$\gamma$}_{q_{g}+q_{p}\times
1},\sigma$ $\displaystyle\sim$
$\displaystyle\textrm{Normal}(\mathbf{0},\mathbf{D}_{\mbox{\boldmath$\gamma$},\sigma}),$
$\displaystyle\gamma_{j}|\omega_{j}$ $\displaystyle\sim$
$\displaystyle\textrm{Bernoulli}(\omega_{j}),\quad\forall
j\in\\{1,...,q_{g}+q_{p}\\},$ $\displaystyle\omega_{j}$ $\displaystyle\sim$
$\displaystyle\textrm{Beta}\Big{(}\mathcal{F}(\boldsymbol{\mathcal{E}}_{j}),\frac{1}{\mathcal{F}(\boldsymbol{\mathcal{E}}_{j})}\Big{)},\quad\forall
j\in\\{1,...,q_{g}+q_{p}\\},$

where
$\mathbf{D}_{\mbox{\boldmath$\gamma$},\sigma}=\sigma^{2}\mathbf{A}_{\mbox{\boldmath$\gamma$}}$,
$\mathbf{A}_{\mbox{\boldmath$\gamma$}\;p\times
p}=\textrm{diag}\\{v_{1}\mathbf{1}_{1+q_{b}},\gamma_{1}v_{1}+(1-\gamma_{1})v_{0},\ldots,\gamma_{q_{g}+q_{p}}v_{1}+(1-\gamma_{q_{g}+q_{p}})v_{0}\\}$
where $v_{1}\geq v_{0}>0$ are respectively the slab and spike variances, and
where $\mathcal{F}$ is a calibration function mapping the evidence
$\boldsymbol{\mathcal{E}}_{j}$ to the prior inclusion probability
$\omega_{j}$.

#### Calibration Functions incorporating External Information

The function $\mathcal{F}$ in the outcome model maps the functional evidence
from the mechanistic models to a parameter scale to inform the Beta hyperprior
for selection of each covariate. The specific mathematical properties and form
of this function depends on the user-specific context of what type, range, and
scale of functional evidence is being used. Further, it is possible to garner
multiple quantities of evidence from different sources for each covariate of
interest, and depending on the scenario, the calibration function may take
inputs in $\mathbb{R}^{K}$ for some $K>1$.

For the suitability of exposition, we describe a specific form of the mapping
function used for our analyses where it is assumed that the evidence measures
$\mathcal{E}_{j}$ are in the form of lBF’s from the mechanistic models for
genes and proteins. Our calibration function then needs to be a map from
$\mathbb{R}\rightarrow\mathbb{R}$ that can reflect the changes in the strength
of evidence in the lBFs in the desired ranges. For lBFs, the strength of
evidence is unidirectional (i.e., higher evidence quantity implies higher
evidence strength). In such a scenario, a mapping function is required to be
non-decreasing. Such a mapping function can either be discrete or continuous,
and the jumps of the discrete curve or the slope of the continuous curve can
be tuned depending on the importance to be put on different scalar values of
the evidence. In particular, we map the standard lBF significance ranges as
described before to a parameter describing a Beta prior distribution for the
covariate-inclusion probabilities $\omega_{j}$. The shape of the lBF to Beta
prior mean curve and the representative densities for some lBF values of
interest in each of the standard lBF significance ranges ($<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), and $>2$ (decisive)) for
this calibration function are summarized in Figure 2. Effectively, we build a
function having the following properties.

  * •

If the mechanistic models do not provide us with enough evidence regarding the
functional utility of a gene/protein, we put a uniform prior (mean prior
probability $=0.5$) for the corresponding selection probability parameter in
the outcome model. This is illustrated by the leftmost point on the X axis of
Figure 2A.

  * •

If the mechanistic models provide us with strong evidence, on the other hand,
we put a strong prior with a large prior mean on the selection probability
parameter. This is illustrated by the increment in the Y axis of Figure 2A.

Specifically, we use a four-parameter logistic function (as shown in Figure
2A) for our applications and the details on choosing the calibration function
are summarized in SM 2.2. For this specification, the prior mean probability
of inclusion for a covariate ranges from $0.5$ (when the lBF is $0$ or
negative) to approximately $1$ (when the lBF is large, say $>5$). The
calibrated prior mean probabilities for some representative points from these
intervals are as follows, with respect to Figure 2B: 0.502 (lBF $=0.25$,
bottom-right), 0.543 (lBF $=0.75$, top-right), 0.726 (lBF $=1.5$, bottom-
left), 0.962 (lBF $=3$, top-left). This construction enables the functional
evidence to inform selection/estimation in the outcome model, while retaining
the flexibility to still allow the model to ignore prior evidence, if
necessary, for the proteogenomic markers without any observed association with
the outcome of interest, as are illustrated by our pan-cancer applications
presented in Section 4.

Figure 2: Calibration function for the outcome model. Panel (A) plots the
calibrated prior Beta means against the lBFs. Panel (B) presents the densities
(means indicated by solid vertical lines, broken vertical lines indicate the
$\pm 1$ sd ranges around the means) for four representative values from the
four ranges of interest in the x-axis of panel (A).

#### Benefits and Utilities of Calibration

The benefits of this calibrated model formulation are multifold. First, by
modifying the calibration function $\mathcal{F}$, our framework allows the
user to incorporate any form of model summary information (such as z-scores,
p-values, etc.) into the outcome model. Unlike existing grouped shrinkage-
based procedures, this eliminates the restriction of only using external
information through categorical covariates. Second, unlike the shrinkage-based
or Bayesian approaches where the external covariates directly inform the final
model, our model only relies on the supply of the summary statistics $s_{j}$
\- rendering the computations less challenging and allowing the use of any
upstream dataset, however large, in learning the functional information from
multiple categorical/continuous sources. In the same vein, any such
calibration function can easily be adapted to a scenario where multiple
sources of evidence are available for a single covariate rather than just a
single summary statistic $s_{j}$. For example, if there are $K$ lBFs
($\boldsymbol{\mathcal{E}}_{j}=(\textrm{lBF}_{1j},\ldots,\textrm{lBF}_{Kj})$)
available from as many mechanistic models for the $j$th biomarker, then we can
think of the calibration function as a composition:
$\mathcal{F}(\boldsymbol{\mathcal{E}}_{j})=\mathcal{F}_{1}(\mathcal{F}_{2}(\boldsymbol{\mathcal{E}}_{j}))$,
where $\mathcal{F}_{1}$ can be similar to the $\mathcal{F}$ we use, and
$\mathcal{F}_{2}:\mathbb{R}^{E}\rightarrow\mathbb{R}$ aggregates the multiple
lines of evidence to a single scalar value $s_{j}$. One possible choice for
this is a linear map, as
$s_{j}=\mathcal{F}_{2}(\boldsymbol{\mathcal{E}}_{j})=\sum_{k=1}^{K}\alpha_{kj}\textrm{lBF}_{kj}$.
Here $\alpha_{kj}$s are convex weights specific to biomarker $j$, interpreted
as quantifications of the importance of each source of evidence. Several
choices of the $\alpha_{kj}$s are possible, as described below.

  1. 1.

Average evidence: $\alpha_{kj}=1/K$ (takes a simple average of all available
evidences).

  2. 2.

Maximal evidence:
$\alpha_{kj}=I\big{(}\textrm{lBF}_{kj}=\underset{k^{\prime}\in\\{1,\ldots,K\\}}{\max}\textrm{lBF}_{k^{\prime}j}\big{)}$
(only takes into account the strongest evidence available from any source).

  3. 3.

Precision-weighted evidence: $\alpha_{kj}=\rho_{kj}/\sum_{k=1}^{K}\rho_{kj}$
(weights the evidences by some metric of reliability of the evidences, such as
$\rho_{kj}=\hat{\tau}^{-2}_{kj}$ where $\hat{\tau}^{2}_{kj}$ is the estimated
noise variance for the source mechanistic model of $\textrm{lBF}_{kj}$).

This provides an additional layer of flexibility that allows the user to
choose how many sources of evidence to use and how best to combine them.

#### Parameters of Interest and Model Dependence Structure

Having described the mechanistic and outcome layers of the modeling framework
and established the connection between them via the evidence calibration
function, we now revisit Equation 1 to interpret the conceptual parameters in
context of our modeling scheme. For the mechanistic models, the quantities
$\boldsymbol{\theta}_{F}$ represent the parameter vector for the Gaussian
process models, namely, $\\{(\tau_{lj},\lambda_{lj}):l$ varies across each
mechanistic model available for each gene/protein $j\\}$. The quantification
of evidence is not performed via direct estimation of these parameters but via
computing and calibrating the Bayes factors corresponding to each model.
Therefore, one component of $\boldsymbol{\theta}_{Y}$ is driven by
$\boldsymbol{\theta}_{F}$ via the evidence calibration. The rest of
$\boldsymbol{\theta}_{Y}$ is specified by
$(\mbox{\boldmath$\beta$}^{T},\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$,
i.e. the parameters of interest from the outcome model. The inferential task
then is to fit the outcome model, estimate the parameters of interest
$(\mbox{\boldmath$\beta$}^{T},\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$,
and perform covariate selection based on the posterior distribution of
$(\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T})^{T}$. The overall
dependence structures in the mechanistic and outcome model settings are
summarized in Supplementary Figure 1 (denoted by SF henceforth). In the next
section, we describe our model-fitting procedures.

#### Generalizations to Non-Gaussian Outcomes

Such generalizations typically involve minimal changes to the outcome model.
As an example, we discuss an extension to survival outcomes briefly. Let us
denote $Y_{i}=\log T_{i}$ as the possibly unobserved log-survival time and
$C_{i}$ as the possibly unobserved censoring time for individual $i$. The
observed response for individual $i$ is $(Z_{i},\delta_{i})$, where
$Z_{i}=\min\\{\log T_{i},\log C_{i}\\}$ and $\delta_{i}=I(T_{i}<C_{i})$. Under
the assumptions that the censoring distribution is free of $\beta$ and
$\gamma$ and that the censoring times are independent of the true outcomes
conditional on $\beta$ and $\gamma$, changing the outcome model from Gaussian
to an accelerated failure time model with log-Normal outcomes only introduces
truncations of the censored (unobserved) outcomes in the log-posterior.
Therefore, the only change in the model fitting procedure occurs in the
expression of the full log-posterior, which now includes the terms
corresponding to the truncations.

###  2.4 Model Fitting and Parameter Estimation

#### Mechanistic Model Fitting

As described in Section 2.2, we are interested in quantifying the functional
evidence for each gene/protein via the GP-based mechanistic model using the
lBF - therefore, we do not require a full Bayesian exploration of the model.
To ensure computational efficiency, we directly compute the lBF following
expressions as in Equation 2. This requires the evaluation of an integral,
which we perform numerically.

#### Outcome Model Fitting

As described in the previous subsection, the complete set of parameters to be
estimated by cBVS is
$(\mbox{\boldmath$\beta$}^{T},\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$.
$\mbox{\boldmath$\beta$}^{T}$ provides estimates of the effect sizes of the
proteogenomic covariates on the outcome, and the
$(\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T})^{T}$ guides their
selection. The parameter estimation is focused on the posterior of
$(\mbox{\boldmath$\beta$}^{T},\mbox{\boldmath$\gamma$}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$.
For large proteogenomic panels, this posterior will be computationally
resource-intensive to directly sample from. The question, therefore, is of a
trade-off between computational simplicity and estimation accuracy. Due to
this reason, we offer three implementations of cBVS in increasing order of
computational efficiency – a standard Markov chain Monte Carlo (MCMC) using
Gibbs sampler to sample from the complete posterior, a selection-only MCMC to
sample from the marginal posterior of $\gamma$, and an expectation-
maximization based variable selection (EMVS) procedure to approximate the
posterior modes of the parameters. Briefly, the selection-only MCMC focuses on
estimating $\gamma$ first and then estimates $\beta$ using a Bayesian model
averaging-type procedure (Hinne et al., 2020), and the EMVS sacrifices the
full posterior along with error estimates to achieve fast point estimation
(Ročková and George, 2014). The exact details of each are in SM 2.3.

#### Model Summaries

The MCMC/EMVS procedures provide us with posterior inclusion probabilities
(PIPs) $\hat{\omega_{j}}$ and regression coefficient estimates
$\hat{\beta_{j}}$ for each covariate in the model. A cut-off on the PIPs is
computed using a false discovery rate adjustment procedure at a specified
level of significance, treating the $1-\hat{\omega_{j}}$ as p-value type
quantities (SM 2.4).

##  3 Simulation Studies

To illustrate the utilities of cBVS, we performed two sets of simulation
studies. Simulation 1 deals with continuous outcomes in synthetic datasets,
comparing metrics of selection and estimation from the cBVS procedure with
existing benchmarks across a class of variable selection methods. We include
both penalized/grouped penalized selection procedures and Bayesian prior-based
selection procedures as benchmarks. Simulation 1 is expected to assess cBVS
for both low and high sample size to dimension ratios and quantify the
improvement and/or preservation of performance along this spectrum compared to
the benchmarks. In Simulation 2, the data generation procedure is informed by
the patient datasets for breast invasive carcinoma from our application. The
next two subsections describe the design, settings, and results from these two
simulation studies.

###  3.1 Simulation 1: Synthetic data-based Simulations

#### Data Generation and Choices of True Effects

To compare performances across a grid of varying sample size/number of
covariates ($n/p$) ratios, we fix the number of proteogenomic biomarkers
generated in the simulated datasets at $p=200$ and vary the sample size across
$n=50,100,200,400,800$, covering a range of $n/p=1/4$ to $4$. For simplicity,
we assume that there is only one upstream covariate for each biomarker. 100
replicates are generated for each n. Briefly, we generate the upstream
covariates (akin to copy number/methylation data) from the standard normal
distribution. The proteogenomic covariates are generated such that there are
12 groups of size 5 each – all possible combinations of effect size (low,
medium, high) and prior functional evidence (no evidence, substantial, strong,
decisive). The rest all have zero true effect sizes. The continuous outcome is
generated from a normal distribution. Further details for each step are
available in SM 3.1.

#### Brief Overview of Benchmark Methods

We compare the performance of the cBVS model against three classes of methods
that perform variable selection based on different approaches. We use one
Bayesian variable selection method without any external evidence, namely, the
expectation-maximization based variable selection (EMVS) (Ročková and George,
2014) - as an uncalibrated counterpart to cBVS. We also include two penalized
selection procedures - additional to LASSO (Tibshirani, 1996), we include a
grouped penalized selection procedure termed Integrative LASSO with Penalty
Factors (IPF-LASSO) (Boulesteix et al., 2017). The reason behind including a
grouped penalized regression procedure is that the prior evidence levels
provide a natural grouping for the covariates to be used in a variable
selection setting. Finally, we implement two recently proposed Bayesian
variable selection methods incorporating external information, namely, graper
(Velten and Huber, 2021), which incorporates categorical external covariates
only, and xtune (Zeng et al., 2021), which can handle continuous covariates as
well. Our simulation scheme (i) allows a comprehensive comparison of
calibrated vs non-calibrated methods and (ii) provided a benchmark for
calibration to evaluate cBVS.

#### Summary of Metrics Used

Each method provides a coefficient estimate $\hat{\mbox{\boldmath$\beta$}}$.
For LASSO, we get a $\hat{\mbox{\boldmath$\gamma$}}$ corresponding to the
best-fit $\lambda$. For each of the other methods, we compute
$\hat{\mbox{\boldmath$\omega$}}$, on which we use the FDR-based adjustment
method as described in SM 2.4 to infer selection. We compute several standard
metrics of selection performance, namely, area under the receiver operating
characteristic curve (AUC), scaled AUC between 0.8 to 1 specificity (AUC20),
true positive rate (TPR), false positive rate (FPR), and Matthew’s correlation
coefficient (MCC). The use of MCC is particularly useful since at a specified
selection threshold it aggregates performance summary from both the TPR and
FPR metrics. Using AUC and AUC20 allows threshold-free evaluation of selection
performance.

#### Results and Discussion

The results from Simulation 1 are summarized in Figure 3A. In terms of AUC,
our calibrated outcome regression model performs the best for the smallest
$n/p$ ratio $=1/4$, and is only second to graper by a very small margin for
the $n/p$ ratios $=2,4$. In terms of MCC, our calibrated outcome regression
model is only second to xtune for $n/p$ ratios $=1/4,1/2$, and is only behind
graper and xtune for the $n/p$ ratios $=1,2$, and $4$. cBVS shows promise
particularly in terms of preventing false positives, by exhibiting the lowest
FPR among the method compared for all $n/p$ ratios $\geq 1$. In particular,
cBVS has lower FPR than IPF-LASSO for all simulation scenarios other than
$n/p=1/4$, and comparable or lower FPR than the calibrated procedure xtune for
the same scenarios as well. These summaries indicate several utilities of the
calibrated outcome regression procedure. First of all, based on the metrics
summarized, calibration based on prior evidence seems to have an evident
benefit, as all the methods utilizing prior evidence (cBVS, graper, xtune)
outperform those not incorporating any prior evidence (EMVS, LASSO, IPF-LASSO)
across all the $n/p$ ratios used and across all metrics. Further, as evident
from the AUC summaries, cBVS offers an improved preservation of selection
performance for $n/p$ ratios $<1$ over the methods incorporating external
covariates, while maintaining a comparable and satisfactory level of
performance in the more favorable scenarios, i.e., $n/p>1$, as well. This is
particularly reassuring, since unlike MCC or other standard metrics of
selection, AUC is a threshold-free summary which enables us to sum up the
performance of the methods across the whole spectrum of lenient to stringent
selection criteria.

Figure 3: Summary of results from (A) synthetic and (B) real data-based
simulations with continuous response. For (A), 100 replicates were generated
for each sample size. The number of covariates was 200, with 60 of them having
non-zero true effect sizes, mixed across a grid of low, medium, and high
values, each with equal proportions of lBFs at the levels no evidence,
substantial, strong, and decisive. For (B), 100 replicates were generated as
well. Covariate matrix, effect sizes, and calibrated evidence were picked from
a calibrated model for BRCA (breast invasive carcinoma).

###  3.2 Simulation 2: Real Data-based Simulations

#### Data Generation and Choices of True Effects

Before the pan-cancer multi-platform proteogenomic analysis using data from
the Cancer Genome Atlas, we perform a modified version of Simulation 1, based
on breast invasive carcinoma (BRCA) patient data. We aggregate and annotate
DNA methylation, copy number alteration, proteogenomic expression, and
censored survival data for BRCA following Section 4. Since the rest of the
simulation scheme is analogous to Simulation 1, we only summarize the changes
briefly. The upstream covariate data, proteogenomic covariate data,
gene/protein evidence quantities (lBFs) are all taken from the BRCA patient
data ($n=790$ and $p=365$). These lBFs are then grouped as before ($<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive)). The true
effect sizes are taken to be estimates from a survival cBVS model using the
real data. These effect sizes are again clustered into four groups based on
their quartiles, meaning that we now have $4\times 4=16$ possible
combinations. The outcome data generation procedure, benchmark methods,
summary metrics, and number of simulations remain the same as before.

#### Results and Discussion

The resulting metrics are summarized in Figure 3B. Notably, cBVS performs the
best among all the methods compared in terms of AUC, AUC20 and FPR, and is at
the second place, only behind graper, for MCC. Again, two calibrated selection
methods (cBVS, graper) outperform uncalibrated methods (EMVS, LASSO, IPF-
LASSO) in terms of all metrics but TPR. cBVS provides the lowest FPR among
all, even lower than the FPRs for both the calibrated competitors - graper and
xtune. Thus, cBVS is well-equipped for the real data scenarios to be
encountered in the pan-cancer setting.

#### Additional Simulations

We perform a simulation study to compare the performances of GPs and linear
models in capturing non-linear associations. The generating models, simulation
procedure and results are described in SM 3.2. Our results show that GPs are
better equipped in identifying nonlinear evidence than linear models. We also
perform comparisons between model outputs from the MCMC and EMVS
implementations to ensure EMVS provides estimates comparable to those from the
MCMC. These details are available in SM 3.3. Our results reassure that
although EMVS sacrifices quantifying uncertainty in lieu of faster
computations, it does not compromise in terms of estimation accuracy.

##  4 Integrative Pan-cancer Proteogenomic Analyses

###  4.1 Data Description, Preprocessing and Analysis

#### Pan-cancer Multi-omics Data

We analyze pan-cancer multi-platform proteogenomic data from the Cancer Genome
Atlas (TCGA). We include 14 cancers across four cancer groups, classified by
commonalities in the sites/tissues of occurrence. For each cancer, DNA
methylation, copy number alteration, transcriptomic and proteomic expression
data were obtained. The details on the cancer groups and types are available
in SM 4.1. The codes used for data collation and cleaning along with the
processed datasets are available in the code archive included in the
submission. The sample size summaries for each cancer are available in the
interactive shiny dashboard. The largest and the smallest sample sizes were
available for breast ovarian carcinoma ($n=1188$) and uterine carcinosarcoma
($n=57$). We only used the genes corresponding to the proteins available in
the TCGA reverse-phase protein array (RPPA) panel, resulting in a number of
covariates generally between $350-400$ across cancers, depending on data
availability and quality.

#### Outcome Data

We investigate two different outcomes. We use the censored survival data
available from TCGA and implement the cBVS model to identify cancer-specific
proteogenomic biomarkers associated with survival. However, survival outcome
alone does not provide biologically interpretable insights into the molecular
progression of the disease. To this end, we use another outcome called mRNA-
based stemness index (SI, in short) - a metric of cancer growth in terms of
its cellular features. Briefly, SIs for TCGA cancers are computed using a one-
class logistic regression model trained on pluripotent stem cell samples from
the Progenitor Cell Biology Consortium dataset (Daily et al., 2017; Malta et
al., 2018). The SIs quantify the stem-cell-like behavior of the tumor of
interest. We build cBVS models selecting proteogenomic biomarkers associated
with SIs.

#### Scientific Questions

Our analyses are driven by three broad scientific aims. First, we intend to
identify cancer-specific and pan-cancer functional drivers and cascades from
the proteogenomic candidates using the functional evidence learned via the
mechanistic models. Second, using the cBVS models, we select cancer-specific
biomarkers associated with changes in the survival and stemness outcomes.
Finally, we assess the mechanistic and outcome model results in conjunction to
evaluate the utility of calibrating outcome models using mechanistic evidence.
We present the numerical results in the following two subsections, followed by
the biological interpretations and discussions of the results in Section 4.4.

#### Modeling

For each cancer, the mechanistic model analyses for each gene and protein are
performed using GPs as described in Section 2.2, and lBFs are computed. As
explained in the previous section, we have a single evidence quantity (the lBF
from the driver gene model) for each gene and a maximum of two evidence
quantities (lBFs from the driver and cascading protein models) for each
protein. We use the maximal approach for evidence aggregation as described in
Section 2.3. Two cBVS models - one using stemness and one using survival data
- are built for each cancer, using the selection-only MCMC procedure as
described in Section 2.4 to estimate
$(\mbox{\boldmath$\beta$}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$ in each.
For each gene and protein, we thus obtain five quantities of interest: lBF,
and $(\hat{\beta}_{j},\hat{\omega}_{j})^{T}$ \- one each from each outcome
model. An overview of the analysis pipeline is presented in SF 4.

###  4.2 Mechanistic Model Results for Pan-Cancer Groups

We summarize the mechanistic model outputs across cancers for the pan-gyne and
pan-GI cancer groups in Figure 4; the rest of the groups are presented in SF
5-10. Figure 4A, C exhibit heatmaps summarizing lBF classes for the gene-
protein pairs which have some evidence across at least three-fourths of the
cancers in at least two out of the three mechanistic model types, along with
corresponding upset plots in Figure 4B, D describing the number of
genes/proteins that are at the decisive level of significance for the three
mechanistic models across the intersections of the different cancers.

For the pan-gyne cancers, 26 gene-protein pairs are at strong/decisive level
of evidence across at least four cancers in at least two out of the three
mechanistic model types, including genes such as RPS6KA1 (protein p90RSK),
YAP1 (proteins YAP and YAPPS127), and DIABLO (protein SMAC) (Figure 4A). The
largest sharing of decisive driver gene signatures is observed across BRCA,
CESC, OV, and UCEC, and that for cascading proteins is observed across BRCA,
OV, and UCEC (Figure 4B).

For the pan-GI cancers, eight gene-protein pairs are at strong/decisive level
of evidence across all three cancers in at least two out of the three
mechanistic model types, including genes such as ERBB2 (proteins HER2 and
HER2PY1248), CCNE1 (protein CYCLINE1), and MAPK9 (protein JNK2) (Figure 4C).
The largest number of decisive driver gene, driver protein, and cascading
protein signatures is observed from CORE, and the largest numbers of pan-
cancer signatures for the driver gene and the cascading protein models are
observed between CORE and ESCA (Figure 4D).

Figure 4: Summary of mechanistic model results for the pan-gyne group cancers:
BRCA (breast invasive carcinoma), CESC (cervical squamous cell carcinoma and
endocervical adenocarcinoma), OV (ovarian serous cystadenocarcinoma), UCEC
(uterine corpus endometrial carcinoma), UCS (uterine carcinosarcoma) (panels A
and B), and the pan-GI group cancers: CORE (colon and rectum adenocarcinoma),
ESCA (esophageal carcinoma), STAD (stomach adenocarcinoma) (panels C and D).
(A/C) The lBF ranges are defined as: $<0.5$ (no evidence), $0.5-1$
(substantial), $1-2$ (strong), $>2$ (decisive). Only the gene-protein pairs
with some evidence across at least 75% of the cancers in at least two out of
the three mechanistic model types are shown here. (B/D) Upset plots exhibiting
the number of genes/proteins at the decisive level of significance for the (i)
driver gene, (ii) driver protein, and (iii) cascading protein mechanistic
models.

###  4.3 Stemness cBVS Results for TCGA Cancers

We summarize the stemness outcome model outputs for each TCGA cancer using
plots showing $-\log_{10}(1-\hat{\omega})$ on the y-axis and $\hat{\beta}$ on
the x-axis for each cancer in Figure 5 and SF 11-15. Further, we present bar
diagrams of the lBFs for the selected genes/proteins and histograms with
density plots of lBFs for the predictors not selected in our shiny app. Here,
we discuss the results from the BRCA and CORE analyses since they are the
cancers with the largest sample sizes in the pan-gyne and pan-GI groups,
respectively.

For BRCA, several genes are selected, such as YAP1, DIABLO, and RAPTOR (Figure
5A); however, no proteins are selected using the 10% FDR cut-off on
$\hat{\mbox{\boldmath$\omega$}}$ from the stemness cBVS model. The genes
selected cover a large span of evidence range from the mechanistic models,
with genes like JUN and COL6A1 having no functional evidence and DVL3, MYH11,
EIF4G1 at decisive evidence (SF 16), indicating that cBVS is capable of
identifying associations even in the absence of prior evidence, as has been
noted in the simulation studies before too (Section 3). The histogram
indicates that a vast majority of the non-selected genes and proteins have
little to no evidence from the mechanistic models, but 194 of them are at the
strong and decisive levels of evidence and yet not selected by cBVS due to the
absence of sufficient association with the outcome data (SF 17).

For CORE, two genes - PEA15 and KDR are selected in the stemness cBVS model,
both negatively associated with the stemness index as can be seen from the
sign of the estimated regression coefficients (Figure 5I). Both the genes
selected are at the decisive level of functional evidence (SF 18). The
histogram again indicates that a vast majority of the non-selected genes and
proteins have little to no evidence from the mechanistic models - however, 179
of them are at the strong/decisive levels of evidence and are still not
selected by cBVS (SF 19). The survival cBVS results are discussed in SM 4.2.

Figure 5: Plots summarizing outcome model results based on stemness indices
(SI) for TCGA cancers. Proteins are represented by triangles and genes by
circles. The shapes are colored red if the estimated $\hat{\beta}_{j}$ from
fiBAG is negative, and green if positive. The x-axis shows the
$\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).

###  4.4 Biological Findings and Implications

While a majority of the proteogenomic biomarkers identified in our analyses
have supporting evidence from past literature in terms of their mechanisms, we
further delineate their functional roles in driving cancer progression and
patient survival. The detailed summary of these biomarkers, along with
relevant references from the literature is presented in SM 4.3. All
associations discussed below are significant at a 10% level of FDR control.

#### RPS6KA1 Gene and p90RSK Kinases in Gynecological Cancer Progression

The gene RPS6KA1 (protein p90RSK) has decisive evidence for all but two
mechanistic models in the pan-gyne cancers (Figure 4A). The p90RSK protein is
selected in the survival cBVS model for OV (SF 22). Both RPS6KA1 and p90RSK
have previously shown significant association with risks of and outcomes in
gynecological cancers (SM 4.3).

#### YAP1 Gene and YAP Proteins in Gynecological Malignancies

The gene YAP1 (proteins YAP, YAPPS127) has decisive evidence in all
mechanistic models for all gynecological cancers except UCS (smallest sample
size in the group, $n=57$) (Figure 4A). Further, YAP1 is negatively associated
with stemness for BRCA (Figure 5A) and with survival for OV and CESC (SF
21-22). Over-activation of YAP is a known cause of malignant transformation in
gynecological cancers, associated with poor prognosis (SM 4.3).

#### DIABLO Gene as a Marker of Gynecological Tumors

We identify higher DIABLO expression to be associated with higher tumor
stemness for BRCA and UCS (Figure 5A and SF 11), and higher survival for
another (OV, SF 22). While neither DIABLO nor its corresponding protein Smac
are identified as top candidates in the mechanistic analyses, the findings of
the cBVS model are in line with existing literature, where DIABLO has been
proposed as a biomarker for breast and uterine cancers (SM 4.3).

#### HER2 as a Therapeutic Target in Gastrointestinal Cancers

The gene ERBB2 (protein HER2) has decisive evidence for all mechanistic models
in the pan-GI group, emerging as the top candidate (Figure 4C). In our cBVS
analyses, both ERBB2/HER2 are positively associated with stemness for ESCA and
STAD (SF 14-15). Past knowledge on HER2 identifies its increased expression as
an agent of poorer progression of gastric cancer, and it has been considered
as a therapeutic target, in line with our findings (SM 4.3).

##  5 Discussion and Future Work

We propose fiBAG, a hierarchical Bayesian framework to perform integration of
multi-omics data and outcomes. Using GPs, we quantify the functional roles of
proteogenomic biomarkers along three axes of interest (driver gene, driver
protein, cascading protein). Then, using a Bayesian variable selection
procedure with a calibrated spike-and-slab prior, we incorporate this
functional evidence in the outcome model to improve covariate selection and
effect size estimation. The framework offers novelty and utility in multiple
directions. First, it offers the user liberty in terms of multi-platform
integration, in the sense that depending on the data available, the
mechanistic and outcome layers can be appended/modified via
additions/alterations of the complete set of parameters of interest. For
example, we use DNA methylation and copy number alterations as upstream
covariates, but other information such as microRNAs/epigenetic factors (e.g.
histone modifications) could easily be incorporated. Second, the calibrated
spike-and-slab prior addresses the more general statistical question of
incorporating external information in a variable selection setting - by
updating the mapping function used to calibrate evidences to a prior
probability scale, the procedure can be adapted to settings where the
numerical evidences are in a different scale and/or sourced from different
models. cBVS is compared with multiple benchmarks via simulation studies in
synthetic and real data-based settings. The benchmarks include uncalibrated
Bayesian variable selection (EMVS), standard and grouped penalized regression
(LASSO and IPF-LASSO), and Bayesian variable selection incorporating external
information (graper and xtune). We use selection metrics like AUC and MCC to
evaluate performances across a broad spectrum of $n/p$ ratios. cBVS
outperforms all the uncalibrated methods across all $n/p$ settings, and
performs comparably with the calibrated methods in the high $n/p$ settings.
For low $n/p$, cBVS offers improvement in selection compared to other
calibrated methods (Figure 3).

Having established the improved utility of our method in an evidence-based
setting, we analyze pan-cancer multiomics data from TCGA, across a total of
four cancer groups (pan-gyne, pan-kidney, pan-squamous and pan-GI) and 14
cancers. Our real data analyses identify both known and novel associations at
cancer-specific and pan-cancer levels, such as the identification of the roles
of RPS6KA1 gene and p90RSK kinases in progression of gynecological cancers,
and the potential utility of the EGFR gene and protein as therapeutic targets
for ESCA where their expressions are negatively associated with survival
outcomes. Our framework is cancer-specific (each mechanistic and outcome model
is built for each cancer separately and then assessed in a pan-cancer fashion)
for several reasons. First of all, non-genomic covariates ($\mathbf{B}$ \-
e.g., demographic information such as gender or age, and clinical information
such as cancer stage) would potentially be entirely different across cancers,
and may have widely different scales of measure. Further, even the outcomes
may not always be normalized across cancers - while a scaled outcome like the
mRNAsi does not pose this problem, survival, for example, may differ
considerably across different cancers. Finally, the mechanistic models need to
be run separately for each gene and protein, and depending on the cancer of
interest, the functional roles of these genes and proteins are potentially
different, effectively rendering the measures of evidence to be cancer-
specific as well. However, our qualitative summarizations across relevant pan-
cancer groups provide valuable insights into the shared biology among these
cancers.

The flexible covariate-specific mechanistic modelling approach offers the
feasibility of utilizing platforms with larger expression pools and deeper
sequencing panels such as the National Cancer Institute Clinical Proteomic
Tumor Analysis Consortium (CPTAC,, 2022). Further, more cancers or cancer
groups can be included in the analysis pipeline as well, since the cBVS model
fitting is cancer-specific and we offer three options in decreasing order of
computational complexity (full MCMC, selection-only MCMC, and EMVS). All these
features make the whole procedure highly parallelizable (covariate level for
mechanistic models and cancer level for outcome models) - the computation
times for the simulation and real data settings described in this paper
reassure the utility of our pipeline in this direction. The median computation
times, rounded to the nearest whole time unit, for each procedure in our
integrative analyses are as follows: 2 minutes(mechanistic model), 5 hours
(cBVS full MCMC implementation), 32 minutes (cBVS selection-only MCMC
implementation), 4 minutes (cBVS EMVS implementation). All computations were
performed on an online cluster with each node having 16G memory.

An important novelty of our integrative analyses is the use of mRNA-based
cancer stemness indices as outcomes. Indeed, tumor stemness is a pertinent
determinant of cancer growth and prognostic outcomes, and only recently there
have been efforts to quantify stemness based on cellular signatures (Malta et
al., 2018). To the best of our knowledge, our study is the first to look at
potential associations of stemness with a pool of proteogenomic biomarkers in
an oncological context. A relevant question, however, would be whether the
associations identified are meaningful or artifacts of the construction of the
index. We further investigated the stemness indices and the inferred
associations to confirm the veracity of our findings. A summary of these
investigations is presented in SM 4.4.

This work indicates a number of potential avenues for future statistical
research. First, the flexibility in choosing the calibration function suggests
that the function can even potentially be data-driven (such as choosing the
tuning parameters via some pre-hoc analysis of the data). The cBVS framework
can also be adapted to settings slightly different from ours - such as
modeling other types of outcomes by using a nonlinear link function for the
mean of the outcome model, or using bivariate outcomes and so on.

To offer seamless interactive visualization of our integrative analysis
results, we have built an R Shiny app, the functionalities of which are
presented in SM 5. All our model and figure codes will be made available
publicly via the app. We hope that fiBAG will be a significant contribution to
current armamentarium of integrative models for multiomic data and can be
profitably employed to enhance prognostic/therapeutic utility in oncological
treatment and research as well as other disease domains.

## References

  * Abrahao-Machado and Scapulatempo-Neto, (2016) Abrahao-Machado, L. F. and Scapulatempo-Neto, C. (2016).  Her2 testing in gastric cancer: An update.  World Journal of gastroenterology, 22(19):4619.
  * Adashek et al., (2020) Adashek, J., Arroyo-Martinez, Y. M., Menta, A. K., Kurzrock, R., and Kato, S. (2020).  Therapeutic implications of epidermal growth factor receptor (egfr) in the treatment of metastatic gastric/gej cancer.  Frontiers in Oncology, 10:1312.
  * Adorno-Cruz et al., (2015) Adorno-Cruz, V., Kibria, G., Liu, X., Doherty, M., Junk, D. J., Guan, D., Hubert, C., Venere, M., Mulkearns-Hubert, E., Sinyuk, M., et al. (2015).  Cancer stem cells: targeting the roots of cancer, seeds of metastasis, and sources of therapy resistance.  Cancer research, 75(6):924–929.
  * Alizadeh et al., (2015) Alizadeh, A. A., Aranda, V., Bardelli, A., Blanpain, C., Bock, C., Borowski, C., Caldas, C., Califano, A., Doherty, M., Elsner, M., et al. (2015).  Toward understanding and exploiting tumor heterogeneity.  Nature medicine, 21(8):846–853.
  * Angel et al., (2020) Angel, P. W., Rajab, N., Deng, Y., Pacheco, C. M., Chen, T., Lê Cao, K.-A., Choi, J., and Wells, C. A. (2020).  A simple, scalable approach to building a cross-platform transcriptome atlas.  PLoS computational biology, 16(9):e1008219.
  * Arbiser, (2018) Arbiser, J. L. (2018).  Diablo: A double-edged sword in cancer?  Molecular Therapy, 26(3):678–679.
  * Arellano-Llamas et al., (2006) Arellano-Llamas, A., Garcia, F. J., Perez, D., Cantu, D., Espinosa, M., De la Garza, J. G., Maldonado, V., and Melendez-Zajgla, J. (2006).  High smac/diablo expression is associated with early local recurrence of cervical cancer.  BMC cancer, 6(1):1–10.
  * Arienti et al., (2019) Arienti, C., Pignatta, S., and Tesei, A. (2019).  Epidermal growth factor receptor family and its role in gastric cancer.  Frontiers in oncology, 9:1308.
  * Aydin et al., (2014) Aydin, K., Okutur, S. K., Bozkurt, M., Turkmen, I., Namal, E., Pilanci, K., Ozturk, A., Akcali, Z., Dogusoy, G., and Demir, O. G. (2014).  Effect of epidermal growth factor receptor status on the outcomes of patients with metastatic gastric cancer: A pilot study.  Oncology letters, 7(1):255–259.
  * Boonstra and Barbaro, (2020) Boonstra, P. S. and Barbaro, R. P. (2020).  Incorporating historical models with adaptive bayesian updates.  Biostatistics, 21(2):e47–e64.
  * Boonstra et al., (2015) Boonstra, P. S., Taylor, J. M., and Mukherjee, B. (2015).  Data-adaptive shrinkage via the hyperpenalized EM algorithm.  Statistics in biosciences, 7(2):417–431.
  * Boulesteix et al., (2017) Boulesteix, A.-L., De Bin, R., Jiang, X., and Fuchs, M. (2017).  Ipf-lasso: Integrative-penalized regression with penalty factors for prediction based on multi-omics data.  Computational and mathematical methods in medicine, 2017.
  * Bradfield et al., (2020) Bradfield, A., Button, L., Drury, J., Green, D. C., Hill, C. J., and Hapangama, D. K. (2020).  Investigating the role of telomere and telomerase associated genes and proteins in endometrial cancer.  Methods and protocols, 3(3):63.
  * Buccitelli and Selbach, (2020) Buccitelli, C. and Selbach, M. (2020).  mrnas, proteins and the emerging principles of gene expression control.  Nature Reviews Genetics, 21(10):630–644.
  * Califano and Alvarez, (2017) Califano, A. and Alvarez, M. J. (2017).  The recurrent architecture of tumour initiation, progression and drug sensitivity.  Nature reviews Cancer, 17(2):116–130.
  * Chatterjee et al., (2020) Chatterjee, D., Maitra, T., and Bhattacharya, S. (2020).  A short note on almost sure convergence of bayes factors in the general set-up.  The American Statistician, 74(1):17–20.
  * Chatterjee et al., (2016) Chatterjee, N., Chen, Y.-H., Maas, P., and Carroll, R. J. (2016).  Constrained maximum likelihood estimation for model calibration using summary-level information from external big data sources.  Journal of the American Statistical Association, 111(513):107–117.
  * Cheng et al., (2015) Cheng, C., Tseng, G., Ghosh, D., and Zhou, X. J. (2015).  From transcription factor binding and histone modification to gene expression: Integrative quantitative models.  Integrating Omics Data, page 380.
  * CPTAC, (2022) CPTAC (2022).  National cancer institute clinical proteomic tumor analysis consortium.  https://proteomics.cancer.gov/programs/cptac.  Accessed: 2022-12-22.
  * Daily et al., (2017) Daily, K., Sui, S. J. H., Schriml, L. M., Dexheimer, P. J., Salomonis, N., Schroll, R., Bush, S., Keddache, M., Mayhew, C., Lotia, S., et al. (2017).  Molecular, phenotypic, and sample-associated data to describe pluripotent stem cell lines and derivatives.  Scientific data, 4(1):1–10.
  * Dobrzycka et al., (2015) Dobrzycka, B., Mackowiak-Matejczyk, B., Terlikowska, K. M., Kulesza-Bronczyk, B., Kinalski, M., and Terlikowski, S. J. (2015).  Prognostic significance of pretreatment vegf, survivin, and smac/diablo serum levels in patients with serous ovarian carcinoma.  Tumor Biology, 36(6):4157–4165.
  * Dobrzycka et al., (2010) Dobrzycka, B., Terlikowski, S. J., Bernaczyk, P. S., Garbowicz, M., Niklinski, J., Chyczewski, L., and Kulikowski, M. (2010).  Prognostic significance of smac/diablo in endometrioid endometrial cancer.  Folia histochemica et cytobiologica, 48(4):678–681.
  * Duan et al., (2021) Duan, R., Gao, L., Gao, Y., Hu, Y., Xu, H., Huang, M., Song, K., Wang, H., Dong, Y., Jiang, C., et al. (2021).  Evaluation and comparison of multi-omics data integration methods for cancer subtyping.  PLoS computational biology, 17(8):e1009224.
  * Espinosa et al., (2021) Espinosa, M., Lizárraga, F., Vázquez-Santillán, K., Hidalgo-Miranda, A., Piña-Sánchez, P., Torres, J., García-Ramírez, R. A., Maldonado, V., Melendez-Zajgla, J., and Ceballos-Cancino, G. (2021).  Coexpression of smac/diablo and estrogen receptor in breast cancer.  Cancer Biomarkers, pages 1–18.
  * Finotello et al., (2020) Finotello, F., Calura, E., Risso, D., Hautaniemi, S., and Romualdi, C. (2020).  Multi-omic data integration in oncology.  Frontiers in oncology, 10:1768.
  * Fulda, (2013) Fulda, S. (2013).  Regulation of apoptosis pathways in cancer stem cells.  Cancer letters, 338(1):168–173.
  * George and McCulloch, (1997) George, E. I. and McCulloch, R. E. (1997).  Approaches for Bayesian variable selection.  Statistica Sinica, pages 339–373.
  * Gevaert et al., (2013) Gevaert, O., Villalobos, V., Sikic, B. I., and Plevritis, S. K. (2013).  Identification of ovarian cancer driver genes by using module network integration of multi-omics data.  Interface focus, 3(4):20130013.
  * Gravalos and Jimeno, (2008) Gravalos, C. and Jimeno, A. (2008).  Her2 in gastric cancer: a new prognostic factor and a novel therapeutic target.  Annals of oncology, 19(9):1523–1529.
  * Guo et al., (2019) Guo, L., Chen, Y., Luo, J., Zheng, J., and Shao, G. (2019).  Yap 1 overexpression is associated with poor prognosis of breast cancer patients and induces breast cancer cell growth by inhibiting pten.  FEBS Open Bio, 9(3):437–445.
  * He et al., (2015) He, C., Mao, D., Hua, G., Lv, X., Chen, X., Angeletti, P. C., Dong, J., Remmenga, S. W., Rodabaugh, K. J., Zhou, J., et al. (2015).  The hippo/yap pathway interacts with egfr signaling and hpv oncoproteins to regulate cervical cancer progression.  EMBO molecular medicine, 7(11):1426–1449.
  * Hinne et al., (2020) Hinne, M., Gronau, Q. F., van den Bergh, D., and Wagenmakers, E.-J. (2020).  A conceptual introduction to bayesian model averaging.  Advances in Methods and Practices in Psychological Science, 3(2):200–215.
  * Jennings et al., (2013) Jennings, E. M., Morris, J. S., Carroll, R. J., Manyam, G. C., and Baladandayuthapani, V. (2013).  Bayesian methods for expression-based integration of various types of genomics data.  EURASIP Journal on Bioinformatics and Systems Biology, 2013(1):1–11.
  * Kandel et al., (2014) Kandel, C., Leclair, F., Bou-Hanna, C., Laboisse, C. L., and Mosnier, J.-F. (2014).  Association of her1 amplification with poor prognosis in well differentiated gastric carcinomas.  Journal of clinical pathology, 67(4):307–312.
  * Kaplan and Lock, (2017) Kaplan, A. and Lock, E. F. (2017).  Prediction with dimension reduction of multiple molecular data sources for patient survival.  Cancer informatics, 16:1176935117718517.
  * Kass and Raftery, (1995) Kass, R. E. and Raftery, A. E. (1995).  Bayes factors.  Journal of the american statistical association, 90(430):773–795.
  * Lawrence et al., (2007) Lawrence, N. D., Sanguinetti, G., and Rattray, M. (2007).  Modelling transcriptional regulation using gaussian processes.  Advances in Neural Information Processing Systems, 19:785.
  * Lee et al., (2019) Lee, J., Franovic, A., Shiotsu, Y., Kim, S. T., Kim, K.-M., Banks, K. C., Raymond, V. M., and Lanman, R. B. (2019).  Detection of erbb2 (her2) gene amplification events in cell-free dna and response to anti-her2 agents in a large asian cancer patient cohort.  Frontiers in oncology, 9:212.
  * Litovkin et al., (2014) Litovkin, K., Joniau, S., Laerut, E., Laenen, A., Gevaert, O., Spahn, M., Kneitz, B., Isebaert, S., Haustermans, K., Beullens, M., et al. (2014).  Methylation of pitx2, hoxd3, rassf1 and tdrd1 predicts biochemical recurrence in high-risk prostate cancer.  Journal of cancer research and clinical oncology, 140(11):1849–1861.
  * Liu et al., (2007) Liu, D., Lin, X., and Ghosh, D. (2007).  Semiparametric regression of multidimensional genetic pathway data: Least-squares kernel machines and linear mixed models.  Biometrics, 63(4):1079–1088.
  * Malaguti et al., (2015) Malaguti, P., D’Aloia, M. M., and Alimandi, M. (2015).  The erbb2 receptor in gastric cancer: the quick-change artist.  Translational Gastrointestinal Cancer, 4(4):282–293.
  * Malta et al., (2018) Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., Kamińska, B., Huelsken, J., Omberg, L., Gevaert, O., et al. (2018).  Machine learning identifies stemness features associated with oncogenic dedifferentiation.  Cell, 173(2):338–354.
  * Mamoor, (2021) Mamoor, S. (2021).  Over-expression of ribosomal protein s6 kinase a1 in human endometrial cancer.  OSF Preprints.
  * Micchelli et al., (2006) Micchelli, C. A., Xu, Y., and Zhang, H. (2006).  Universal kernels.  Journal of Machine Learning Research, 7(12).
  * Morris and Baladandayuthapani, (2017) Morris, J. S. and Baladandayuthapani, V. (2017).  Statistical contributions to bioinformatics: Design, modelling, structure learning and integration.  Statistical modelling, 17(4-5):245–289.
  * Plummer et al., (2016) Plummer, M., Stukalov, A., Denwood, M., and Plummer, M. M. (2016).  Package ‘rjags’.  Vienna, Austria.
  * Rattray et al., (2019) Rattray, M., Yang, J., Ahmed, S., and Boukouvalas, A. (2019).  Modelling gene expression dynamics with gaussian process inference.  Handbook of Statistical Genomics: Two Volume Set, pages 879–20.
  * Richardson et al., (2016) Richardson, S., Tseng, G. C., and Sun, W. (2016).  Statistical methods in integrative genomics.  Annual review of statistics and its application, 3:181–209.
  * Ročková and George, (2014) Ročková, V. and George, E. I. (2014).  EMVS: The EM approach to Bayesian variable selection.  Journal of the American Statistical Association, 109(506):828–846.
  * Shareefi et al., (2020) Shareefi, G., Turkistani, A. N., Alsayyah, A., Kussaibi, H., Had, M. A., and Alkharsah, K. R. (2020).  Pathway-affecting single nucleotide polymorphisms (snps) in rps6ka1 and mbip genes are associated with breast cancer risk.  Asian Pacific Journal of Cancer Prevention: APJCP, 21(7):2163.
  * Solvang et al., (2011) Solvang, H. K., Lingjærde, O. C., Frigessi, A., Børresen-Dale, A.-L., and Kristensen, V. N. (2011).  Linear and non-linear dependencies between copy number aberrations and mrna expression reveal distinct molecular pathways in breast cancer.  BMC bioinformatics, 12(1):1–12.
  * Song et al., (2019) Song, X., Ji, J., Gleason, K. J., Yang, F., Martignetti, J. A., Chen, L. S., and Wang, P. (2019).  Insights into impact of dna copy number alteration and methylation on the proteogenomic landscape of human ovarian cancer via a multi-omics integrative analysis.  Molecular & Cellular Proteomics, 18(8):S52–S65.
  * Stephens and Balding, (2009) Stephens, M. and Balding, D. J. (2009).  Bayesian statistical methods for genetic association studies.  Nature Reviews Genetics, 10(10):681–690.
  * Subramanian et al., (2020) Subramanian, I., Verma, S., Kumar, S., Jere, A., and Anamika, K. (2020).  Multi-omics data integration, interpretation, and its application.  Bioinformatics and biology insights, 14:1177932219899051.
  * Sun et al., (2021) Sun, G., Rong, D., Li, Z., Sun, G., Wu, F., Li, X., Cao, H., Cheng, Y., Tang, W., and Sun, Y. (2021).  Role of small molecule targeted compounds in cancer: Progress, opportunities, and challenges.  Frontiers in Cell and Developmental Biology, page 2043.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * Torchiaro et al., (2016) Torchiaro, E., Lorenzato, A., Olivero, M., Valdembri, D., Gagliardi, P. A., Gai, M., Erriquez, J., Serini, G., and Di Renzo, M. F. (2016).  Peritoneal and hematogenous metastases of ovarian cancer cells are both controlled by the p90rsk through a self-reinforcing cell autonomous mechanism.  Oncotarget, 7(1):712.
  * Tseng et al., (2015) Tseng, G., Ghosh, D., and Zhou, X. J. (2015).  Integrating omics data.  Cambridge University Press.
  * Tsujiura et al., (2014) Tsujiura, M., Mazack, V., Sudol, M., Kaspar, H. G., Nash, J., Carey, D. J., and Gogoi, R. (2014).  Yes-associated protein (yap) modulates oncogenic features and radiation sensitivity in endometrial cancer.  PloS one, 9(6):e100974.
  * Tu et al., (2015) Tu, Z., Zhang, B., and Zhu, J. (2015).  Network integration of genetically regulated gene expression to study complex diseases.  Integrating Omics Data, 88:88–109.
  * Velten and Huber, (2021) Velten, B. and Huber, W. (2021).  Adaptive penalization in high-dimensional regression and classification with external covariates using variational bayes.  Biostatistics, 22(2):348–364.
  * Wang et al., (2020) Wang, D., He, J., Dong, J., Meyer, T. F., and Xu, T. (2020).  The hippo pathway in gynecological malignancies.  American journal of cancer research, 10(2):610.
  * Wang et al., (2013) Wang, W., Baladandayuthapani, V., Morris, J. S., Broom, B. M., Manyam, G., and Do, K.-A. (2013).  ibag: integrative bayesian analysis of high-dimensional multiplatform genomics data.  Bioinformatics, 29(2):149–159.
  * Zeng et al., (2021) Zeng, C., Thomas, D. C., and Lewinger, J. P. (2021).  Incorporating prior knowledge into regularized regression.  Bioinformatics, 37(4):514–521.

Supplementary Materials

##  1 fiBAG Method

###  1.1 Bayes Factor Expression for Gaussian Process Mechanistic Models

We describe the Bayes factor derivation for the Gaussian process specification
corresponding to the driver gene mechanistic model in this section. The
expressions for the other models can be derived similarly. To recall the
notations introduced, main manuscript, the $j^{\textrm{th}}$ driver gene model
is written as
$G_{ij}=f_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T})+e_{1ij}$, where
we assume $e_{1ij}\overset{\textrm{iid}}{\sim}\textrm{N}(0,\tau_{1j}^{2})$.
Let us also denote
$f_{1j}^{(i)}=f_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T})$ for all
$i$. Then, the GP prior on $f_{1j}$ is placed as follows:

GP prior:
$\displaystyle(f_{1j}^{(1)},\ldots,f_{1j}^{(n)})^{T}\sim\mathbf{N}(\mathbf{0},\mathbf{K}_{1j}),$
Covariance matrix:
$\displaystyle\mathbf{K}_{1j(i,k)}=K_{1j}((\mathbf{C}_{ij}^{T},\mathbf{M}_{ij}^{T})^{T},(\mathbf{C}_{kj}^{T},\mathbf{M}_{kj}^{T})^{T}),$
Kernel function: $\displaystyle
K_{1j}(\mathbf{u},\mathbf{v})=g\tau_{1j}^{2}\textrm{exp}\big{(}-\frac{||\mathbf{u}-\mathbf{v}||^{2}}{\lambda_{1j}^{2}}\big{)}.$

The hyperpriors specify $\tau_{1j}^{2}\sim\textrm{Inverse-
Gamma}(\frac{\nu_{01j}}{2},\frac{\nu_{01j}\tau_{01j}^{2}}{2})$ and
$\lambda_{1j}\sim\textrm{exp}(\lambda_{01j})$. For all our analyses, we set
$g=n$. To test the hypothesis of interest in this setting, we need to compare
this model with an intercept model, specified by $G_{ij}=\mu_{j}+e_{1ij}$. We
assume a Zellner’s g-prior on $\mu_{j}$ as
$N\left(0,\tau_{1j}^{2}\frac{g}{1+g}\right)$. We use $g=n$ for all our models.
Under this specification, the marginal likelihood is the following.

$\displaystyle\int_{0}^{\infty}\int_{-\infty}^{\infty}\textrm{P}(\mathbf{G}_{\bullet
j}|\mu_{j},\tau_{1j}^{2})\textrm{P}(\mu_{j}|g)\textrm{P}(\tau_{1j}^{2}|\nu_{01j},\tau_{01j})d\mu_{j}d\tau_{1j}^{2}\
$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{-\infty}^{\infty}(2\pi)^{-n/2}(\tau_{1j}^{2})^{-n/2}\exp(-\sum_{i=1}^{n}(G_{ij}-\mu_{j})^{2}/2\tau_{1j}^{2}).(2\pi)^{-1/2}(\tau_{1j}^{2}\frac{g}{1+g})^{-1/2}$
$\displaystyle\exp(-\mu_{j}^{2}/(2\tau_{1j}^{2}\frac{g}{1+g})).\frac{1}{\Gamma(\nu_{01j}/2)}(\nu_{01j}\tau_{01j}^{2}/2)^{\nu_{01j}/2}(\tau_{01j}^{2})^{-\nu_{01j}/2-1}\exp(-\nu_{01j}\tau_{01j}^{2}/2\tau_{1j}^{2})d\mu_{j}d\tau_{1j}^{2}$
$\displaystyle=$ $\displaystyle
C.\int_{0}^{\infty}(\tau_{1j}^{2})^{-(n+\nu_{01j}+3)/2}\exp(-\nu_{01j}\tau_{01j}^{2}/2\tau_{1j}^{2})\int_{-\infty}^{\infty}\exp\big{[}-\big{\\{}(n-\frac{1+g}{g})\mu_{j}^{2}-2\mu_{j}\sum_{i=1}^{n}G_{ij}\big{\\}}/2\tau_{1j}^{2}\big{]}d\mu_{j}$
$\displaystyle\exp(-\sum_{i=1}^{n}G_{ij}^{2}/2\tau_{1j}^{2})d\tau_{1j}^{2}$
$\displaystyle=$ $\displaystyle
C.\int_{0}^{\infty}(\tau_{1j}^{2})^{-(n+\nu_{01j}+3)/2}\exp\big{\\{}-(\nu_{01j}\tau_{01j}^{2}+\sum_{i=1}^{n}G_{ij}^{2})/2\tau_{1j}^{2}\big{\\}}\sqrt{2\tau_{1j}^{2}\pi/(n-1-1/g)}$
$\displaystyle\exp\big{[}(\sum_{i=1}^{n}G_{ij}/2\tau_{1j}^{2})^{2}\big{\\{}2\tau_{1j}^{2}/(n-1-1/g)\big{\\}}\big{]}d\tau_{1j}^{2}$
$\displaystyle=$ $\displaystyle
C^{{}^{\prime}}.\int_{0}^{\infty}(\tau_{1j}^{2})^{-(n+\nu_{01j})/2-1}\exp(-\frac{\nu_{01j}\tau_{01j}^{2}+\sum_{i=1}^{n}G_{ij}^{2}-(\sum_{i=1}^{n}G_{ij})^{2}/(n-1-1/g)}{2\tau_{1j}^{2}})d\tau_{1j}^{2}$
$\displaystyle=$ $\displaystyle
C^{{}^{\prime}}.\Gamma((n+\nu_{01j})/2).\big{\\{}\frac{\nu_{01j}\tau_{01j}^{2}+\sum_{i=1}^{n}G_{ij}^{2}-(\sum_{i=1}^{n}G_{ij})^{2}/(n-1-1/g)}{2}\big{\\}}^{\big{\\{}-(n+\nu_{01j})/2\big{\\}}}.$

Then, for the intercept model, the expression for the log-marginal likelihood
simplifies to the following:
$a^{{}^{\prime}}_{n}+b_{n}\ln\Big{(}a+\sum_{i=1}^{n}G_{ij}^{2}-\frac{(\sum_{i=1}^{n}G_{ij})^{2}}{c_{n}}\Big{)}$,
where $a^{{}^{\prime}}_{n},b_{n},c_{n}$, and $a$ are constants dependent on
data dimensions and hyperparameters. Following the same procedures, we can
derive the log-marginal likelihood for the Gaussian process model. The only
difference in the calculation above occurs due to the fact that the joint
distribution of $\mathbf{G}_{\bullet j}$ does not have a diagonal covariance
matrix anymore. Therefore, the terms involving the data inside the double
integral does not anymore simplify to functions of $\sum_{i=1}^{n}G_{ij}^{2}$
and $(\sum_{i=1}^{n}G_{ij})^{2}$. The integral with respect to
$d\tau_{1j}^{2}$ can still be simplified as before, but the one involving
$d\lambda_{01j}$ does not have a closed-form solution. The final expression
for the log-marginal likelihood is then:
$a^{{}^{\prime\prime}}_{n}+\ln\int_{0}^{\infty}\exp(-\lambda_{01j}\lambda_{1j})\frac{2^{b_{n}}|\mathbf{K}_{1j}/\tau_{1j}^{2}+I|^{-\frac{1}{2}}}{\big{\\{}\mathbf{G}_{\bullet
j}^{T}(\mathbf{K}_{1j}/\tau_{1j}^{2}+I)\mathbf{G}_{\bullet
j}+a)\big{\\}}^{b_{n}}}d\lambda_{1j}$. Subtracting the previous log-marginal
likelihood from this expression provides us the expression of the log-marginal
likelihood as presented, main manuscript.

###  1.2 Choice of the Evidence Calibration Function

We perform the following steps to obtain the shape parameters of the final
Beta prior as $\mathcal{F}(\boldsymbol{\mathcal{E}}_{j})=\mathcal{F}(s_{j})$
and $1/\mathcal{F}(\boldsymbol{\mathcal{E}}_{j})=1/\mathcal{F}(s_{j})$. Note
that due to this formulation, the prior mean for $\omega_{j}$ takes the form
$\mathcal{F}(s_{j})\big{\\{}\mathcal{F}(s_{j})+1/\mathcal{F}(s_{j})\big{\\}}^{-1}$,
and the corresponding variance is
$\big{[}\big{\\{}\mathcal{F}(s_{j})+1/\mathcal{F}(s_{j})\big{\\}}^{2}\big{\\{}\mathcal{F}(s_{j})+1/\mathcal{F}(s_{j})+1\big{\\}}\big{]}^{-1}$.
We first ensure that all small positive and non-positive lBFs from the
mechanistic model are effectively truncated to zero (no evidence, to be mapped
to a uniform prior), by setting $s^{*}_{j}=\max(s_{j},10^{-6})$. Then, our
mapping function evaluates
$\mathcal{G}(s_{j})=\frac{1}{2}\big{[}\big{\\{}1+(s^{*}_{j}/3)^{-2.75}\big{\\}}^{-1}+1\big{]}$.
Finally, we compute
$\mathcal{F}(s_{j})=(2\mathcal{G}(s_{j}))^{4}=16\mathcal{G}^{4}(s_{j})$ to
ensure sharp increase in prior mean probability of inclusion when lBF
increases from 1 to 2. Effectively, $\mathcal{F}(\bullet)$ is a composition of
three functions.

  * •

A four-parameter logistic
($s^{*}_{j}\rightarrow\frac{1}{1+(s^{*}_{j}/3)^{-2.75}}$) allowing us to
control the shape and scale of the prior probability calibration curve along
with its asymptotic behavior and lower bound.

  * •

A linear map ($x\rightarrow\frac{1}{2}(x-1)$) that takes $[0,1]$ to $[0.5,1]$
enabling us to make the prior noninformative when the lBF is small or
negative.

  * •

A polynomial map ($c\rightarrow(2c)^{4}$) allowing additional control over the
steepness of the resulting prior probability calibrated means.

The tuning parameters with values 3 and 2.75 are chosen via computational
checks across ranges of $[2,4]$ and $[2,3]$, respectively. For the
specification chosen, the prior mean probability of inclusion for a covariate
ranges from $0.5$ (when the lBF is $0$ or negative) to approximately $1$ (when
the lBF is large, say $>5$).

###  1.3 Outcome Model Fitting Procedures

As described in the main manuscript, the complete set of parameters to be
estimated by cBVS is
$(\boldsymbol{\beta}^{T},\boldsymbol{\gamma}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$.
$\boldsymbol{\beta}^{T}$ provides estimates of the effect sizes of the
proteogenomic covariates on the outcome, and the
$(\boldsymbol{\gamma}^{T},\boldsymbol{\omega}^{T})^{T}$ guides their
selection. The parameter estimation is focused on the posterior of
$(\boldsymbol{\beta}^{T},\boldsymbol{\gamma}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$.
For large proteogenomic panels, this posterior will be computationally
resource-intensive to directly sample from. The question, therefore, is of a
trade-off between computational simplicity and estimation accuracy. Due to
this reason, we offer three implementations of cBVS in increasing order of
computational efficiency – a standard Markov chain Monte Carlo (MCMC) using
Gibbs sampler to sample from the complete posterior, a selection-only MCMC to
sample from the marginal posterior of $\boldsymbol{\gamma}$, and an
expectation-maximization based variable selection (EMVS) procedure to
approximate the posterior modes of the parameters. Briefly, the selection-only
MCMC focuses on estimating $\boldsymbol{\gamma}$ first and then estimates
$\boldsymbol{\beta}$ using a Bayesian model averaging-type procedure (Hinne et
al., 2020), and the EMVS sacrifices the full posterior along with error
estimates to achieve fast point estimation (Ročková and George, 2014). The
exact details of each are as follows.

  1. 1.

Full MCMC: The simplest approach is to perform a complete Markov chain Monte
Carlo (MCMC) procedure to sample from the joint posterior – we utilize a Gibbs
sampler for this. For both continuous and survival outcomes, we use the rjags
(Plummer et al., 2016) package in R. The rjags model descriptions are
available in the code archive included in this submission.

  2. 2.

Selection-only MCMC: One way to mitigate the time- and resource-intensiveness
of the full MCMC is to focus on $\boldsymbol{\gamma}$ only - by integrating
$(\boldsymbol{\beta}^{T},\boldsymbol{\omega}^{T},\sigma)^{T}$ out of the joint
posterior. This results in $\Pi(\boldsymbol{\gamma}|\textrm{Data,
Hyperparameters})$, to be then approximated via MCMC. The final estimates for
$\boldsymbol{\omega}$ are computed by taking the average of the traversed MCMC
path for $\boldsymbol{\gamma}$ post burn-in. The components of
$\boldsymbol{\beta}$ are estimated using a Bayesian model averaging-type
procedure (Hinne et al., 2020), taking a convex combination of draws from
their conditional posterior at every step, weighted proportionally to the
negative log posterior evaluated at that step. Since this MCMC only performs a
search across a lattice space of size $2^{p}$, it results in significant
improvements in computation times. The overview of the selection-only MCMC
procedure is presented in Supplementary Algorithm 1.1. The codes are available
in the code archive included in this submission.

  3. 3.

E-M based variable selection (EMVS): The fastest alternative to a full MCMC is
to sacrifice approximating the full posterior (along with error margins) and
to only focus on point estimates of the parameters of interest. For this
purpose, we adapt the EMVS procedure by Ročková and George, (2014) for our
continuous and survival settings. The EMVS procedure estimates the posterior
mode instead of approximating the full posterior, resulting in faster
iterations. The model codes for the implementation of EMVS in our setting are
available in the code archive included in this submission.

Supplementary Algorithm 1.1. Selection-only MCMC algorithm 0. Initiate the
MCMC with $\boldsymbol{\gamma}=\boldsymbol{\gamma}_{\textrm{old}}$, where
$\gamma_{\textrm{old,}j}\sim
Ber(\frac{a(s_{j})}{a(s_{j})+\frac{1}{a(s_{j})}}),\forall j$. 1. With
probability $1/3$ each, select one of the following steps and perform the
corresponding task. 1.1. Add: Randomly pick one of the
$\boldsymbol{\gamma}_{\textrm{old}}$ indices with value $0$, change it to $1$,
and call the resulting vector $\boldsymbol{\gamma}_{\textrm{new}}$. If there
is no index $j$ where $\gamma_{\textrm{old,}j}=0$, perform 1.2. 1.2. Delete:
Randomly pick one of the $\boldsymbol{\gamma}_{\textrm{old}}$ indices with
value $1$, change it to $0$, and call the resulting vector
$\boldsymbol{\gamma}_{\textrm{new}}$. If there is no index $j$ where
$\gamma_{\textrm{old,}j}=1$, perform 1.1. 1.3. Swap: Randomly pick one of the
$\boldsymbol{\gamma}_{\textrm{old}}$ indices with value $0$, and
independently, randomly pick one of the $\boldsymbol{\gamma}_{\textrm{old}}$
indices with value $1$. Swap the values of these two indices, and call the
resulting vector $\boldsymbol{\gamma}_{\textrm{new}}$. If there is no index
$j$ where $\gamma_{\textrm{old,}j}=1$, perform 1.1. If there is no index $j$
where $\gamma_{\textrm{old,}j}=0$, perform 1.2. 2. Compute
$\log(p^{*}):=\min\\{0,\log(\Pi(\boldsymbol{\gamma}_{\textrm{new}}|.))-\log(\Pi(\boldsymbol{\gamma}_{\textrm{old}}|.))\\}$.
The proposed $\boldsymbol{\gamma}_{\textrm{new}}$ is then accepted with
probability $p^{*}$. 3. Iterate between 1-2 until convergence.

###  1.4 FDR Control for Variable Selection using Posterior Inclusion
Probabilities

The MCMC/EMVS procedures provide us with posterior inclusion probabilities
(PIPs) $\hat{\omega_{j}}$ and regression coefficient estimates
$\hat{\beta_{j}}$ for each covariate in the model. A cut-off on the PIPs is
computed using a false discovery rate adjustment procedure at a specified
level of significance, treating the $1-\hat{\omega_{j}}$ as p-value type
quantities. Suppose the estimated posterior probabilities of inclusion are
denoted by $\\{\hat{\omega}_{j},j\in\\{1,\ldots,q_{g}+q_{p}\\}\\}$. Then, we
define p-value type quantities $p_{j}=1-\hat{\omega}_{j},\forall
j\in\\{1,\ldots,q_{g}+q_{p}\\}$, and sort them in the increasing order of
magnitude as $\\{p^{*}_{j},j\in\\{1,\ldots,q_{g}+q_{p}\\}\\}$, with the
understanding that for each $j$, $p^{*}_{j}=p_{k_{j}}$ for some
$k_{j}\in\\{1,\ldots,q_{g}+q_{p}\\}$. Let the cumulative sums of these ordered
quantities be denoted by $r_{j}=\sum_{l=1}^{j}p^{*}_{j}$ for each $j$. Let
$j^{*}:=\min\\{j:r_{j}\geq\alpha\\}$, where $\alpha$ is a pre-specified level
for the false discovery rate control. Then we infer that the covariates with
indices $k_{1},\ldots,k_{j^{*}}$ are selected in the outcome model.

##  2 Simulation Studies

###  2.1 Data Generation for Simulation 1

To compare performances across a grid of varying sample size/number of
covariates ($n/p$) ratios, we fix the number of proteogenomic biomarkers
generated in the simulated datasets at $p=200$ and vary the sample size across
$n=50,100,200,400,800$, covering a range of $n/p=1/4$ to $4$. For simplicity,
we assume that there is only one upstream covariate for each biomarker. 100
replicates are generated for each n. We follow the steps described below to
generate one such replicate.

  * •

Upstream covariates $\boldsymbol{U}_{n\times p}$ (comparable to copy
number/methylation data in the biological setting) are generated such that
each $U_{ij}\overset{\textrm{iid}}{\sim}\textrm{N}(0,1)$.

  * •

Let $\boldsymbol{X}_{n\times p}$ denote the design matrix of the proteogenomic
expression data for the outcome model. Then the generating model for the
$j^{\textrm{th}}$ expression is
$X_{ij}\overset{\textrm{ind}}{\sim}\textrm{N}(\xi_{j}U_{ij},1)$, $\forall
i\in\\{1,\ldots,n\\}$. The mechanistic effect size $\xi_{j}$ controls the
level of evidence reflected by an lBF for the mechanistic model of
$\boldsymbol{X}_{\bullet j}$ on $\boldsymbol{U}_{\bullet j}$. The
correspondence between values of $\xi_{j}$ and the four levels of lBF: no
evidence, substantial, strong, and decisive, is established numerically.

Among the 200 $\boldsymbol{X}_{\bullet j}$s, the first 60 are distributed into
four groups of 15 each, in the order of no evidence, substantial, strong, and
decisive, followed by two groups of five at the levels strong and decisive,
respectively. Those first 60 are distributed further into varying outcome
effect sizes in the later steps to cover the complete range of combinations of
prior evidence and effect size. The latter two groups will be assigned no true
effects to include an in-built checkpoint for false positive evidences.

  * •

Among each group of 15 for the four levels of prior functional evidence, we
put three groups of five $\boldsymbol{X}_{\bullet j}$s with effect sizes
$\beta_{j}$s generated respectively at the low ($U(0,0.2)$), medium
($U(0.4,0.6)$), and high ($U(0.9,1.1)$) levels. All the other $\beta_{j}$s are
set $=0$.

This results in 12 groups of covariates of size five each, because there are
12 possible combinations from the three levels of effect sizes and the four
levels of evidences. Additionally, we have two groups of size five each with
strong/decisive level of evidence but no true effect. All other 130 covariates
have no evidence and no true effect.

  * •

Finally, we generate the continuous outcome data as
$\mathbf{Y}_{i}\overset{\textrm{ind}}{\sim}\mathcal{N}_{p}(\boldsymbol{\beta}^{T}\mathbf{X}_{i\bullet},1)$.

###  2.2 Comparing Linear Models and Gaussian Processes in Capturing Non-
linear Evidence

We perform a simple simulation study to compare the performances of Gaussian
process models and linear models in capturing significance based on data
generated from varying orders of non-linear associations. Hypothetically, with
increasing proportion of non-linearity in the generating model, the Bayes
factors from the Gaussian process models should be better able to quantify
significance than those from the linear models. To perform this, our
generating models are defined as the following.

  * •

For a given sample size $n=100$ and number of covariates $p=5$, we generate
$\mathbf{X}_{n\times p}$ where each $X_{ij}\overset{iid}{\sim}U(0,1)$.

  * •

The fully nonlinear model is defined as the following.

$Y_{i}\sim
N(10\cos(X_{i1})-15X_{i2}^{2}+10\exp(-X_{i3})X_{i2}-8\sin(X_{i3})\cos(X_{i4})+20X_{i1}X_{i5},1)$,
i.e., $Y_{i}\sim N(\sum_{j=1}^{p}\beta_{j}f_{ij},1)$. For
$l\in\\{0,1,\ldots,5\\}$, a $20l\%$ nonlinear model is then generated from the
following distribution.

$Y_{i}\sim
N(\sum_{j=1}^{l}\beta_{j}f_{ij}+\sum_{j=l+1}^{p}\beta_{j}X_{ij},1)$.

  * •

For each non-linearity scenario, 100 independent replicates are generated. For
each replicate, the Bayes factors corresponding to the Gaussian process model
(as described in Section 2.2, main manuscript) and the linear model are
computed.

As seen in Supplementary Figure 2, for 0% and 20% levels of non-linearity, the
median lBF from the linear models is larger than that from the Gaussian
processes, whereas this pattern is reversed for all the higher levels of non-
linearity. Noticeably, the difference between the two sets of lBFs increase
steadily with increasing non-linearity. This supports our claim that the
Gaussian process models are better equipped in identifying evidence based on
data originating from nonlinear generating models.

###  2.3 Comparing Coefficient Estimation from MCMC and EMVS Algorithms

To perform a comparison between the estimated coefficients from the full MCMC
and EMVS-based cBVS implementations, we used the data generated in Simulation
1. For each sample size and each replicate in Simulation 1, we estimate the
coefficients using both procedures. For a given pool of $p$ covariates
($p=200$) in Simulation 1, let us call these two sets of coefficient estimates
$\hat{\boldsymbol{\beta}}_{M}$ and $\hat{\boldsymbol{\beta}}_{E}$, each having
dimension $p\times 1$. For each replicate within each sample size scenario
then, we compute the mean squared difference between these as
$||\hat{\boldsymbol{\beta}}_{M}-\hat{\boldsymbol{\beta}}_{E}||^{2}/p$. These
quantities are summarized via box and violin plots in Supplementary Figure 3.

As can be seen in Supplementary Figure 3, the highest mean squared difference
between the two sets of estimates for the smallest sample size ($n=50$) is
lower than $1$. The mean difference across $100$ replicates decreases
consistently with increase in the sample size. In particular, all differences
are in the order (or less) of $10^{-1}$ for the sample sizes $\geq p=200$.
These summaries indicate that the compromise in approximating the full
posterior in lieu of estimating the posterior mode via EMVS does not affect
the estimation performance of cBVS substantially.

##  3 Integrative Proteogenomic Analyses

###  3.1 List of TCGA Cancer Groups and Cancers

  * •

Pan-gyne: breast invasive carcinoma (BRCA), cervical squamous cell carcinoma
and endocervical adenocarcinoma (CESC), ovarian serous cystadenocarcinoma
(OV), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma
(UCS).

  * •

Pan-kidney: kidney chromophobe (KICH), kidney renal clear cell carcinoma
(KIRC), kidney renal papillary cell carcinoma (KIRP).

  * •

Pan-squamous: esophageal carcinoma (ESCA, squamous), head and neck squamous
cell carcinoma (HNSC), lung squamous cell carcinoma (LUSC).

  * •

Pan-gastrointestinal (pan-GI): colon and rectum adenocarcinoma (CORE),
esophageal carcinoma (ESCA, adeno), stomach adenocarcinoma (STAD).

###  3.2 Survival Outcome Model Results

We summarize the survival cBVS outputs for each TCGA cancer using plots
showing $-log_{10}(1-\hat{\omega})$ on the y-axis and $\hat{\beta}$ on the
x-axis for each cancer, in Supplementary Figures 20-33. Further, we present
bar diagrams exhibiting the lBFs for the selected genes and proteins and
histograms with density plots of lBFs for the predictors not selected. These
figures are available in the shiny app for interactive exploration. Here,
similar to the stemness outcome model results, main manuscript, we briefly
discuss the results from the BRCA and CORE analyses.

For BRCA, no gene is selected using the 10% FDR cut-off (Supplementary Figure
20). The only protein selected is Collagen VI, which has decisive functional
evidence. Again, a vast majority of the non-selected predictors had little to
no functional evidence. For CORE, no gene/protein is selected at the 10% FDR
threshold (Supplementary Figure 31).

###  3.3 Past Evidences Related to Biological Findings

#### RPS6KA1 gene and p90RSK kinases in gynecological cancer progression

The gene RPS6KA1 (protein p90RSK) has decisive evidence for all but two
mechanistic models in the pan-gyne cancers, as seen in Figure 5A, main
manuscript. The gene has been known to be differentially expressed in
endometrial cancer tissue as opposed to benign endometrial tissue (Mamoor,
2021). Specifically, it is known to have a favorable prognostic effect on
clinical outcomes in endometrial cancers (Bradfield et al., 2020). The p90RSK
protein is selected in the survival cBVS model for OV (Supplementary Figure
22). p90RSK has been known to impact metastatic seeding of ovarian cancer
cells, effecting the invasiveness of the cancer via activating a self-
reinforcing cell autonomous circuit (Torchiaro et al., 2016). The RPS6KA1 gene
is also significantly associated with the increased risk of breast cancer
(Shareefi et al., 2020).

#### YAP1 gene and YAP proteins in gynecological malignancies

The gene YAP1 (proteins YAP, YAPPS127) has decisive evidence in all
mechanistic models for all gynecological cancers except UCS (smallest sample
size in the group, $n=57$), as seen in Figure 5A, main manuscript. YAP is a
crucial agent impacting gynecological cancers. As a transcriptional co-
activator within the Hippo pathway, over-activation of YAP leads to
uncontrolled cell growth and malignant transformation in gynecological
malignancies, including cervical, ovarian, and endometrial cancers (Wang et
al., 2020). Further, YAP expression is associated with a poor prognosis for
gynecological cancers - activation of YAP induces cancer cell proliferation
and migration (breast: Guo et al., (2019), cervical: He et al., (2015),
endometrial: Tsujiura et al., (2014)). This aligns with our identification of
YAP1 as negatively associated with stemness for BRCA (Figure 6A, main
manuscript) and with survival for OV and CESC (Supplementary Figures 21-22).

#### DIABLO gene as a marker of gynecological tumors

Another interesting candidate emerging from the mechanistic models is the gene
DIABLO (protein Smac), which has been proposed as a biomarker for
gynecological tumors, so far with little knowledge about its cellular
mechanism. A recent study shows some evidence in favor of a positive
association of Smac/DIABLO expression levels with estrogen receptor positivity
in breast cancer (Espinosa et al., 2021). Our cBVS analyses identify higher
DIABLO expression to be associated with higher tumor stemness for two
gynecological cancers (BRCA and UCS, Figure 6, main manuscript), which is in
line with prior knowledge (Arellano-Llamas et al., 2006; Arbiser, 2018). On
the other hand, DIABLO is associated with higher survival for OV
(Supplementary Figure 22), supported by earlier evidence - higher expression
of DIABLO is a good prognostic sign for ovarian and endometrial cancers
(Dobrzycka et al., 2010, 2015). However, neither DIABLO nor its corresponding
protein Smac are identified as top candidates in the mechanistic analyses.

#### HER2 as a therapeutic target in gastrointestinal cancers

The gene ERBB2 (protein HER2) has decisive evidence for all mechanistic models
in the pan-GI group, emerging as the top candidate (Figure 5C, main
manuscript). HER2 overexpression has been known to be a frequent molecular
abnormality in gastric cancers via gene amplification (Gravalos and Jimeno,
2008). HER2 has also been considered as a molecular therapeutic target for
patients with advanced gastric cancer (Abrahao-Machado and Scapulatempo-Neto,
2016). In our cBVS analyses, both ERBB2/HER2 are positively associated with
stemness for ESCA and STAD (Supplementary Figures 14-15), which aligns with
past knowledge - increased expression of HER2 leads to quicker growth and
poorer progression of gastric and esophageal cancers (Malaguti et al., 2015;
Lee et al., 2019).

#### EGFR exhibiting contradictory associations with gastrointestinal cancers

The EGFR gene/protein emerge jointly on top in the pan-GI mechanistic models,
with consistent decisive evidence (Figure 5C, main manuscript). Among gastric
cancer patients, 2–35% are reported to have EGFR protein overexpression and/or
gene amplification (Adashek et al., 2020). However, the utility of EGFR as a
therapeutic target has been questionable at best so far, with no general
consensus on its prognostic value in gastric cancer (Arienti et al., 2019).
While some studies have indicated that high EGFR gene amplification is
associated with poor outcome (Kandel et al., 2014), others have suggested the
opposite (Aydin et al., 2014). In line with this, cBVS found contradictory
associations - EGFR gene expression is negatively associated with both
stemness and survival for ESCA but positively associated with both for STAD
(Supplementary Figures 14-15, 32-33).

###  3.4 Investigating Biomarker Associations with Stemness Indices

One concern related to the stemness outcome analysis arises due to the fact
that the mRNA-based stemness indices (mRNAsi) are constructed by Malta et al.,
(2018) using a one-class logistic regression (OCLR) model that uses a large
pool of genes as the covariates. Potentially, these genes could overlap with
the panel of genes used in our integrative proteogenomic analyses and
resultingly, some of the associations detected by our stemness cBVS models may
be spurious and caused only by the bias towards these genes in the stemness
index training models. To reassure that this is not the case, we investigate
the overlap of the gene panel in the OCLR training model with that of our cBVS
models, and assess the weights of the common genes in the mRNAsi formula.

All the $160+$ genes that are included in our cBVS models are included in the
training phase of the OCLR model for stemness, which uses a total of $12,945$
genes. We first ordered the $12,945$ genes in decreasing order of the absolute
value of their weights in mRNAsi construction. None of the cBVS genes feature
in the top 50 of this list, only two of them feature in the top 100, with
three more each in the top 500 and the top 1000. For the overlapping genes,
the largest absolute value of the weights is $7.1\times 10^{-3}$ and the
median weight is $3.7\times 10^{-5}$. This indicates that there is not a
considerable effect of these genes in the mRNAsi. Further, we assessed the
correlations of mRNAsi (logit-transformed, as used in cBVS) for BRCA (the
largest sample size among the cancers of our interest) with the corresponding
gene expressions from the TCGA data. In summary, the largest absolute value
among these correlations is $<0.65$ and the median correlation is $0.03$.
Again, as before, this suggests that the mRNAsi is not spuriously correlated
with the expressions of the genes of our interest, and that the associations
detected by the stemness cBVS models are reliable.

##  4 Interactive R Shiny Dashboard

We provide an interactive R shiny dashboard summarizing all the results from
our pan-cancer integrative proteogenomic analyses. The main page of the
dashboard provides a brief overview of the fiBAG framework, as seen in
Supplementary Figure 34. Broadly, the shiny dashboard contains pages providing
interactive visualizations on the sample sizes, mechanistic model results, and
stemness and survival outcome model results for the different cancers. The
sample size summary page enables a user to choose a specific cancer, alongside
an outcome data type, in order to provide summaries in the form of Venn
diagrams presenting shared sample sizes (Supplementary Figure 35). The
mechanistic model summary page allows the user to explore the results for each
of the four cancer groups via heatmaps, upset plots, or wordclouds based on
the lBFs (Supplementary Figure 36). The stemness and survival outcome model
results are presented in two separate tabs, which work in similar fashion. For
example, the stemness outcome models are summarized for each cancer using
volcano-type plots, upset plots, and histograms of the evidence quantities
(Supplementary Figure 37). In its online de-anonymized version, the shiny
dashboard will have an additional page providing the user with the options to
(1) download the software codes used for the analyses and (2) visit the GitHub
page of the project. Due to the anonymity requirements of this submission, we
remove the said tab from the shiny dashboard and separately submit the codes
and the processed datasets in a compressed file.

##  5 Supplementary Figures

Supplementary Figure 1: Plate diagram summarizing the dependence structure
across the variables, parameters, and hyperparameters in fiBAG. Panels (A) and
(B) respectively summarize the mechanistic and outcome models.
$\mathbf{C},\mathbf{M},\mathbf{B},\mathbf{G}$, and $\mathbf{P}$ each have $n$
rows and are matched across samples, along with the outcome vector
$\mathbf{Y}$. $\mathbf{B},\mathbf{G}$, and $\mathbf{P}$ respectively have
$q_{b}$, $q_{g}$ and $q_{p}$ columns. $\mathbf{C}$ and $\mathbf{M}$
respectively represent copy number and methylation data. All the parameters
and data structures are described in full detail in Section 2 in the main
manuscript. Supplementary Figure 2: Results from Simulation comparing Gaussian
processes and linear models. The x-axis is in increasing order of non-
linearity (in terms of the mean function of the outcome based on the
covariates) in the generating model. The y-axis presents box plots across 100
replicates of the lBFs from Gaussian processes and linear models.
Supplementary Figure 3: Results from Simulation comparing MCMC and EMVS
implementations of cBVS. The x-axis presents the sample sizes for generating
the replicates, with 200 covariates in each case. The y-axis presents box
plots across 100 replicates of the mean squared differences of the
coefficients estimated via the two implementations.

Supplementary Figure 4: Integrative analysis using the fiBAG pipeline. In each
mechanistic model presented in the middle panel, upstream covariates are
matched to each proteogenomic biomarker. The lBFs computed from each model are
used in the cBVS outcome models to perform calibrated variable selection, as
outlined in the right panel. The outcome models provide posterior inclusion
probabilities for each biomarker along with estimated coefficient for
association with the outcome of choice. The entire pipeline is independently
applied on each cancer. We build one survival and one stemness outcome model
for each cancer.

Supplementary Figure 5: Mechanistic model heatmaps for the pan-kidney cancers.
Each cancer column consists of three sub-columns, one each for the three
mechanistic models (driver gene, driver protein and cascading protein). The
lBF ranges are defined as: $<0.5$ (no evidence), $0.5-1$ (substantial), $1-2$
(strong), $>2$ (decisive). Only the gene-protein pairs which are at the
decisive level of significance across all four cancers in at least two out of
the three mechanistic model types are shown here. Supplementary Figure 6:
Mechanistic model upset plots for the pan-kidney cancers. Upset plots exhibit
the number of genes (panel A) or proteins (panels B-C) that are at the
decisive level of significance (lBF $>2$) for the (A) driver gene, (B) driver
protein, and (C) cascading protein mechanistic models respectively, stratified
by intersections across cancers. Supplementary Figure 7: Mechanistic model
word clouds for the pan-kidney cancers. Word clouds summarize pan-cancer
mechanistic model results for genes (panel A) or proteins (panels B-C) for the
(A) driver gene, (B) driver protein, and (C) cascading protein mechanistic
models. The size of the gene/protein names are proportional to $(\textrm{no.
of cancers where the gene/protein is at the decisive level of
significance})^{3}$. Here, decisive is defined as lBF $>2$. Supplementary
Figure 8: Mechanistic model heatmaps for the pan-squamous cancers. Each cancer
column consists of three sub-columns, one each for the three mechanistic
models (driver gene, driver protein and cascading protein). The lBF ranges are
defined as: $<0.5$ (no evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$
(decisive). Only the gene-protein pairs which are at the decisive level of
significance across all four cancers in at least two out of the three
mechanistic model types are shown here. Supplementary Figure 9: Mechanistic
model upset plots for the pan-squamous cancers. Upset plots exhibit the number
of genes (panel A) or proteins (panels B-C) that are at the decisive level of
significance (lBF $>2$) for the (A) driver gene, (B) driver protein, and (C)
cascading protein mechanistic models respectively, stratified by intersections
across cancers. Supplementary Figure 10: Mechanistic model word clouds for the
pan-kidney cancers. Word clouds summarize pan-cancer mechanistic model results
for genes (panel A) or proteins (panels B-C) for the (A) driver gene, (B)
driver protein, and (C) cascading protein mechanistic models. The size of the
gene/protein names are proportional to $(\textrm{no. of cancers where the
gene/protein is at the decisive level of significance})^{3}$. Here, decisive
is defined as lBF $>2$. Supplementary Figure 11: Plot summarizing outcome
model results based on stemness indices (SI) for TCGA cancer UCS. Proteins are
represented by triangles and genes by circles. The shapes are colored red if
the estimated $\hat{\beta}_{j}$ from fiBAG is negative, and green if positive.
The x-axis shows the $\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).
Supplementary Figure 12: Plot summarizing outcome model results based on
stemness indices (SI) for TCGA cancer KICH. Proteins are represented by
triangles and genes by circles. The shapes are colored red if the estimated
$\hat{\beta}_{j}$ from fiBAG is negative, and green if positive. The x-axis
shows the $\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).
Supplementary Figure 13: Plot summarizing outcome model results based on
stemness indices (SI) for TCGA cancer ESCA-squamous. Proteins are represented
by triangles and genes by circles. The shapes are colored red if the estimated
$\hat{\beta}_{j}$ from fiBAG is negative, and green if positive. The x-axis
shows the $\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).
Supplementary Figure 14: Plot summarizing outcome model results based on
stemness indices (SI) for TCGA cancer ESCA-adeno. Proteins are represented by
triangles and genes by circles. The shapes are colored red if the estimated
$\hat{\beta}_{j}$ from fiBAG is negative, and green if positive. The x-axis
shows the $\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).
Supplementary Figure 15: Plot summarizing outcome model results based on
stemness indices (SI) for TCGA cancer STAD. Proteins are represented by
triangles and genes by circles. The shapes are colored red if the estimated
$\hat{\beta}_{j}$ from fiBAG is negative, and green if positive. The x-axis
shows the $\hat{\beta}_{j}$s, and the y-axis shows the
$-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for multiple
comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value type
quantity at the 10% FDR level. Only the selected biomarkers are marked in non-
gray colors and labeled. The sizes of the points are in the increasing order
of evidence from the mechanistic models: lBF ranges are defined as: $<0.5$ (no
evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).
Supplementary Figure 16: Bar diagram summarizing lBFs for the selected
proteogenomic biomarkers in the stemness outcome model for TCGA cancer BRCA.
An FDR check to adjust for multiple comparisons is performed treating
$1-\hat{\omega}_{j}$ as a p-value type quantity at the 10% FDR level. Only the
selected biomarkers are presented. Supplementary Figure 17: Histogram of lBFs
for the non-selected proteogenomic biomarkers in the stemness outcome model
for TCGA cancer BRCA. An FDR check to adjust for multiple comparisons is
performed treating $1-\hat{\omega}_{j}$ as a p-value type quantity at the 10%
FDR level. Only the non-selected biomarkers are included. Supplementary Figure
18: Bar diagram summarizing lBFs for the selected proteogenomic biomarkers in
the stemness outcome model for TCGA cancer CORE. An FDR check to adjust for
multiple comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value
type quantity at the 10% FDR level. Only the selected biomarkers are
presented. Supplementary Figure 19: Histogram of lBFs for the non-selected
proteogenomic biomarkers in the stemness outcome model for TCGA cancer CORE.
An FDR check to adjust for multiple comparisons is performed treating
$1-\hat{\omega}_{j}$ as a p-value type quantity at the 10% FDR level. Only the
non-selected biomarkers are included. Supplementary Figure 20: Plot
summarizing outcome model results based on overall survival for TCGA cancer
BRCA. Proteins are represented by triangles and genes by circles. The shapes
are colored red if the estimated $\hat{\beta}_{j}$ from fiBAG is negative, and
green if positive. The x-axis shows the $\hat{\beta}_{j}$s, and the y-axis
shows the $-\log_{10}(1-\hat{\omega}_{j})$s. An FDR check to adjust for
multiple comparisons is performed treating $1-\hat{\omega}_{j}$ as a p-value
type quantity at the 10% FDR level. Only the selected biomarkers are marked in
non-gray colors and labeled. The sizes of the points are in the increasing
order of evidence from the mechanistic models: lBF ranges are defined as:
$<0.5$ (no evidence), $0.5-1$ (substantial), $1-2$ (strong), $>2$ (decisive).