Joint estimation of multiple related biological networks

Consider a space $\mathcal{G}$ of graphs on the vertex set $\mathcal{P}=\{1,\ldots,P\}$. To keep the presentation general, we do not specify the type of graph or restrictions on $\mathcal{G}$ at this stage (the special case of DBNs for time-course data is described below). As shown in Figure 1, each individual network $G^{j}\in\mathcal{G}$ is modeled with dependence on a latent network $G\in\mathcal{G}$ that in turn depends on a prior network $G^{0}\in\mathcal{G}$ (Section 2.2). In this way, estimates of the individual networks $G^{j}$ are regularized by shrinkage toward the common latent network $G$ that, in turn, may be constrained by an informative network prior. As in any graphical model, observations $\mathbf{Y}^{j}$ on individual $j$ are dependent upon a graph $G^{j}$ and parameters $\bolds{\theta}^{j}$. Here $Z^{j}$ denotes any ancillary information available on individual $j$. The model is specified by

and a suitably chosen graphical model likelihood $p(\mathbf{Y}^{j}|G^{j},\bolds{\theta}^{j},Z^{j})$. Equation (1) follows the “network prior” approach of Mukherjee and Speed (2008) that was proposed for biological applications where subjective prior structural information is available. The functionals $d^{j},d\dvtx\mathcal{G}\times\mathcal{G}\rightarrow\mathbb{R}$ and hyperparameters $\eta,\lambda^{j}$ must be specified (Section 2.2). This paper restricts attention to exchangeable models, in particular, we consider functionals $d^{j}$ that are independent of the index $j$. We refer to the above formulation as joint network inference (JNI).

The network prior [equation (1)] requires a penalty functional $d\dvtx\mathcal{G}\times\mathcal{G}\rightarrow\mathbb{R}$ and a prior network $G^{0}\in\mathcal{G}$, with the former capturing how close a candidate network $G\in\mathcal{G}$ is to the latter. We discuss choice of $G^{0}$ below. Given $G^{0}$, a simple choice of penalty function $d$ is the structural Hamming distance (SHD) given by $d(G,G^{0})=\|G-G^{0}\|$, where $\|M\|=\sum_{i,j}|m_{i,j}|$ is the $\ell_{1}$-norm of an adjacency matrix and the differential network $G-G^{0}$ is defined to have edges that occur in exactly one of the networks $G$, $G^{0}$ [see also Ibrahim and Chen (2000), Imoto et al. (2003)]. The hyperparameter $\eta$ controls the strength of the prior network $G^{0}$ [equation (1)]. Motivated by an application in cancer biology where prior structural information $G^{0}$ is available, we follow Penfold et al. (2012) by restricting attention to SHD priors, however, our statistical formulation is general (see below) and compatible with other penalty functionals. Alternatively, one could employ a beta-binomial prior as described in, for example, Dondelinger, Lèbre and Husmeier (2013), that allows for the hyperparameters of the binomial to be integrated out [see also Oyen and Lane (2013)]. Note that in the latter case it is not possible to integrate specific prior structural information, making beta-binomial priors unsuitable for the application that this paper considers.

Given a latent network $G$, individual networks $G^{j}$ are regularized in a similar way, as $d^{j}(G^{j},G)=\|G^{j}-G\|$. In their work on combining multiple data sources, Werhli and Husmeier (2008) allow the $\lambda^{j}$ to vary over individuals $j\in\mathcal{J}$. Likewise, Penfold et al. (2012) learn the $\lambda^{j}$ on a per-individual basis. However, in both studies, hyperparameter elicitation is nontrivial (see Section 3.3). In the present paper, we consider only the special case where $\lambda^{1}=\lambda^{2}=\cdots=\lambda^{J}:=\lambda$.

