 Research
 Open Access
 Published:
Bayesian graphical models for computational network biology
BMC Bioinformatics volume 19, Article number: 63 (2018)
Abstract
Background
Computational network biology is an emerging interdisciplinary research area. Among many other network approaches, probabilistic graphical models provide a comprehensive probabilistic characterization of interaction patterns between molecules and the associated uncertainties.
Results
In this article, we first review graphical models, including directed, undirected, and reciprocal graphs (RG), with an emphasis on the RG models that are curiously underutilized in biostatistics and bioinformatics literature. RG’s strictly contain chain graphs as a special case and are suitable to model reciprocal causality such as feedback mechanism in molecular networks. We then extend the RG approach to modeling molecular networks by integrating DNA, RNA and proteinlevel data. We apply the extended RG method to The Cancer Genome Atlas multiplatform ovarian cancer data and reveal several interesting findings.
Conclusions
This study aims to review the basics of different probabilistic graphical models as well as recent development in RG approaches for network modeling. The extension presented in this paper provides a principled and efficient way of integrating DNA copy number, DNA methylation, mRNA gene expression and protein expression.
Background
This article starts with a comprehensive review of graphical models and recent development in constructing biological networks using reciprocal graphs (RG’s). In order to integrate multiomics data including proteomics, we extend the work in [1] utilizing fundamental biological knowledge and the factorization of the joint distribution.
Computational network biology (CNB) is an emerging research field that encompasses theory and applications of network models to systematically study different molecules (DNA, RNA, proteins, metabolites and small molecules) and their complex interactions in living cells. CNB provides new insights in diverse research areas including system biology, molecular biology, genetics, pharmacology and precision medicine. The development of highthroughput profiling technologies allows for the interrogation of the status of many molecules in a cell at the same time, which makes it possible to jointly model these cell components with network modeling approaches. CNB is an interdisciplinary research area and has received contributions from researchers with various background. In this paper, we focus on network modeling approaches from a statistical perspective, namely, probabilistic graphical models (PGM) which allow a complete probabilistic description of interaction patterns between molecules and the associated uncertainties when estimating such interactions from noisy genomic and/or proteomic data.
PGM are probabilistic models for multivariate random variables whose conditional independence (also known as Markov property) structure is characterized by an underlying graph. PGM provides a concise, complete and explicit representation of joint distribution and allows for convenient Gibbs factorization of the density (if it exists) and hence local computations. The Markov property can be directly read off the graph through the notion of graph separation. PGM is closely related to causal inference, path analysis and expert system and finds a wide range of applications in political sciences, economics, genetics, biology, physics and psychology.
The most commonly used graphs are undirected graphs (UG) and directed acyclic graphs (DAG). UG represents symmetric interactions between variables whereas DAG allows for asymmetric ones that can be potentially interpreted as causeeffect relationships. Both UG and DAG are a special case of a more general graphical model known as chain graphs (CG) [2, 3] which allow for both symmetric and asymmetric conditional independence relationships. Despite of the flexibility of CG, the directed relationships are still constrained to be (block) recursive, that is, reciprocal causality is prohibited. However, cyclic causal relationship are fundamental and ubiquitous in science, including, for example, the law of supply and demand in economics or feedback mechanism in gene regulatory networks. Spirtes [4] and Koster [5], in their respective seminal work, independently proposed coherent nonrecursive graphical models including directed cycles of which the joint distribution was previously thought illdefined ([6], p. 72). Spirtes [4] developed directed cyclic graphical (DCG) models which include DAG’s as a special case, and [5] proposed an even more general class of graphical models, reciprocal graphical (RG) models, which include CG (hence UG and DAG) and DCG as special cases. The relationships between the aforementioned graphs are depicted as a Venn’s diagram in Fig. 1.
In this paper, we focus on RG’s due to their statistical generality, the capability of modeling feedback loops and their still underappreciated popularity in the biostatistics and bioinformatics literature. The rest of this paper is structured as follows. We first review the basics of UG, DAG and RG in this section. In the next section, we provide a comprehensive review of novel Bayesian approaches in modeling molecular networks using RG. Then, we extend the RG approach in [1] to modeling The Cancer Genome Atlas (TCGA) ovarian cancer data that are measured on three different levels: DNA, RNA and proteins. The extension is nontrivial: we factorize an RG into two separate RG’s (one for DNA and RNA and the other for RNA and protein) and make coherent joint inference by exploiting two known biochemical processes: DNA transcribes to RNA and RNA in turn translates to protein, the latter also allowing us to estimate the direction of the relationship between RNA and protein which informs the directionality of the regulatory relationships between proteins. To the best of our knowledge, we are among the first to reconstruct multiomics functional networks that consider mass spectrometry proteomics data generated by the CPTAC.
Overview of directed, undirected and reciprocal graphical models. Graphical models are a class of statistical models for a set of random variables Z=(Z_{1},…,Z_{ p })^{T} that provide a graphic representation of conditional independence relationships among the variables. Using graphical models allow practitioners to simplify inferences, obtain parsimonious solutions via variable selection, enable local computations, and define a framework for causal inference.
In a graphical model, variables are represented by a set of nodes and their associated interactions are represented by edges. A missing edge usually represents the conditional independence between the corresponding pair of random variables conditional on a certain set of other variables. The conditioning set depends on the chosen type of graph. Here, we provide a minimal set of terminologies required to describe the Markov properties in later sections.
A graph \(\mathcal {G} = (V, E)\) is defined by a set of nodes V={1,…,p} and a set of directed and undirected edges E=E^{d}∪E^{u}. For a pair of nodes i,j∈V, we denote an undirected edge by i−j or {i,j}∈E^{u} and a directed edge by i→j or (i,j)∈E^{d}. We write i↦j if there exists a path from i to j. A path is an ordered sequence (i_{0},…,i_{ K }) of distinct nodes except possibly i_{0}=i_{ K } such that {i_{k−1},i_{ k }}∈E or (i_{k−1},i_{ k })∈E for k=1,…,K. A cycle is then a path with i_{0}=i_{ K }. A path component is a set of nodes that are all connected by an undirected path. The boundary of a node i is bd(i)={j∣{j,i}∈E or (j,i)∈E} and the boundary of a subset A⊆V is \(\text {bd}(A)=\bigcup _{i\in A}\text {bd}(i)\backslash A\). A node i is said to be a parent of node j, if (i,j)∈E. The set of parents of j is denoted by pa(j). The ancestors of node i are an(i)={i}∪{jj↦i}. The descendants of node i, denoted by de(i), are the nodes i such that there is a path from i to j. The nondescendants of i are nd(i)=V∖(de(i)∪{i}). For A⊆V, the induced subgraph by A is \(\mathcal {G}_{A}=(A,E_{A})\) where E_{ A }=E∩(A×A). A subset A⊆V is anterior if bd(A)=∅ and the minimal anterior set of A is denoted by an(A). Note that we use an() to denote both ancestors and minimal anterior set; it should not cause any confusion because in fact an({i})=an(i). Notice that (minimal) anterior set [5] is equivalent to (minimal) ancestral set in [7]. A graph is complete if all nodes are joined by an edge. A complete subgraph that is maximal with respect to the subset relation is called clique. A useful procedure in defining Markov properties of graphs with directed edges is moralization. The moralization is performed in two steps. Make the boundary of each path component complete by adding undirected edges and then replace all directed edges by undirected ones. We call the resulting undirected graph moral graph, denoted by \(\mathcal {G}^{m}=(V,E^{m})\).
We say that graph \(\mathcal {G}\) is a undirected graph (UG) if E=E^{u} contains only undirected edges. And graph \(\mathcal {G}\) is said to be a directed acyclic graph (DAG) if E=E^{d} contains only directed edges and no cycles. A reciprocal graph (RG) is a graph \(\mathcal {G}\) that can have both directed and undirected edges and cycles, with the only restriction being no directed edges between nodes in the same path component [5]. In the following, we provide a concise overview of each type of graph. More comprehensive and technical details of the contents can be found in [7] and [6].
Undirected graphical models
Let P denote the joint distribution of Z=Z_{ V }=(Z_{1},…,Z_{ p })^{T}. We say P has a Gibbs factorization (or simply factorization) with respect to a UG \(\mathcal {G}\) if its density f(Z) can be written as
where \(\mathcal {C}\) is the collection of cliques of \(\mathcal {G}\) and ψ_{ C } is an arbitrary nonnegative function on the domain of Z_{ C } for \(C\in \mathcal {C}\).
Not surprisingly, the factorization property of P is closely related to its Markov properties; in fact, the former implies the latter. The Markov property of a UG relies on a notion of graph separation. For three disjoint subsets A,B,C∈V of nodes, A is said to be separated from B by C if every path between A and B includes a node in C. The global Markov property of a UG is given in the following definition.
Definition 1
(Global Markov property) Given a UG \(\mathcal {G}\), the joint distribution P of Z_{ V } is global Markov with respect to \(\mathcal {G}\) if for any three disjoint subsets A, B and C of V, Z_{ A } is independent of Z_{ B } given Z_{ C }, written as Z_{ A } ⊥ ⊥Z_{ B }∣Z_{ C } whenever C separates A and B.
The global Markov property always implies the local Markov property.
Definition 2
(Local Markov property) Given a UG \(\mathcal {G}\), the joint distribution P of Z_{ V } is local Markov with respect to \(\mathcal {G}\) if for any node i∈V, Z_{ i } ⊥ ⊥Z_{V∖{bd(i),i}}∣Z_{bd(i)}.
The local Markov property in turn always implies the pairwise Markov property.
Definition 3
(Pairwise Markov property) Given a UG \(\mathcal {G}\), the joint distribution P of Z_{ V } is pairwise Markov with respect to \(\mathcal {G}\) if for any pair of nodes i,j∈V, Z_{ i } ⊥ ⊥Z_{ j }∣Z_{V∖{i,j}}.
In summary, factorization in (1) ⇒ global Markov property ⇒ local Markov property ⇒ pairwise Markov property. A natural question to ask is: is the relation also true backwards? The answer is yes when P has a positive and continuous density, which is summarized in the following theorem, often referred to as HammersleyClifford theorem [8].
Theorem 1
(HammersleyClifford) Suppose a probability distribution P has a positive and continuous density f. Then P satisfies the pairwise Markov property with respect to a UG \(\mathcal {G}\) if and only if f factorizes according to \(\mathcal {G}\).
When P=N(0,Λ^{−1}) is multivariate Gaussian with precision matrix Λ, undirected graphical models are defined by zero constraints on Λ. There is a onetoone correspondence between the graph structure \(\mathcal {G}\) and the zero patterns in Λ. Specifically, the precision matrix Λ is constrained to the cone of positive definite matrices with offdiagonal entry Λ_{ ij }=0 if and only if there is a missing edge in \(\mathcal {G}\) between nodes i and j. Therefore, in the Gaussian case, learning the graph structure is equivalent to recovering the zero patterns in Λ. Interested readers are referred to recent Bayesian approaches [9–11] and nonBayesian approaches [12, 13].
Directed graphical models
We say the distribution P factorizes with respect to a DAG \(\mathcal {G}\) if its density f(Z) can be written as
Similar to UG’s, the Markov property of a DAG can be also defined by the concept of graph separation. However, unlike UG’s, the separation is more complicated in DAG’s. There are two different but equivalent approaches in describing separation in a DAG. Pearl [14] introduced the notion of dseparation characterized by several conditions and then defined the global Markov property of a DAG in the same way as that of an UG except that separation is replaced by dseparation. Another approach [3] is based on the moralization and graph separation in UG’s. We elaborate more on the latter because it will be used again in defining Markov properties of RG’s. Similar to the previous section, we introduce global, local and pairwise Markov properties of DAG and then state the relationships between them. Recall that \(\mathcal {G}^{m}\) denotes the moralized graph.
Definition 4
(Directed global Markov property) Given a DAG \(\mathcal {G}\), the joint distribution P of Z_{ V } is global Markov with respect to \(\mathcal {G}\) if for any three disjoint subsets A, B and C of V, Z_{ A } is independent of Z_{ B } given Z_{ C }, written as Z_{ A } ⊥ ⊥Z_{ B }∣Z_{ C } whenever C separates A and B in \(\mathcal {G}_{an(A\cup B\cup C)}^{m}\).
Definition 5
(Directed local Markov property) Given a DAG \(\mathcal {G}\), the joint distribution P of Z_{ V } is local Markov with respect to \(\mathcal {G}\) if for any node i∈V, Z_{ i } ⊥ ⊥Z_{nd(i)}∣Z_{pa(i)}.
Definition 6
(Directed pairwise Markov property) Given a DAG \(\mathcal {G}\), the joint distribution P of Z_{ V } is pairwise Markov with respect to \(\mathcal {G}\) if for any pair of nodes i,j∈V, Z_{ i } ⊥ ⊥Z_{ j }∣Z_{nd(i)∖{j}}.
Similar to UG’s, it is easy to show that factorization in (2) ⇒ directed global Markov property ⇒ directed local Markov property ⇒ directed pairwise Markov property. In fact, directed global and local Markov properties are equivalent [15]. Moreover, a theorem similar to HammersleyClifford’s theorem holds for DAG’s as well.
Theorem 2
([7]) Suppose a probability distribution P has a continuous density f. Then P satisfies the local Markov property with respect to a DAG \(\mathcal {G}\) if and only if f factorizes according to \(\mathcal {G}\). Furthermore, if the density f is positive, local and pairwise Markov properties are equivalent.
In a Gaussian DAG, each conditional distribution in (2) is a linear model
where N(·∣μ,σ^{2}) is a Gaussian density with mean μ and variance σ^{2}. Essentially, a DAG is decomposed as a system of recursive and independent linear regressions for which the graph structure inference can be performed using standard model selection technique.
However, there is an important caveat in using DAG’s for statistical inference. That is the the notion of Markov equivalence – different DAG’s may induce the same set of conditional independence relationships. For example, the graphs i←j→k and i←j←k induce the same conditional independence relationship, i ⊥ ⊥k∣j. Both are equivalent to the UG i−j−k. In fact, any perfect DAG (i.e. a DAG without configuration i→j←k, that is, all parents are married) is Markov equivalent to a decomposable UG (i.e. UG in which all cycles of length four or more have a chord, which is an edge that is not part of the cycle but connects two nodes of the cycle). Conversely, given an undirected decomposable graph, there exists a Markovequivalent perfect DAG; see [7] or [16]. One approach to avoid repeatedly analyzing Markov equivalent DAG’s is adopting a prior ordering of the nodes. The ordering may be obtained from subjective prior knowledge such as a reference network [17] or more objectively from an ordering learning algorithm such as orderMCMC [18]. For other recent DAG approaches, see [19–21] or [22].
Reciprocal graphical models, simultaneous equation models and path diagrams
We introduce the global Markov property for RG’s in the following definition and omit the local and pairwise Markov properties since they are not as practically useful as in UG and DAG.
Definition 7
(Reciprocal global Markov property) Given a RG \(\mathcal {G}\), the joint distribution P of Z_{ V } is global Markov with respect to \(\mathcal {G}\) if for any three disjoint subsets A, B and C of V, Z_{ A } is independent of Z_{ B } given Z_{ C }, written as Z_{ A } ⊥ ⊥Z_{ B }∣Z_{ C } whenever C separates A and B in \(\mathcal {G}_{an(A\cup B\cup C)}^{m}\).
Similarly to UG and DAG, there is also a HammersleyCliffordtype theorem for RG’s ([5], Theorem 3.4).
However, while it is a beautiful theoretical result on the equivalence between the global Markov property and the Gibbs factorization in RG, the latter is difficult to work with in practice due to its complex form. Fortunately, in the Gaussian case, there is an easy connection between RG and simultaneous equation models (SEMs). To describe this connection, let Z=(Y,X) be divided into two sets of variables Y and X and consider an SEM,
where we assume the following five conditions:

