Introduction

Causal modeling is a useful analytical tool in various fields due to its applicability in action planning, prediction and diagnosis [1,2,3,4,5]. However, learning a causal Bayesian network (CBN) solely from data is a challenging task [6,7,8]. CBN learning can be thought of as model selection problem in which model is a directed acyclic graph (DAG), where problem is to find the graph 𝐺 that maximizes some objective (score) function of dataset 𝐷. In some CBN learning methods, score function is likelihood 𝑝(𝐷 | 𝐺) representing overall fit of 𝐺 to 𝐷 in context of a generative model for the data. For a dataset with 𝑛 observables (features), the number of possible DAGs-and thus requirement for data-grows super-exponentially with 𝑛 [9]. In most network learning applications, prior knowledge exists about causal (or suspected causal) relationships among the observables; such prior knowledge can be valuable resource for network structure learning [10]. Supposing that prior knowledge can be represented as prior probability 𝑝(𝐺) on network structure, one can alternatively choose as basis for CBN scoring function, the posterior probability 𝑝(𝐺 | 𝐷) = 𝑝(𝐷 | 𝐺) 𝑝(𝐺)/𝑝(𝐷). In contrast to substantial amount of work done on variety of marginal likelihood and scoring methods, less attention has been given to functional form (and associated parameterization) of the prior 𝑝(𝐺) for application contexts where structured prior knowledge is available. Without expert knowledge, standard network inference approaches, by default, assume uniform (uninformative) prior which can lead to erroneous relationships or relationship orientations both due to (i) size of space of networks and (ii) degeneracy of Markov-equivalent networks. Proper incorporation of informative priors can enhance model efficiency [11] and can also overcome weakness of smaller dataset.

For most applications of causal modeling, some prior knowledge is available. For example, in medicine, most cases have prior knowledge about etiology, symptoms, and treatment of underlying diseases or conditions which can be obtained from biomedical literature or knowledgebases. Although there is in general large scale availability of structured prior knowledge (for example, ontologies) in various scientific domains, these mostly comprise disparate information sources in various standards and formats, which poses a challenge to integrate them into single structure. These problems motivated building of large multi-graphs called knowledge graphs (KG) [12] that incorporate structured knowledge from multiple sources within a consistent schema. Knowledge graph is a term of art to mean a large graph-structured model to store interlinked relationships between nodes representing concepts [13]. These large-scale networks accommodate structural information which can be leveraged for reasoning, recommendation or decision making. We hypothesized that combining information from structured databases of general prior knowledge with causal modeling based on context-specific multivariate measurements will improve accuracy of learned network compared to result of datadriven causal modeling without incorporating prior knowledge. In this work, we propose a method, "Kg2Causal", for extracting relations as pairs of nodes from a knowledge graph, and for incorporating them as priors on corresponding edges in a score-based, data-driven causal network learning method. In this study, prior edges from knowledge graph are accounted for in prior probability of the graph, using the method of Castelos and Siebes [11]. We used a large-scale biomedical knowledge graph (KG) 1 that we and collaborators (see Acknowledgments) had previously constructed (see Sec. 2.5) containing millions of nodes representing drugs, genes, diseases, or phenotypes (Fig. 1b), as well as edges between nodes representing various types of predicate relationships. For the measurement-based network learning component of KG2Causal, we used an optimizing method combining Tabu search algorithm [14] with Bayesian Dirichlet uniform (BDeu) [15,16] network score. Using five different multivariate molecular biology datasets for which groundtruth networks were available (see Sec. 4), we empirically analysed network learning accuracy of "Kg2Causal" along with various types of uninformative prior to test the usefulness of adding KG-based priors. We provide a comparative benchmark of our methods performance over five real-world biological datasets and two synthetic datasets of varying sizes where we found that Kg2Causal had superior network learning accuracy to methods that do not use general knowledge-base as network structure prior. Finally, we demonstrate (Sec. 4.3) the application of Kg2Causal if data are collected under different experimental conditions including interventions on the observed variables. We implemented "Kg2Causal" in the R programming language (leveraging the bnLearn package [17]) and provide the code as open-source software2 .

Related work and Background