A graph $G\in\mathcal{G}$ can be characterized by (i) its adjacency matrix or (ii) its parent sets as $G=(\pi_{1},\ldots,\pi_{P})$, where $\pi_{p}\subseteq\mathcal{P}=\{1\cdots P\}$ are the parents of vertex $p$ in $G$. We write $\mathcal{G}_{p}$ for the set of possible parent sets for $p$, such that formally $\mathcal{G}=\mathcal{G}_{1}\times\cdots\times\mathcal{G}_{P}$. Although we focus on SHD priors, the inference procedures presented in this paper apply to the more general class of “modular” priors, that may be factored over $p\in\mathcal{P}$ and written in the form

for some functionals $d_{p},d_{p}^{j}\dvtx\mathcal{G}_{p}\times\mathcal{G}_{p}\rightarrow\mathbb{R}$. Here $\pi_{p}^{0}$ and $\pi_{p}^{j}$ are parent sets for variable $p$, corresponding to $G^{0}$ and $G^{j}$, respectively.

In general, inference for the JNI model [equations (1), (2)] will be computationally intensive, as demonstrated in Werhli and Husmeier (2008), Penfold et al. (2012). In Section 3 below we show that efficient, exact inference is nevertheless possible within the DBN class of graphical models.

A DBN is a graphical model based on a DAG on the vertex set $\mathcal{P}\times\mathcal{T}$, where $\mathcal{T}$ is a set of time indices [Figure 8(a); see Murphy (2002)]. This DAG with $PT$ vertices is known as the “unrolled” DAG. Here, following Hill et al. (2012) and others, we use DBNs that permit only edges forward in time and that are stationary in the sense that neither the network nor parameters change with time. For such DBNs, the network can be described by a directed graph $G$ with exactly $P$ vertices, with edges understood to go forward in time in the unrolled DAG [see Appendix B and Figure 8(b)]. Note that $G$ may have cycles. In what follows, all graphs (prior, latent and individual) describe the latter $P$-vertex representation.

Under a modular network prior, structural inference for DBNs can be carried out efficiently as described in Hill et al. (2012). In brief, the posterior $G^{j}|\mathbf{y}$ factorizes into a product of local posteriors $\pi_{p}^{j}|\mathbf{y}$, one factor for each target variable $p$. Background and assumptions for DBNs are summarized in Appendix B; for general background on DBNs we refer the interested reader to Murphy (2002) and for relevant details concerning the class of DBNs used here to Hill et al. (2012).

Write $\mathbf{y}(t)$ for the matrix of observed data at time $t$ for all individuals $j$ and variables $p$. In order to simplify notation, we define a data-dependent functional

that implicitly conditions upon observed history. Let $y_{p}^{j}(t)$ denote the observed value of variable $p$ in individual $j$ at time $t$. The above notation allows us to conveniently summarize the product

as $\mathfrak{P}(\mathbf{y}_{p}^{j}|\pi_{p}^{j})$. Thus, we have that, for DBNs, the full likelihood also satisfies

where $\mathbf{y}$ denotes the complete data (for all individuals, variables and times). In other words, the parent sets $\pi_{p}^{j}$ for $p\in\mathcal{P}$, $j\in\mathcal{J}$ are mutually orthogonal in the Fisher sense, so that inference for each may be performed separately.

We carry out exact inference in this setting using belief propagation [Pearl (1982)]. Belief propagation is an iterative procedure in which messages are passed between variables in such a way as to compute exact marginal distributions; in this respect belief propagation belongs to a more general class of iterative algorithms known as “sum-product” algorithms [Kschischang, Frey and Loeliger (2001)]. Our algorithm is summarized as follows (for simplicity we suppress dependence upon ancillary information $Z^{j}$):

[(3)]

We begin by marginalizing over parameters $\bolds{\theta}^{j}$ and caching the local scores $\mathfrak{P}(\mathbf{y}_{p}^{j}|\pi_{p}^{j})$ for all parent sets $\pi_{p}^{j}\in\mathcal{G}_{p}$, all variables $p\in\mathcal{P}$ and all individuals $j\in\mathcal{J}$; these could be obtained using any DBN likelihood. In this paper we exploited conjugate priors to obtain exact expressions for marginal likelihoods [equation (C), see Appendix C for details].