C1.
(Y,X)∼N(0,Σ);

C2.
E ⊥ ⊥X;

C3.
I−A is invertible;

C4.
Cov(E) is diagonal;

C5.
Ψ=Cov(X) is block diagonal with full diagonal blocks.
To link an SEM to an RG, we draw a path diagram \(\mathcal {G}=(V,E)\) of SEM by the following rules:

(i)
define nodes V={1,…,p,p+1,…,p+q} which represent (Y,X)=(Y_{1},…,Y_{ p },X_{1},…,X_{ q });

(ii)
draw directed edges E^{d}={(j,i)∣a_{ ij }≠0 or b_{i,j−p}≠0}; and

(iii)
draw undirected edges E^{u}{{i,j}∣Ψ_{i−p,j−p}≠0}.
In words, (i) we introduce a node for each variable in (Y,X) with nodes j=1,…,p corresponding to Y_{ j } and p+j corresponding to X_{ j } for j=1,…,q; (ii) nodes i=1,…,p (i.e. Y_{ i } nodes) become targets of directed edges from node j if the corresponding a_{ ij }≠0 or b_{i,j−p}≠0; (iii) we introduce undirected edges between X_{ i } and X_{ j } (i.e. nodes i+p and j+p) if Ψ_{i,j}≠0. Figure 2a shows an example of an RG with p=2 and q=3,
with ∗ indicating nonzero elements. With conditions C1C5, the Markov properties of an SEM can be read off its path diagram.
Theorem 3
([5]) Let Y=AY+BX+E be an SEM satisfying conditions C1C5. The joint distribution P of (Y,X) is global Markov with respect to the path diagram \(\mathcal {G}\).
For example, the SEM with configuration in (4) implies X_{1} ⊥ ⊥X_{2},X_{3} which is evident by looking at the moral subgraph \(\mathcal {G}_{an(X_{1},X_{2},X_{3})}^{m}=\mathcal {G}_{\{X_{1},X_{2},X_{3}\}}^{m}\) in Fig. 2b. Clearly, the empty set ∅ separates X_{1} and X_{2},X_{3}. We can also find X_{1} ⊥ ⊥X_{2},X_{3}∣Y_{1},Y_{2} by moralizing the entire path diagram (Fig. 2c) since set {Y_{1},Y_{2}} separates X_{1} and X_{2},X_{3} in \(\mathcal {G}^{m}\). However, we remark that the class of SEMs satisfying conditions C1C5 is a subclass of RG. For instance, i→j−k←l is a CG (hence an RG) but does not correspond to an SEM. Similarly to DAG’s, DCG’s (or more generally RG’s) can also be divided into Markov equivalence classes; see [23–28] for characterization and learning of Markov equivalence classes of DCG’s.
Methods
In this section, we will review recent RG approaches in modeling molecular networks including [29, 30] and [1].
Modeling signaling pathways through latent variables
Telesca et al. [29] considered gene expressions y_{ ij } for p genes measured over n samples, i=1,…,n and j=1,…,p. They follow [31] and introduce trinary latent variables e_{ ij }∈{−1,0,1} defining three possible categories of each expression y_{ ij },
The e_{ ij } can be interpreted as an underlying biologic signal of under, over or normal expression. The motivation of the proposed inference is that only e_{ ij } should be modeled without overinterpreting the additional noise in y_{ ij } as biologically meaningful. Given e_{ ij }, a Gaussianuniform mixture model is assumed for gene expression y_{ ij },
where \(f_{1,j}=U(\kappa _{j}^{},0),f_{1,j}=U\left (0,\kappa _{j}^{+}\right),f_{0,j}=N\left (0,\sigma _{j}^{2}\right)\) and \(\tilde {y}_{ij}=y_{ij}(\alpha _{i}+\mu _{j})\) is the normalized relative abundance, corrected for a samplespecific effect α_{ i } and a genespecific effect μ_{ j }. Conditional conjugate priors are assigned to \(\alpha _{i},\mu _{j},\kappa _{j}^{},\kappa _{j}^{+}\): \(\alpha _{i}\sim N\left (0,\tau _{\alpha }^{2}\right), \mu _{j}\sim N\left (m_{\mu },\tau _{\mu }^{2}\right),\kappa _{j}^{}\sim IG\left (a_{\kappa }^{},b_{\kappa }^{}\right),\kappa _{j}^{+}\sim IG\left (a_{\kappa }^{+},b_{\kappa }^{+}\right),\sigma _{j}^{2}\sim IG(a_{\sigma },b_{\sigma })\). Note that α_{ i } and μ_{ j } are not identifiable because α_{ i }+μ_{ j }=(α_{ i }−c)+(μ_{ j }+c) for any constant c. If desired, the identifiability problem can be resolved by setting \(\sum _{i=1}^{n}\alpha _{i}=0\). For posterior simulation, the constraint can be implemented by setting α_{ i }←α_{ i }−c and μ_{ j }←μ_{ j }+c with \(c=\sum _{i=1}^{n}\alpha _{i}\) at each iteration.
The dependence of gene expressions is modeled through the prior model for the latent variables e_{ ij }. In [29], they further introduced another set of latent Gaussian variables z_{ ij } as probit scores for each trinary e_{ ij },
The continuous latent probit scores are then modeled by an SEM
with \(\epsilon _{ij}\sim N\left (0,s_{j}^{2}\right)\) and \(m_{ij}=\boldsymbol {x}_{i}^{T}\beta _{j}\) for covariates x_{ i }. Conjugate normal priors are assumed for regression coefficients β_{ j } and a_{ jk }. By (3) this defines an RGM. Let A be the p×p matrix of which the diagonal elements are zeros and the offdiagonal (j,k) entries are a_{ jk }. Then the joint distribution of Z_{ i }=(z_{i1},…,z_{ ip })^{T} is given by
where m_{ i }=(m_{i1},…,m_{ ip })^{T} and S=diag(s_{1},…,s_{ p }). The structure of the path diagram \(\mathcal {G}\) for this SEM is determined by the zero patterns in A, i.e. a_{ jk }=0 if and only if there is a missing edge from k to j. They define the prior over graph \(\mathcal {G}\) to be a subgraph of some fixed graph \(\mathcal {G}_{0}=(V,E_{0})\). The idea is to specify \(\mathcal {G}_{0}\) as some known maximum pathway. That is, \(p(\mathcal {G}\psi)=\psi ^{E}(1\psi)^{E_{0}E}\) where ψ∼Beta(a_{ ψ },b_{ ψ }). Restricting \(\mathcal {G}\) to subgraphs of \(\mathcal {G}_{0}\) very effectively incorporates the prior information that might be available about an established biological pathway, using, for example, a public database such as KEGG. This restriction to subgraphs of \(\mathcal {G}_{0}\) is important to mitigate nonidentifiability due to Markov equivalence between distinct graphs. Moreover, it also greatly reduces the graph search space, which allows for efficient posterior simulation through transdimensional Markov chain Monte Carlo (MCMC). Further discussion and detailed treatment of this topic can be found in [29].
The network structure in [29] is built on the latent Gaussian variables which are linked to the observed gene expressions through another layer of latent trinary variables. Therefore it defines an indirect dependence structure on the actual gene expression. [30] extend this work to allow for a more direct characterization of the dependence structure of the observed protein expressions. Particularly, they proposed a fully general dependence prior between latent binary variables and developed a parsimonious framework for Bayesian model determination that allows for local computation.
Telesca et al. [30] start with the same sampling model as in (5) and (6), but then deviate substantially from the earlier model by reducing the trinary variable e_{ ij } to a binary indicator variable z_{ ij }=2I(e_{ ij }=1)−1. That is,
The motivation of this reduction is the simplicity of the prior model on the resulting binary vector Z_{ i }=(z_{ ij }). In fact, it is possible [8] to describe all possible probability models p(Z_{ i }) that comply with a given conditional independence structure. See below. Also, investigators usually focus on protein activations in RPPA studies. The joint distribution of Z_{ i }=(z_{i1},…,z_{ ip })^{T} conditioned on the graph \(\mathcal {G}\) is by an Ising model,
The model is completed with normal priors for α_{ j } and β_{ jk }. In essence, the approach first maps an RG \(\mathcal {G}=(V,E)\) to its moral graph \(\mathcal {G}^{m}=(V,E^{m})\) and then defines the joint probability model for Z_{ i } based on \(\mathcal {G}^{m}\). The prior model for the graph \(\mathcal {G}\) is again hinged upon a known reference network \(\mathcal {G}_{0}=(V,E_{0})\). However, it is not constrained to be a subgraph of the reference network as before. Instead, they assume \(p(\mathcal {G}\psi)\propto \psi ^{d(\mathcal {G},\mathcal {G}_{0})}\) with ψ∼Beta(a_{ ψ },b_{ ψ }) where d(·,·) is a discrepancy measure. They choose \(d(\cdot,\cdot)=E^{c}\cap E_{0}+\delta E\cap E_{0}^{c}\) with δ>1 which is a weighted sum of edges dropped from and added to the reference graph \(\mathcal {G}_{0}\). With δ>1, the discrepancy measure allows for parsimonious inference by imposing a heavier penalty for adding edges than removing edges. A distinctive feature of this construction is that the prior decreases to zero exponentially fast as discrepancy measure increases, which allows the posterior simulation to be concentrated around the reference network. See [30] for more modeling and inference details of this topic.
Integrative network analysis from multiplatform genomic data
Different types of molecules in a living cell do not act in isolation; in fact, DNA, RNA and proteins work closely together to carry out various functions for each cell. Studies that consider only one specific type of molecules remain unnecessarily restricted. In [1] they integrate DNA and RNA molecules for a more robust and biologically interpretable estimate of a gene network. In particular, they consider p gene expressions Y=(Y_{1},…,Y_{ p })^{T}, together with corresponding copy number and methylation, collectively denoted by X=(X_{1},…,X_{2p})^{T}. They then exploit the central dogma of molecular biology that gene expression is produced by transcription from segments of DNA on which the copy number and methylation are measured, but the reverse processes are rare and biologically uninterpretable. The importance of this restriction is that in a network of Z=(Y,X) some edges have fixed predetermined directions corresponding this observation. This in turn allows us to report meaningful inference on the directions of the remaining edges. Similar idea but in a different context has been proposed by [32–34]. With this assumption, they then model the regulatory relationships between copy number, methylation and gene expression with the following SEM
where \(\boldsymbol {A}=(a_{ij})\in \mathbb {R}^{p\times p}\) with zeros on the diagonal, \(\boldsymbol {B}=(b_{ij})\in \mathbb {R}^{p\times 2p}\) and E=(ε_{1},…,ε_{ p })^{T}∼N_{ p }(0,Σ). The SEM in (7) immediately prohibits gene expression from regulating copy number or methylation since there is no X on the righthand side of the equation. It also explicitly allows for feedback loops between the genes which are quite common motifs in molecular networks and have key functional roles in many cellular processes such as regulating gene expressions and acting as bistable switches [35]. Such feature is missing in DAG’s, UG’s and CG’s.
Next they assign thresholded priors [36] for A and B,
for i=1,…,p, j≠i and k=2i−1,2i. The threshold parameter t_{ i } controls the minimum effect sizes of a_{ ij } and b_{ ik } and is assigned a uniform hyperprior t_{ i }∼Uniform(0,t_{0}). Normal priors are given to \(\tilde {a}_{ij}\) and \(\tilde {b}_{ik}\), \(\tilde {a}_{ij}\sim \textrm {N}\left (0,\tau _{ij}^{(a)}\right)\) and \(\tilde {b}_{ik}\sim \textrm {N}\left (0,\tau _{ik}^{(b)}\right)\) and conjugate hyperpriors \(\tau _{ij}^{(a)}\sim \text {IG}(\alpha _{\tau },\beta _{\tau })\) and \(\tau _{ik}^{(b)}\sim \text {IG}(\alpha _{\tau },\beta _{\tau })\). The thresholded priors enjoy nice theoretical properties as they are closely related to spikeandslab prior and nonlocal prior as shown in [36–38]. We refer readers to [1] for a more comprehensive description of the work.
Reverse engineering gene networks with heterogeneous samples
Genomic data are often heterogeneous in the sense that the samples can be naturally divided into observed or hidden groups. For example, in a pancancer study, patients are normally grouped by the known cancer types whereas when studying a specific cancer, patients can be split into homogeneous latent subtypes. The inferential goals of these two cases are typically quite different. In the case with observed groups, much of the statistical work has been devoted to building models that allow information to be shared across groups. Such approaches improve the statistical power in graph estimation especially for groups with very limited samples by borrowing strength from other groups. For hidden groups, the underlying task becomes a clustering problem which aims to partition subjects into homogeneous clusters. In this section, we introduce two hierarchical RG models that are suitable for heterogeneous datasets. Interested readers can find more details in [38].
Suppose there are C known groups in the data. Let {s_{ i }=c} be the group indicator for patient i in group c∈{1,…,C}. Denote Y_{ c }={Y_{ i }; s_{ i }=c} and X_{ c }={X_{ i }; s_{ i }=c} the set of observations for patients in group c. Then [38] entertain a groupspecific SEM
where A_{ c } = (a_{c,ij}),B_{ c } = (b_{c,ij})and E_{ c } = (ε_{c,1},…,ε_{c,p})^{T} ∼ N_{ p }(0,Σ_{ c }) with \(a_{c,ij}=\tilde {a}_{c,ij}\textrm {I}(\mid \tilde {a}_{c,ij}\mid >t_{i})\) and \(b_{c,ik}=\tilde {b}_{c,ik}\textrm {I}(\mid \tilde {b}_{c,ik}\mid >t_{i})\). To link the graphs across groups, they impose multivariate normal priors on the edge strength
The covariance matrix \(\boldsymbol {\Omega }\!\!\!\,=\!(\omega _{cc^{\prime }})\) connects edge strength and graph structure across groups. When \(\omega _{cc^{\prime }}\) is close to ±1, \(\tilde {a}_{c,ij}\) and \(\tilde {a}_{c^{\prime },ij}\) tend to have similar absolute values, which gives rise to a high probability that the edge i←j is included or excluded in groups c and c^{′} simultaneously since the threshold t_{ i } is shared across groups. Therefore, graphs from group c and c^{′} are more likely to share common edges. On the other hand, when \(\omega _{cc^{\prime }}\) is close to 0, the association between groups c and c^{′} is negligible. They do not fix Ω. Instead, they use an inverseWishart prior, Ω∼IW(ν,Ψ) and let the data dictate how strong the associations should be.
The Bayesian hierarchical model above allows for easy extension to the case where the groups are unknown. Ni et al. [38] simply augment the model with a Dirichletmultinomial (DM) allocation model on the group indicator cπ,C∼Multinomial(1,π_{1},…,π_{ C }) and πC∼Dir(η,…,η). The number of clusters C is usually unknown; they assume a geometric prior for C∼Geo(ρ), which eliminates the need to fix C a priori. They set the covariance matrix Ω to be diagonal since the goal of clustering is to split samples into groups with disparate networks rather than encourage similarity across the clusters. Alternatively, one can also set Ω to be a Stieltjes matrix (i.e. a positive definite matrix with nonpositive offdiagonal entries) to induce repulsive graphs.
Results
In this section, we generalize the RG approach of [1] by incorporating protein expressions in modeling the molecular network. Proteins provide important and orthogonal information in reverse engineering the network because they represent the downstream cumulative effect of changes that happen at the DNA and RNA levels and are more directly related to the phenotypical changes in cancer cells. Clinical utilization of genomic data alone show limited benefits, which is partly due to the poor concordance between gene and protein expressions [39, 40]. Many factors such as complex interactions between different types of molecules are responsible for the discrepancy between gene and protein abundance. The goal of this analysis is to explicitly model such complex interactions. We choose RG because it allows for efficient computation and great interpretability and is suitable to model reciprocal causality such as feedback mechanism.
We use the same notations Y=(Y_{1},…,Y_{ p })^{T} and X=(X_{1},…,X_{2p})^{T} to denote gene expressions and the matched copy number and methylation. Importantly, additionally, for each gene Y_{ j }, there is also the associated protein expression, denoted by Z_{ j } for j=1,…,p. We extend the SEM in (7) to
where \(\boldsymbol {A}=(a_{ij})\in \mathbb {R}^{p\times p}\) and \(\boldsymbol {C}=(c_{ij})\in \mathbb {R}^{p\times p}\) have zeros on the diagonal, \(\boldsymbol {B}=(b_{ij})\in \mathbb {R}^{p\times 2p}\), \(\boldsymbol {D}=(d_{ij})\in \mathbb {R}^{p\times p}\) and independent errors E_{ Y }∼N_{ p }(0,Σ_{ Y }) and E_{ Z }∼N_{ p }(0,Σ_{ Z }) with diagonal covariance matrices Σ_{ Y }=diag(σ_{1},…,σ_{ p }) and Σ_{ Z }=diag(λ_{1},…,λ_{ p }). In essence, model (8) defines the following joint probability model for (Y,Z) given X,
The network structure is embedded in the parameters A,B,C and D. In this analysis, we restrict our attention to cisregulatory relationships between DNAs and RNAs by constraining b_{ ij } to 0 for j≠2i−1 or 2i for i=1,…,p. That is, copy number and methylation of gene i can only affect gene expression i. The interpretation of model (9) is given as follows: b_{i,2i−1}≠0 if and only if copy number i is associated with its gene expression; b_{i,2i}≠0 if and only if methylation i is associated with its gene expression; a_{ ij }≠0 if and only if gene j regulates gene i; c_{ ij }≠0 if and only if protein j regulates protein j; and d_{ ij }≠0 if and only if gene j regulates protein i. We use the thresholded prior for A,B,C and D. For error variances, we assume inversegamma priors, σ_{ j },λ_{ j }∼IG(a,b) for j=1,…,p which imply an inverse gamma conditional posterior distribution which in turn allows the use of a Gibbs sampling transition probability by sampling from the respective inverse gamma. The complete MCMC algorithm is provided in Algorithm 1.
Using TCGAAssembler ([41], version 2), we acquired, processed and combined TCGA ovarian cancer data, including copy number, methylation and gene expression from the Genomic Data Commons (GDC) and mass spectrometry proteomic samples generated by the Clinical Proteomic Tumor Analysis Consortium (CPTAC). The resulting dataset contains 104 samples. In this case study, we focus on p=10 genes that are core members of the PI3K pathway [42], as shown in Fig. 3a, by playing a critical role in cell cycle progression, survival, motility, angiogenesis and immune surveillance [43]. It is one of the most frequently altered pathways in ovarian cancer [44, 45]. For each gene, we have matched copy number, methylation, gene expression and protein expression (i.e. 40 nodes in total). We report inference based on 100,000 MCMC samples (thinned out to every 10th iteration) after discarding the first 50,000 iterations as burnin.
The recovered network is shown in Fig. 3b where the blue solid lines represent activation whereas the red dashed lines are inactivations. The edge width is proportional to the posterior probabilities of inclusion p(a_{ ij }≠0X,Y,Z), p(b_{ ij }≠0X,Y,Z), p(c_{ ij }≠0X,Y,Z) and p(d_{ ij }≠0X,Y,Z), which can be approximated by MCMC sample averages, for example, \(p(a_{ij}\neq 0\boldsymbol {X},\boldsymbol {Y},\boldsymbol {Z})=\frac {1}{L}\sum _{l=1}^{L}I\left (a_{ij}^{(l)}\neq 0\right)\) where the superscript (l) denotes the lth MCMC sample and L is the total number of samples. The suffixes represent four different type of molecules: c=copy number, m=methylation, g=gene and p=protein. For clarity we do not show disconnected molecules. The degree of each node in Fig. 3b is reported in Table 1. The top three highly connected nodes are the protein AKT1 which is involved in 15 regulatory activities, the protein mTOR with 10 activities and the protein AKT3 with 9 activities. AKT1 and AKT3 belong to the same AKT family which is known to have shown strong oncogenic function and is a key mediator of PI3K pathway function [46]. AKT isoforms are often phosphorylated in ovarian cancers and may play a role in mediating the progression of latestage serous ovarian carcinomas [47]. mTOR plays a critical role in regulating cell growth and proliferation by integrating signals including growth factors, nutrients, energy and stress [48, 49]. We also capture the wellknown positive regulatory relationship between AKT1 and its downstream target mTOR [45, 49]. Other findings that are consistent with the biological literature include the inhibitory relationship between proteins PTEN and AKT1 and stimulatory relationship between proteins PIK3CA and AKT1 [50, 51]. Interestingly, the recovered network includes a negative feedback from proteins of mTOR to AKT1. Recent studies [45, 52–54] confirmed that such feedback mechanism can be either positive or negative depending on many factors including cell types and conditions.
Inference includes some surprising findings. For example, PIK3R1 is expected to activate rather than inactivate AKT1. Although the unexpected result may simply be a false positive, it deserves further experimental validation.
Discussion and conclusions
We have reviewed the basics of UG’s, DAG’s and RG’s and discussed recent RGbased approaches in modeling molecular networks. The extension of the approach in [1] allows to integrate multiplatform data including copy number, methylation, gene expression and protein expression. Our approach can be further be generalized to incorporating other data types such as single nucleotide polymorphisms (SNPs), mutation status and microRNAs. SEMbased RG approaches, however, cannot be directly applied to SNPs and mutation status because they have discrete support and the condition C1 requires the data to be multivariate normal. This limitation can be addressed by imputing latent normal variables. Although we do not find serious violation of normality assumption in our data, if desired, one could adopt the approach of [29] and introduce additional layers of hidden variables when the data are seemingly nonnormal.
Another important assumption of the discussed approaches in this paper is homogeneity across samples. It is an unrealistic assumption in many applications, especially in oncology where the consensus is that tumors are extremely heterogeneous. Characterizing tumor heterogeneity can fundamentally improve our understanding of the cancer biology and practically allow us to divide heterogeneous patient population into homogeneous subpopulations so that refined personalized treatments can be developed targeting specific subgroups of cancer patients. Several approaches have been proposed for modeling networks with heterogeneous samples [38, 55–57].
References
Ni Y, Ji Y, Müller P. Reciprocal Graphical Models for Integrative Gene Regulatory Network Analysis. Bayesian Anal. 2017.
Lauritzen SL, Wermuth N. Graphical models for associations between variables, some of which are qualitative and some quantitative. Ann Stat. 1989; 17(1):31–57.
Frydenberg M. The chain graph markov property. Scand J Stat. 1990; 17(4):333–53.
Spirtes P. Directed Cyclic Graphical Representations of Feedback Models. In: Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence. San Francisco: Morgan Kaufmann Publishers Inc.: 1995. p. 491–8.
Koster JT. Markov properties of nonrecursive causal models. Ann Stat. 1996; 24(5):2148–77.
Whittaker J. Graphical models in applied multivariate statistics. New York: Wiley Publishing; 2009.
Lauritzen SL, Vol. 17. Graphical models. Oxford: Clarendon Press; 1996.
Besag J. Spatial interaction and the statistical analysis of lattice systems. J R Stat Soc Series B. 1974; 36(2):192–236.
Dobra A, Lenkoski A, Rodriguez A. Bayesian inference for general gaussian graphical models with application to multivariate lattice data. J Am Stat Assoc. 2011; 106(496):1418–33.
Green PJ, Thomas A. Sampling decomposable graphs using a markov chain on junction trees. Biometrika. 2013; 100(1):91.
Wang H, Li SZ, et al. Efficient gaussian graphical model determination under gwishart prior distributions. Electron J Stat. 2012; 6:168–98.
Meinshausen N, Bühlmann P. Highdimensional graphs and variable selection with the lasso. Ann Stat. 2006; 34(3):1436–62.
Friedman J, Hastie T, Tibshirani R. Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 2008; 9(3):432–41.
Pearl J. Probabilistic inference in intelligent systems. San Mateo: Morgan Kaufmann; 1988.
Lauritzen SL, Dawid AP, Larsen BN, Leimer HG. Independence properties of directed markov fields. Networks. 1990; 20(5):491–505.
Ni Y, Stingo FC, Baladandayuthapani V. Sparse MultiDimensional Graphical Models: A Unified Bayesian Framework. J Am Stat Assoc. 2017; 112(518):779–93.
Ni Y, Stingo FC, Baladandayuthapani V. Bayesian nonlinear model selection for gene regulatory networks. Biometrics. 2015; 71(3):585–95.
Friedman N, Koller D. Being bayesian about network structure. a bayesian approach to structure discovery in bayesian networks. Mach Learn. 2003; 50(12):95–125.
Shojaie A, Michailidis G. Penalized likelihood methods for estimation of sparse highdimensional directed acyclic graphs. Biometrika. 2010; 97:519–38.
Stingo FC, Chen YA, Vannucci M, Barrier M, Mirkes PE. A bayesian graphical modeling approach to microrna regulatory network inference. Ann Appl Stat. 2010; 4(4):2024.
Altomare D, Consonni G, La Rocca L. Objective bayesian search of gaussian directed acyclic graphical models for ordered variables with nonlocal priors. Biometrics. 2013; 69(2):478–87.
Yajima M, Telesca D, Ji Y, Müller P. Detecting differential patterns of interaction in molecular pathways. Biostatistics. 2015; 16(2):240–51.
Richardson T. A Discovery Algorithm for Directed Cyclic Graphs. In: Proceedings of the Twelfth International Conference on Uncertainty in Artificial Intelligence. San Francisco: Morgan Kaufmann Publishers Inc.: 1996. p. 454–61.
Richardson T. A Polynomialtime Algorithm for Deciding Markov Equivalence of Directed Cyclic Graphical Models. In: Proceedings of the Twelfth International Conference on Uncertainty in Artificial Intelligence. San Francisco: Morgan Kaufmann Publishers Inc.: 1996. p. 462–69.
Richardson T. A characterization of markov equivalence for directed cyclic graphs. Int J Approx Reason. 1997; 17(23):107–62.
Mooij JM, Janzing D, Heskes T, Schölkopf B. On Causal Discovery with Cyclic Additive Noise Models. In: Proceedings of the 24th International Conference on Neural Information Processing Systems. USA: Curran Associates Inc.: 2011. p. 639–47.
Mooij JM, Heskes T. Cyclic Causal Discovery from Continuous Equilibrium Data. In: Proceedings of the TwentyNinth Conference on Uncertainty in Artificial Intelligence. Arlington: AUAI Press: 2013. p. 431–9.
Lacerda G, Spirtes P, Ramsey J, Hoyer P. Discovering Cyclic Causal Models by Independent Components Analysis. In: Proceedings of the TwentyFourth Conference Annual Conference on Uncertainty in Artificial Intelligence (UAI08). Corvallis: AUAI Press: 2008. p. 366–74.
Telesca D, Müller P, Parmigiani G, Freedman RS. Ann Appl Stat. 2012; 6(2):542–60.
Telesca D, Müller P, Kornblau SM, Suchard MA, Ji Y. Modeling protein expression and protein signaling pathways. J Am Stat Assoc. 2012; 107(500):1372–84.
Parmigiani G, Garrett ES, Anbazhagan R, Gabrielson E. A statistical framework for expressionbased molecular classification in cancer. J R Stat Soc Series B. 2002; 64(4):717–36.
Cai X, Bazerque JA, Giannakis GB. Inference of gene regulatory networks with sparse structural equation models exploiting genetic perturbations. PLoS Comput Biol. 2013; 9(5):1003068.
Ni Y, Stingo FC, Baladandayuthapani V. Integrative bayesian network analysis of genomic data. Cancer Inf. 2014; 13(Suppl 2):39.
Zhang L, Kim S. Learning gene networks under snp perturbations using eqtl datasets. PLoS Comput Biol. 2014; 10(2):1003420.
Shin SY, Rath O, Zebisch A, Choo SM, Kolch W, Cho KH. Functional roles of multiple feedback loops in extracellular signalregulated kinase and wnt signaling pathways that regulate epithelialmesenchymal transition. Cancer Res. 2010; 70(17):6715–24.
Ni Y, Stingo F, Baladandayuthapani V. Bayesian Graphical Regression. J Am Stat Assoc. 2018.
Rossell D, Telesca D. Nonlocal Priors for HighDimensional Estimation. J Am Stat Assoc. 2017; 112(517):254–65.
Ni Y, Müller P, Zhu Y, Ji Y. Heterogeneous reciprocal graphical models. Biometrics. 2018.
Gygi SP, Rochon Y, Franza BR, Aebersold R. Correlation between protein and mrna abundance in yeast. Mol Cellular Biol. 1999; 19(3):1720–30.
Akbani R, Ng PKS, Werner HM, Shahmoradgoli M, Zhang F, Ju Z, Liu W, Yang JY, Yoshihara K, Li J, et al. A pancancer proteomic perspective on The Cancer Genome Atlas. Nat Commun. 2014; 5:3887.
Zhu Y, Qiu P, Ji Y. Tcgaassembler: opensource software for retrieving and processing tcga data. Nat Methods. 2014; 11(6):599–600.
Cantley LC. The phosphoinositide 3kinase pathway. Science. 2002; 296(5573):1655–57.
Liliac L, Amalinei C, Balan R, Grigoras A, Caruntu ID. Ovarian cancer: insights into genetics and pathogeny. Histol Histopathol. 2012; 27(6):707–19.
Kotsopoulos IC, Papanikolaou A, Lambropoulos AF, Papazisis KT, Tsolakidis D, Touplikioti P, Tarlatzis BC. Serous ovarian cancer signaling pathways. Int J Gynecol Cancer. 2014; 24(3):410–7.
Cheaib B, Auguste A, Leary A. The pi3k/akt/mtor pathway in ovarian cancer: therapeutic opportunities and challenges. Chin J Cancer. 2015; 34(1):4.
Bast RC, Hennessy B, Mills GB. The biology of ovarian cancer: new opportunities for translation. Nat Rev Cancer. 2009; 9(6):415–28.
Hanrahan AJ, Schultz N, Westfal ML, Sakr RA, Giri DD, Scarperi S, Janikariman M, Olvera N, Stevens EV, She QB. Genomic complexity and akt dependence in serous ovarian cancer. Cancer Discov. 2012; 2(1):56–67.
Wullschleger S, Loewith R, Hall MN. Tor signaling in growth and metabolism. Cell. 2006; 124(3):471–84.
Liu P, Cheng H, Roberts TM, Zhao JJ. Targeting the phosphoinositide 3kinase pathway in cancer. Nat Rev Drug Discov. 2009; 8(8):627–44.
Vara JÁF, Casado E, de Castro j, Cejas P, BeldaIniesta C, GonzálezBarón M. Pi3k/akt signalling pathway and cancer. Cancer Treat Rev. 2004; 30(2):193–204.
Song MS, Salmena L, Pandolfi PP. The functions and regulation of the pten tumour suppressor. Nat Rev Mol Cell Biol. 2012; 13(5):283–96.
Hay N. The aktmtor tango and its relevance to cancer. Cancer Cell. 2005; 8(3):179–83.
Dormond O, Madsen JC, Briscoe DM. The effects of mtorakt interactions on antiapoptotic signaling in vascular endothelial cells. J Biol Chem. 2007; 282(32):23679–86.
Mead H, Zeremski M, Guba M. mTOR Signaling in Angiogenesis In: Polunovsky VA, Houghton PJ, editors. mTOR Pathway and mTOR Inhibitors in Cancer Therapy. Totowa: Humana Press: 2010. p. 49–74.
Ickstadt K, Bornkamp B, Grzegorczyk M, Wieczorek J, Sheriff MR, Grecco HE, Zamir E. Nonparametric bayesian networks. Bayesian Stat 9. 2011; 9:283.
Rodriguez A, Lenkoski A, Dobra A. Sparse covariance estimation in heterogeneous samples. Electron J Stat. 2011; 5:981.
Mukherjee C, Rodriguez A. Gpupowered shotgun stochastic search for dirichlet process mixtures of gaussian graphical models. J Comput Graph Stat. 2016; 25(3):762–88.
Acknowledgements
The authors appreciate the instructive suggestions from the Editor and anonymous reviewers.
Funding
This work and the publication costs of this article were supported by NIH/NCI grant RO1 CA 132897.
Availability of data and materials
The TCGA ovarian cancer dataset is publicly available in Genomic Data Commons (https://gdc.cancer.gov/).
About this supplement
This article has been published as part of BMC Bioinformatics Volume 19 Supplement 3, 2018: Selected original research articles from the Fourth International Workshop on Computational Network Biology: Modeling, Analysis, and Control (CNBMAC 2017): bioinformatics. The full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume19supplement3.
Author information
Authors and Affiliations
Contributions
YN, PM and YJ conceived the idea of this article. YN carried out the literature reviews, extended the model and drafted the manuscript. PM and YJ made critical revisions. LW and YJ provided the data. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Ni, Y., Müller, P., Wei, L. et al. Bayesian graphical models for computational network biology. BMC Bioinformatics 19 (Suppl 3), 63 (2018). https://doi.org/10.1186/s128590182063z
Published:
DOI: https://doi.org/10.1186/s128590182063z
Keywords
 Directed graph
 Undirected graph
 Chain graph
 Reciprocal graph
 Causality