In this section, we describe Kg2Causal's conceptual foundations including CBNs, score-based causal modeling, interventions, and knowledge graph-based priors in network learning.

Causal network: Brief Overview

A causal Bayesian Network [1,2] is a DAG 𝐺 = (𝑉 , 𝐸), where 𝑉 = {𝑉 1 , … , 𝑉 𝑛 } denotes the set of variables (nodes) and 𝐸 ⊂ 𝑉 × 𝑉 denotes the causal relationships (edges). For an edge (𝑉 𝑖 , 𝑉 𝑗 ), we say that 𝑉 𝑖 is a parent (cause) of 𝑉 𝑗 , and 𝑉 𝑗 is a child (effect) of 𝑉 𝑖 . We will use Pa(𝑉 𝑖 ) to denote the set of parents of 𝑉 𝑖 . The conditional probability distribution 𝑝 𝑖 defines the probability of 𝑉 𝑖 given the states of its parents Pa(𝑉 𝑖 ). A causal network represents a joint distribution 𝑝 over variables 𝑉 as long as it satisfies two main assumptions: a) Causal Markov: Any given variable 𝑉 𝑖 is independent of its non-descendants, conditioned on all of its direct causes. The assumption implies that the joint distribution 𝑝(𝑉 ) can be factored as: 𝑝(𝑉 ) = ∏ 𝑛 𝑖=1 𝑝 𝑖 (𝑉 𝑖 | Pa(𝑉 𝑖 )). b) Faithfulness: The joint distribution 𝑝(𝑉 1 , … , 𝑉 𝑛 ) is faithful to 𝐺 if every conditional independence relation in 𝑝 is entailed by the causal Markov assumption applied to 𝐺 [18].

Constructing a causal network

Let us assume, we have a dataset 𝐷 having observations over set of 𝑛 variables. One of the main classes of causal learning approaches is the Score-based approach which is derived from classic Bayesian method where a scoring function evaluates the fit of graph 𝐺 to data 𝐷 [16,6] with a higher value indicating better fit. A search algorithm is used to explore the space of all possible graphs, to maximize the scoring function. Typical heuristic algorithms used for this purpose include hill-climbing or Tabu search approaches [14]. Other common score based methods are GDS [19] and Gies [6]. According to standard Bayesian rule, a causal graph 𝐺, is a 𝐷𝐴𝐺 learned from given data 𝐷 as 𝑝(𝐺 | 𝐷) ∝ 𝑝(𝐺)𝑝(𝐷 | 𝐺), where 𝑝(𝐺) is prior distribution over space of all possible DAGs reflecting prior knowledge and 𝑝(𝐷 | 𝐺) is marginal likelihood of the data 𝐷. As described in Sec. 3, the Kg2Causal method incorporates a score-based approach.

Learning with interventions