Following marginalization, the JNI graphical model collapses to the discrete Bayesian network shown in Figure 2, whose nodes are themselves graphs.

Posterior marginal distributions $p(\pi_{p}|\mathbf{y}_{p},\pi_{p}^{0})$ and $p(\pi_{p}^{j}|\mathbf{y}_{p},\pi_{p}^{0})$ are then computed using belief propagation on this discrete Bayesian network. Pseudocode for this step is provided in Algorithm 1 in Appendix D.

Let $\mathbf{A}_{\mathrm{JNI}}$ denote the $P\times P$ matrix of marginal posterior inclusion probabilities for edges in the latent network $G$, that is, $(\mathbf{A}_{\mathrm{JNI}})_{ip}:=p(i\in\pi_{p}|\mathbf{y},G^{0})$. These quantities are analogous to posterior inclusion probabilities in Bayesian variable selection and are computed, using Bayesian model averaging, as

where $\mathbf{1}_{A}$ is the indicator of the event $A$ and similarly for individual networks $(\mathbf{A}_{\mathrm{JNI}}^{j})_{ip}:=p(i\in\pi_{p}^{j}|\mathbf{y},G^{0})$.

Following Hill et al. (2012), we reduced the space of parent sets $\mathcal{G}_{p}$ using an in-degree sparsity restriction of the form $|\pi_{p}^{j}|\leq c$ for all $\pi_{p}^{j}\in\mathcal{G}_{p}$, $p\in\mathcal{P}$, $j\in\mathcal{J}$. Thus, the cardinality of the space of parent sets $|\mathcal{G}_{p}|=\mathcal{O}(P^{c})$ is polynomial in $P$, where it was previously super-exponential. As in variable selection, the bound $c$ should be chosen large enough that $\mathcal{G}_{p}$ includes the true data-generating model with high probability.

Caching of selected probabilities is used to avoid redundant recalculation. Pseudocode is provided in Algorithm 1 in Appendix D, consisting of three phases of computation. Storage costs are dominated by phases I and II, each requiring the caching of $\mathcal{O}(JP^{1+c})$ terms. Phase II dominates computational effort, with total (serial) algorithmic complexity $\mathcal{O}(J^{2}P^{1+2c})$. However, within-phase computation is “embarrassingly parallel” in the sense that all calculations are independent (indicated by square parentheses notation in the pseudocode). In practice, we have found that problems of size $P\leq 20$, $J\leq 20$, $c\leq 3$ can be solved within minutes using serial computation on a standard laptop computer. We provide serial and parallel MATLAB R2014a implementations in Supplement B [Oates et al. (2014b)].

Elicitation of hyperparameters for network priors is an important and nontrivial issue. Here we specify the hyperparameters $\lambda,\eta$ in a subjective manner. We do so due to reported difficulties in estimation of hyperparameters for related models [Werhli and Husmeier (2008), Dondelinger, Lèbre and Husmeier (2013), Penfold et al. (2012)]. We present three criteria below that, for the special case of SHD, are simple to implement and can be used for expert elicitation. These heuristics seek to relate the hyperparameters to more directly interpretable measures of the similarity and difference that they induce between prior, latent and individual networks: (i) First, we note the following formula for the probability of maintaining the status (present/absent) of a candidate parent $i\in\mathcal{P}$ between the latent network $G$ and an individual network $G^{j}$:

This probability provides an interpretable way to consider the influence of $\lambda$. For example, a prior confidence of $h_{\lambda}\approx 0.73$ that a given edge status in $G$ is preserved in a particular individual $G^{j}$ translates into an odds ratio $h_{\lambda}/(1-h_{\lambda})\approx 2.7$ and a hyperparameter $\lambda\approx 1$ (see SFigure 1 in the supplementary material [Oates et al. (2014a)]). An analogous equation relates $\eta$ and $h_{\eta}:=p(i\notin\pi_{p}\Delta\pi_{p}^{0})$, allowing prior strength to be set in terms of the probability that an edge status in the prior network $G^{0}$ is maintained in the latent network $G$. (ii) A second, related approach is to consider the expected total SHD between an individual network $G^{j}$ and the latent network $G$:

This can be interpreted as the average number of edge changes needed to obtain $G^{j}$ from $G$. An analogous equation holds for $\eta$ and $h_{\eta}$. (iii) Third, in certain applications, the latent network $G$ may not have a direct scientific interpretation, in which case the criteria presented above may be unintuitive. Then, hyperparameters can be elicited by consideration of (a) similarity between individual networks $G^{j},G^{k}$ and (b) concordance of individual networks $G^{j}$ with the prior network $G^{0}$ (see Supplement A [Oates et al. (2014a)] for further discussion).

Below we consider the extent to which information can be shared between individuals within JNI, providing an upper bound that is attained as the number of individuals $J$ grows large. To formalize the contribution to inference from information sharing, we consider the case in which no data is available on a specific individual (without loss of generality, individual $j=1$) and analytically quantify the extent to which JNI can estimate the true network $\overline{G^{1}}$ by “borrowing strength” from the data $\mathbf{Y}^{2},\ldots,\mathbf{Y}^{J}$ that represent observations on the remaining individuals. (Over-lines will be used to signify the “true” data-generating networks.) As a baseline, write $\mathbf{A}_{0}^{j}=p(i\in\pi_{p}^{j}|\mathbf{Y}^{j})$ for the (naive) estimator that prohibits the sharing of information between individuals. For simplicity we restrict attention to the case where no network prior is used ($\eta=0$), the data-generating hyperparameter $\lambda$ is known and in-degree restrictions are not in place ($c=P$). Then, with neither data nor prior information available on individual $1$, it trivially follows that

where the expectation is taken over all possible data-generating networks and corresponding data.

From standard, independent network inference we know that consistent estimation requires unbiasedness of the likelihood function $p(\mathbf{Y}^{j}|G^{j})$, in the sense that $\mathbb{E}_{\mathbf{Y}^{j}|\overline{G^{j}}}p(\mathbf{Y}^{j}|G^{j})$ is maximized by $G^{j}=\overline{G^{j}}$. We therefore begin by constructing the analogous regularity condition for joint estimation: Write $\mathbf{R}$ for the matrix that encodes the prior metric on $\mathcal{G}$ as $(\mathbf{R})_{G,G^{\prime}}=\exp(-\lambda\|G-G^{\prime}\|)/C(\lambda)$, where $C(\lambda)=\sum_{G\in\mathcal{G}}\exp(-\lambda\|G\|)$. Write $\mathbf{S}$ for the matrix of expected Bayesian scores $(\mathbf{S})_{G^{j},\overline{G^{j}}}=\mathbb{E}_{\mathbf{Y}^{j}|\overline{G^{j}}}p(\mathbf{Y}^{j}|G^{j})$.

[(Joint regularity)] For each column of the matrix $\mathbf{M}=(\mathbf{R}\mathbf{S}\mathbf{R})_{G,\overline{G}}$, the nondiagonal entries are strictly smaller than the diagonal entry, that is, $M_{G,\overline{G}}<M_{\overline{G},\overline{G}}$ for all $G\neq\overline{G}$.

To gain intuition for the joint regularity assumption, consider the special case where $\lambda\rightarrow\infty$; here $\mathbf{R}=\mathbf{I}$ and we only require that the expected local Bayesian score $(\mathbf{S})_{G^{j},\overline{G^{j}}}$ is maximized by $G^{j}=\overline{G}$, that is, we recover the unbiasedness condition from standard network inference.

Under the joint regularity assumption, there exists $0<\varepsilon<1$ such that

where $f(J)\leq 2P^{2}\varepsilon^{J-1}\rightarrow 0$ as $J\rightarrow\infty$.