Interventions-external manipulations of nodes ("targets") in a network-are important to detect causal relations that can help disambiguate Markov equivalent sub-networks [16]. Let 𝐼 𝑒 represent set of target nodes that are altered in interventional experiment 𝑒 and 𝑂 𝑒 = 𝑉 \𝐼 𝑒 be the complementary set of observational variables. Each intervention can have one or more targets whose conditional probabilities are changed (so that, conditioned on intervention, target variable's distribution may depend only on a (possibly empty) subset of its parent observables). Hence, each intervention results in deletion of arrows pointing towards the intervened nodes. The joint distribution of 𝑝 after intervention is 𝑝(𝑉 1 , ..., 𝑉 𝑛 ) = ∏ 𝑉 𝑖 ∈𝑂 𝑒 𝑝 𝑖 (𝑉 𝑖 | 𝑃𝑎(𝑉 𝑖 ))⋅∏ 𝑉 𝑗 ∈𝐼 𝑒 𝑝 ′ 𝑗 (𝑉 𝑗 | 𝑃𝑎(𝑉 𝑗 )), where 𝑝(𝑉 𝑖 | 𝑃𝑎(𝑉 𝑖 )) is conditional probability similar to 𝑉 𝑖 , given that 𝑉 𝑖 is not a target node, and 𝑝 ′ (𝑉 𝑗 | 𝑃𝑎 ′ (𝑉 𝑗 )) is post-intervention conditional probability of 𝑉 𝑗 given its new set of parents 𝑃𝑎 ′ (𝑉 𝑗 ). For a so-called "perfect" intervention, one would set 𝑃𝑎 ′ (𝑉 𝑗 ) = ∅ [1]. Score-based approaches are well-suited to mixed interventional-observational datasets, in contrast to constraint-based approaches which are applicable to observational data.

Incorporation of Priors

In this subsection we introduce three types of uninformative priors on the network structure 𝑝(𝐺), uniform prior, marginal prior, and Bayesian variable selection prior (VSP). We then describe the knowledge graph-based prior that we use in Kg2Causal method. In cases for lack of prior knowledge, default choice for prior 𝑝(𝐺) is a uniform prior distribution, as follows:

𝑝(𝐸 ∪ {𝑉 𝑖 , 𝑉 𝑗 } | 𝐷) 𝑝(𝐸 | 𝐷) = 𝑝(𝐸 ∪ {𝑉 𝑖 , 𝑉 𝑗 }) 𝑝(𝐸) 𝑝(𝐷 | 𝐸 ∪ {𝑉 𝑖 , 𝑉 𝑗 }) 𝑝(𝐷 | 𝐸)

where nodes 𝑉 𝑖 and 𝑉 𝑗 can have three possible cases 𝑉 𝑖 ⇒V 𝑗 (representing (𝑉 𝑖 , 𝑉 𝑗 ) ∈ 𝐸),V 𝑖 ⇐V 𝑗 (representing (𝑉 𝑗 , 𝑉 𝑖 ) ∈ E) or 𝑉 𝑖 ⇎V 𝑗 (no arc) and each have equal probability of occurrence. So the probability for these edges are assigned as 𝑝(𝑉 𝑖 ⇒V 𝑗 ) = 𝑝(𝑉 𝑖 ⇐V 𝑗 ) = 𝑝(𝑉 𝑖 ⇎ 𝑉 𝑗 ) = 1/3, since we know that 𝑝(𝑉 𝑖 ⇒V 𝑗 ) + 𝑝(𝑉 𝑖 ⇐ 𝑉 𝑗 ) + 𝑝(𝑉 𝑖 ⇎ 𝑉 𝑗 ) = 1. This implies 𝑝(𝑉 𝑖 ⇒ 𝑉 𝑗 ) + 𝑝(𝑉 𝑖 ⇐ 𝑉 𝑗 ) = 2/3, which means a higher promotion for the inclusion of new arcs and favouring the propagation of false positives in 𝐺. Hence, its not always a good idea to use uniform prior specially for cases where data is not too supportive of the DAG learned and where 𝑛 is large. A better version of uniform prior is to use marginal probabilities instead, where an independent prior can be assumed for each arc with same independent marginal probabilities as uniform priors, also called marginal uniform [20]. In this case, the probability of inclusion of each edge is assigned as 𝑝(𝑉 𝑖 ⇒ 𝑉 𝑗 ) = 𝑝(𝑉 𝑖 ⇐ 𝑉 𝑗 ) = 1/4 and 𝑝(𝑉 𝑖 ⇎ 𝑉 𝑗 ) = 1/2. Compared to the uniform prior, the marginal uniform prior is less prone to false-positive edges in the posterior-probability-maximizing graph. The Bayesian variable selection prior (VSP) assigns a probability of inclusion of possible parent nodes, with the default being 1/𝑛. The heart of Kg2Causal is the use of an informative prior based on a general-purpose knowledge graph; for this purpose we use an edge decomposition technique described by Castelo and Siebes [11]. For any pair of vertices (𝑉 𝑖 , 𝑉 𝑗 ) for which an edge 𝑉 𝑖 ⇒V 𝑗 exists in the general-purpose knowledge graph, we assign a prior probability (𝛽 = 1/2) on those edges, with probability 1/4 for 𝑉 𝑗 ⇒V 𝑖 and probability 1/4 for 𝑉 𝑖 ⇎V 𝑗 , since the later two are alternate edges that have no corresponding edge in the general knowledge graph we use the uniform probability distribution as shown in Fig. 2. In this way we can create a complete prior probability (from partial knowledge) over the network 𝐺; on log scale, we define 𝑝(𝐺) as log 𝑝(𝐺) = ∑ 𝑉 𝑖 ⇔𝑉 𝑗 ∈𝐸, 𝑖≠𝑗 log 𝑝(𝑉 𝑖 ⇔ 𝑉 𝑗 ) + ∑ 𝑉 𝑖 ⋯𝑉 𝑗 ∈𝐸, 𝑖≠𝑗 log 𝑝(𝑉 𝑖 ⇎ 𝑉 𝑗 ). A "knowledge graph" [13] is a multigraph consisting of nodes and edges (labeled by relationship type or description of instance attributes) between them. Although most relationships in knowledge graphs are between entities and context-based associations, these do not always imply causal relationship. Nevertheless, such links are strong association that can strengthen causal relationships that we seek to discover. The key idea of Kg2Causal is to use links from large knowledge graphs as generalised prior information to aid in data-driven network learning in highly specific application contexts. For this work, we leveraged a general biomedical knowledge graph that we and collaborators (see Acknowledgments) had constructed, KG1 3 . KG1 has 130,443 nodes, 3.5M edges, 11 node semantic types, and 17 edge relation types, and was compiled from 20 different biomedical knowledge-bases (Monarch, COHD, ChEMBL, DGIdb, DisGeNet, Disease Ontology, GeneProf, HMDB, KEGG, miRBase, miRGate, mychem.info, mygene.info, NCBI Gene, OMIM, Pathway Commons, Pharos, PubChem, Reactome, and UniprotKB). We hosted KG1 in a Neo4j database (ver. 3.5.13) and used the Cypher query language to search for concept mappings between ground-truth network variables and concept nodes in the KG1 knowledge graph, and for edge connections between mapped concepts within KG1.

Knowledge Graphs

Our Approach

We developed KG2Causal to leverage a general-purpose biomedical knowledge graph (see Sec. 2.5) in order to improve context-specific, data-driven network learning from multivariate observations; such observations could consist of gene expression measurements, proteomics measurements, or electronic health records. The key ideas of our approach are (i) mapping each variable in the dataset to a node in the knowledge graph, and querying relationships between them; (ii) extracting a subgraph containing the connected variables with edges between them; and (iii) use this edge set as our prior knowledge to guide the optimizing scoring step for inferring causal network. Mathematically, given a dataset 𝐷, with a set 𝑉 of observable variables and given a general-purpose prior knowledge graph Γ as a multigraph, we want to learn a causal graph 𝐺 𝑓 = (𝑉 , 𝐸) that approximately maximizes the posterior probability, i.e., argmax 𝐺 (𝑝(𝐺 | 𝐷, Γ)), given a prior 𝑝(𝐺 | Γ). As a comparison, we used three uninformative prior distributions, namely uniform, marginal and Bayesian variable selection priors with each dataset in order to understand whether or not-and to what extent-using an informative network prior improves accuracy of causal network learning in a biomedical context. The Kg2Causal network discovery workflow, illustrated in Figure 3, consists of the following steps:

• Map variables 𝑉 to nodes in Γ, and extract a list 𝛽 of edges from Γ among the nodes (collapsing same-direction multiedges to single edges). Step 2 -Cleaning and discretizing the data.

Step 3 -Generating 100 random DAGs using the observed nodes.

Step 4 -Optimizing each of the 100 DAGs using Tabu search. In this step, we also extract edges present in the KG and incorporate them as a prior.

Step 5 -Calculating probability of occurrence for each possible arc in the 100 optimized DAGs.

Step 6 -Constructing the final network by selecting arc strengths above a threshold.

• Generate 100 random DAGs with nodes 𝑉 . We empirically determined, based on our previous study [21], that this number is adequate for the medium-to-large datasets 4.

• In the score function, we include edge probability contributions from the prior knowledge graph (we assign probability 0.5 for every edge in 𝛽). For each DAG, we used the stochastic algorithm Tabu [14] to find a DAG that maximizes standard Bayesian Dirichlet equivalent uniform scoring function (BDeu) [15,16].

• The previous step yields 100 optimized networks. Using these we compute the probability of each possible directed edge as its empirical frequency of occurrence among the DAGs. For example, if an edge (𝑋 , 𝑌 ) appears in 80 out of 100 optimized DAGs, we assign it an empirical probability of 0.80. We store the edge probabilities in a list.

• We threshold the edge probabilities in order to obtain the set of edges 𝐸 for 𝐺 𝑓 . Based on empirical studies, we chose a threshold of 0.85.

We chose Tabu for its robustness, simplicity (uses few parameters) and history-dependent ("memory"), although Kg2Causal is in principle compatible with any optimizing method.

Observational experiment

In the case where the dataset 𝐷 is purely observational (i.e., no interventions) from a single experiment, Kg2Causal can be implemented algorithmically as described above; we provide a pseudocode description of the "observational" formulation of Kg2Causal in Algorithm 1. Next, from each DAG in randomDAG, we learn an optimized network and store the 𝑛 networks in a list optimizedDAG. In this step, we also pass 𝛽 to createRandomDAG. Next, using the list of networks in optimizedDAG, we compute the probabilistic arc strength for each ordered pair of nodes as its empirical frequency, using procedure edgeStrength and store them as edgeProb. Finally, we use learnCausalDAG where we select the edges with weight above a predefined Threshold. Algorithm 2: Kg2Causal for mixed observational-interventional data.

Mix of Observational and Interventional experiments

With causal network learning based on a single observational dataset, it is difficult to differentiate between compatible Markov equivalent models [22]. In the simple case of three variables 𝑉 𝑖 , 𝑉 𝑗 and 𝑉 𝑘 , there are three possible causal models 𝑉 𝑖 ⇒ 𝑉 𝑗 ⇒ 𝑉 𝑘 , 𝑉 𝑖 ⇐ 𝑉 𝑗 ⇐ 𝑉 𝑘 , and 𝑉 𝑖 ⇐ 𝑉 𝑗 ⇒ 𝑉 𝑘 ; all three structures are Markov equivalent. This ambiguity can be resolved by incorporating measurements from interventional experiments, causing the Markov equivalent structures to have different likelihoods. However, in real-world settings, it is difficult to obtain such interventional measurements as compared to observational measurement [23]. Even when interventional datasets are available, learning a causal network from mixed observational and interventional data is challenging, for two reasons: (i) datasets collected from different experiments under different environmental conditions or batches are not identically distributed, in which case their underlying causal structures may differ leading to errors if network inference is applied to the combined set of measurements; and (ii) in real-world settings interventions are not "perfect" but rather "uncertain" (i.e., "imperfect" or "fat-hand"), meaning that the interventions have other unknown targets, which if ignored would likely yield spurious interactions in network discovery. To deal with such cases, based on our previous study demonstrating the effectiveness of the Learn and Vote algorithm [21,24], we extended Kg2Causal to include learning from a multi-experiment dataset using a voting-based integration method where experiment-specific causal networks are learned and combined by weighted averaging into a consensus causal network. The additional steps in Algorithm 2 are as follows:

Analysis and Results

In this section, we describe the observational datasets and ground-truth networks (Sec. 4.1) and the simulated mixed interventional-observational datasets (Sec. 4.2) that we analyzed. We present (Sec. 4.3) the results of empirical studies of network learning performance of Kg2Causal on these datasets in comparison to other types of network structure priors.

Observational datasets that were analyzed

To assess performance of Kg2Causal on biological network inference problems, we empirically analyzed five real world datasets for which published ground-truth networks were available:

Hepatic encephalopathy: This is a clinical study about a serious liver complication called hepatic encephalopathy (HE) [25] with conditions like electrolyte disorders, infections, poor spirits. It is a categorical dataset with eight nodes and ground-truth containing ten edges.

Sachs et al. T cell signaling:

This is a study on mixed observational and interventional experiments to infer causal connections between eleven protein and phospholipids in the intracellular signaling network of individual human CD4+ T-cells [26]. The dataset contains measurement of gene expressions with ground truth network containing twenty edges.

Hematopoietic Stem Cell Differentiation (HSC): This is a real-world gene regulatory network to study underlying myeloid differentiation from multipotent myeloid progenitors to megakaryocytes, erythrocytes, granulocytes and monocytes [27] in mammals [28]. The dataset contains measurement of gene expressions with ground-truth network having thirty edges.

Gonadal Sex Determination (GSD): This a real-world model which represents the gonadal differentiation circuit which monitors the transformation of the bipotential gonadal primordium (BGP) into either female or male gonads [29]. The network consists of eighteen genes and one node for the urogenital ridge. The dataset contains measurement of gene expressions with ground-truth network containing seventy nine edges [28].

Yeast cell cycle: This is a dataset derived from a network model of thirty genes participating in cell-cycle regulation of yeast [30]. The dataset was created by integrating gene expression data with transitive protein-protein interaction. The ground-truth network has 317 edges.

Mixed observational-interventional datasets

We tested Kg2Causal using Sachs et al. interventional dataset and simulated observational and interventional measurement data from synthetic networks using the bnlearn package. For observational data, we drew random samples and for interventional data, we set some target nodes in the network to fixed values in order to create mutilated networks [31] before drawing samples from them. To simulate an uncertain intervention (or "fat-hand") [32] we intervened one or more of child nodes of the intervention's target node.

Cancer: This is a synthetic network [33] on causes and consequences of lung cancer. We simulated data from one observational and one interventional experiment with equal number of samples (500) from each experiment to avoid bias. For interventional experiment we generated a mutilated network: cancer_mut with one intervention (node Smoker).

Asia: This is a synthetic network [34], about occurrence of lung disease and their epidemiological connection a prior visit to Asia. We simulated one observational and two interventional experiment from the synthetic network with equal number of samples (500) from each experiment to avoid bias. For the interventional experiments we conducted experiments to generate two mutilated networks: asia_mut1 with one intervention (node "Lung Cancer") and asia_mut2 with two intervention (at nodes "Lung Cancer" and "Tuberculosis").

Analysis of results

In this section we present results of empirical studies of network learning performance on the five observational datasets (see Sec. 4.1) and three mixed observational-interventional datasets (see Sec. 4.2), for Kg2Causal in comparison with three other types of network structure priors.

To quantify the performance, we considered presence of an edge in the ground truth network as a "true positive" and absence of an edge as a "true negative" causal arc. For the observational datasets, we used Algorithm 1 with indicated prior (KG, uniform, marginal, or Bayesian VSP) as described in Sec. 3. For mixed interventional-observational datasets, we used Algorithm 2 with the indicated prior. For each dataset, we found (Fig. 5) that using general knowledge graph as prior improves performance, by ROC, precision/recall, F1, and accuracy. Quantitatively, Kg2Causal had higher area under ROC curve (AUROC) and area under precision-recall curve (AUPR) scores than network learning with three non-KG priors tested, for the five observational (Table 1) and three mixed interventional-observational (Table 2) datasets. Moreover, the results of comparative analysis of Kg2Causal performance on mixed datasets (Table 2) show effect of pooling data from different experiments (Algorithm 1) as compared to voting (Algorithm 2) for such cases: pooling is better for small network (Cancer) (consistent with our previous findings [21]), whereas voting is better for medium-sized networks (Asia and Sachs).

Discussion and Conclusion

A limitation of this study is that due to lack of availability large ground-truth causal networks, all datasets analyzed in this work are for small to medium sized networks (8-30 nodes); due to scalability issue of score-based methods, Kg2Causal method as described here would be challenging to apply to larger networks (many hundreds to thousands of nodes and beyond), which is an area of future work. Further, we plan to explore ways to incorporate a network structure prior in constraint based algorithms (for example, PC algorithm [2]), given (in general) more favorable scalability of constraint-based algorithms and given the overwhelming preponderance of observational-only datasets that are available. We also want to evaluate alternative methods (other than the method [11] that we are using) for incorporating priors and compare them. Present work clearly demonstrates, for the case of causal network learning from small-to medium-sized biomedical or biological datasets, the importance of aggregating and leveraging structured prior knowledge in order to maximize network learning accuracy.