See Appendix A.

Comparing equation (11) to (10), we see that information sharing offers gains in estimation, agreeing with intuition, with larger gains when the true networks are almost homogeneous ($\lambda$ large). Moreover, the statistical power of JNI to estimate $\overline{G^{1}}$ converges to its maximum exponentially quickly as $J\rightarrow\infty$.

The proposed methodology addresses three questions, some or all of which may be of scientific interest depending on the application: (i) estimation of the latent network $G$, (ii) estimation of individual networks $G^{1},\ldots,G^{J}$, and (iii) estimation of differences between individual networks [“differential networks”; Ideker and Krogan (2012)]. We quantify performance for each task using the area under the receiver operating characteristic (ROC) curve (AUR). This metric, equivalent to the probability that a randomly chosen true edge is preferred by the inference scheme to a randomly chosen false edge, summarizes, across a range of thresholds, the ability to select edges in the data-generating network. AUR may be computed relative to the true latent network $G$ or relative to the true individual networks $G^{j}$, quantifying performance on tasks (i) and (ii), respectively. Both sets of results are presented below, in the latter case averaging AUR over all individual networks. For (iii), in order to assess ability to estimate differential networks, we computed AUR scores based on the statistics $F_{ip}^{j}=|p(i\in\pi_{p}^{j}|\mathbf{y},G^{0},Z^{j})-p(i\in\pi_{p}|\mathbf{y},G^{0},Z^{1},\ldots,Z^{J})|$ that should be close to one if $i\in\pi_{p}^{j}\Delta\pi_{p}$, otherwise $F_{ip}^{j}$ should be close to zero.

It is easy to show that inference for the latent network, under only the prior (i.e., $\hat{G}=G^{0}$), attains mean AUR equal to $h_{\eta}$. Similarly, prior inference for the individual networks (i.e., $\hat{G}^{j}=G^{0}$) attains mean AUR equal to $1-h_{\eta}-h_{\lambda}+2h_{\eta}h_{\lambda}$. This provides a baseline for the proposed methodology at tasks (i) and (ii) and allows performance to be decomposed into AUR due to prior knowledge and AUR contributed through inference.

Using a systematic variation of data-generating parameters, we defined 15 distinct data-generating regimes described below. For all 15 regimes we considered 50 independent data sets; standard errors accompany average AUR scores. Results presented below use a computationally favorable in-degree restriction $c=3$. In order to check robustness to $c$, a subset of experiments were repeated using $c=4$, with close agreement observed (SFigure 4 in the supplementary material [Oates et al. (2014a)]).

We consider experimental data derived from human breast cancer cell lines, focusing on protein signaling networks within which many (wild type) causal relationships are well understood from extensive biochemistry (Figure 3). The investigation presented below serves three purposes: First, it allows investigation of the applicability of the proposed joint approaches to experimental data. Second, it allows investigation of the use of ancillary information, in the form of mutational status and histological information. Finally, the results and approach are relevant to the topical question of exploring signaling heterogeneity across cancer cell lines.

Data were obtained using reverse-phase protein arrays [Hennessy et al. (2010)] from $J=6$ breast cancer cell lines (AU 565, HCC 1569, MCF 7, MDA MB 231, SKBR3 and SUM 190PT; experimental protocol is described in brief in Supplement A [Oates et al. (2014a)]). Data comprised observations for the $P=17$ proteins shown in Figure 3 (see also STable 1 in the supplementary material [Oates et al. (2014a)]; we note that these data form part of a larger study including further cell lines and proteins). Specifically, $\mathbf{y}$ contains the logarithms of the measured concentrations. Data were acquired under treatment with an EGFR/HER2 inhibitor Lapatinib (“EGFRi”), an Akt inhibitor (“Akti”), EGFRi and Akti in combination, and without inhibition (“DMSO”) at 0.5, 1, 2, 4, 8 and 24 hours following Serum stimulation, giving a total of $n_{j}=24$ observations of each variable in each individual cell line.