Investigating ADR mechanisms with Explainable AI: a feasibility study with knowledge graph mining

Abstract

Background

Adverse drug reactions (ADRs) are statistically characterized within randomized clinical trials and postmarketing pharmacovigilance, but their molecular mechanism remains unknown in most cases. This is true even for hepatic or skin toxicities, which are classically monitored during drug design. Aside from clinical trials, many elements of knowledge about drug ingredients are available in open-access knowledge graphs, such as their properties, interactions, or involvements in pathways. In addition, drug classifications that label drugs as either causative or not for several ADRs, have been established.

Methods

We propose in this paper to mine knowledge graphs for identifying biomolecular features that may enable automatically reproducing expert classifications that distinguish drugs causative or not for a given type of ADR. In an Explainable AI perspective, we explore simple classification techniques such as Decision Trees and Classification Rules because they provide human-readable models, which explain the classification itself, but may also provide elements of explanation for molecular mechanisms behind ADRs. In summary, (1) we mine a knowledge graph for features; (2) we train classifiers at distinguishing, on the basis of extracted features, drugs associated or not with two commonly monitored ADRs: drug-induced liver injuries (DILI) and severe cutaneous adverse reactions (SCAR); (3) we isolate features that are both efficient in reproducing expert classifications and interpretable by experts (i.e., Gene Ontology terms, drug targets, or pathway names); and (4) we manually evaluate in a mini-study how they may be explanatory.

Results

Extracted features reproduce with a good fidelity classifications of drugs causative or not for DILI and SCAR (Accuracy = 0.74 and 0.81, respectively). Experts fully agreed that 73% and 38% of the most discriminative features are possibly explanatory for DILI and SCAR, respectively; and partially agreed (2/3) for 90% and 77% of them.

Conclusion

Knowledge graphs provide sufficiently diverse features to enable simple and explainable models to distinguish between drugs that are causative or not for ADRs. In addition to explaining classifications, most discriminative features appear to be good candidates for investigating ADR mechanisms further.

Background

Molecular mechanisms behind harmful or beneficial effects of drugs are largely unknown. For instance, the molecular mechanism of highly prescribed drugs acetaminophen, lithium, and metformin is not completely understood. Indeed, drug development process relies mainly on randomized clinical trials and postmarketing pharmacovigilance that evaluate drug efficacy and safety, independently from any mechanistic investigation [1]. However, understanding a drug’s mechanism is fruitful: it can guide drug development, improve drug safety, and enable precision medicine, through better dosing or combination of drugs [2]. Aside from this partial ignorance, many elements of knowledge about drug ingredients are available in open-access knowledge graphs, such as their chemical and physical properties, their interactions with biomolecules such as their targets, or their involvements in biological pathways or molecular functions [3]. Knowledge graphs can be broadly defined as graphs of data with the intent to compose knowledge. Here, we consider knowledge graphs represented using Semantic Web technologies, including RDF (Resource Description Framework) and URI (Uniform Resource Identifier) [4, 5]. In such knowledge graphs, nodes represent entities, also named individuals, of a domain (e.g., acetaminophen), classes of individuals (e.g., analgesics), or literals (e.g., strings, dates, numbers). Literals are purposely discarded in this study. Nodes are connected by directed edges that are labeled with predicates (e.g., transportedBy).

We propose in this article to leverage elements of knowledge about drugs that lie in biomedical knowledge graphs to investigate ADR molecular mechanisms. To this aim, we experiment with knowledge graphs as an input to machine learning approaches (i.e., methods of Artificial Intelligence, commonly denoted AI) that are natively explainable. Indeed, Explainable AI usually refers to research on methods that provide explanatory elements to results (i.e., a classification) of sub-symbolic approaches (e.g., ensembles or Deep Neural Networks) [6]. In a broadly manner, we consider symbolic approaches that provide models that are interpretable by humans, and investigate if features of these models may be explanatory for biomolecular processes involved in ADRs. We particularly consider a knowledge graph named PGxLOD, which encompasses and connects drug, pathway, and biomolecule data [7]; and two particular types of ADRs: drug-induced liver injuries (DILI) and severe cutaneous adverse reactions (SCAR). We choose these types of ADRs first because hepatic or skin toxicities are commonly monitored during drug development for their importance in pharmacovigilance [8]. Indeed, drugs cause frequently hepatic and skin events, and those are severe enough to potentially lead to drug withdrawal in Phase IV. Second, it exists good quality expert classifications that label sets of drugs as either causative or not for these types of ADRs [9, 10]. First, our work identifies biomolecular features from our knowledge graph that enable an automatic reproduction of expert classifications. In particular, we mine the graph for neighbors of drugs, paths and path patterns (i.e., paths composed of general classes) rooted by drugs and passing by at least one entity of the following types: pathway, gene/protein, Gene Ontology (GO) term or MeSH term. Second, we isolate both predictive and interpretative features hypothesizing that, in addition to be explanatory for the classification, those may also be explanatory for ADR mechanisms. To this second aim, we consider simple, but explanatory classification techniques, i.e., Decision Tree and propositional rule learner over extracted features, because they provide human-readable models in the form of rules. Finally, we ask three human experts if they consider isolated features as possibly explanatory for ADRs.

A first family of related works can be described as explanatory, where known Drug-ADR associations, such as those found in SIDER, are used to highlight molecular mechanisms that may be impacted in ADRs. A second family of works is predictive, where data about molecular mechanisms (e.g., GO molecular processes or KEGG pathways) are associated with drugs and used as features to predict ADRs [11]. Boland et al. [12] survey existing works for both predicting ADRs and elucidating their mechanisms; they interestingly list data and knowledge resources that may support these efforts.

In the explanatory family, Lee et al. [13] associate SIDER side effects and GO biological processes through drugs, by the combination of a Drug-Side effect and a Drug-Biological process networks. Highlighted relationships are obtained using statistical approaches (i.e., enrichment and t-score) and evaluated in regard to co-occurrences in PubMed abstracts. Wallach et al. link drugs to proteins through molecular docking, then to pathways through the KEGG database, and use logistic regression and feature selection approaches to select pathways most probably impacted in side effects [14]. Bresso et al. [15] group frequently associated drug reactions and propose elements of explanations of their grouping using Inductive Logic Programming. Also, Chen et al. [16] propose a computational algorithm to infer Protein-ADR relationships from a network of protein-protein interactions, ADR-ADR similarities and known protein-ADR relations.

In the predictive family, Bean et al. [17] build a network of drugs, targets, indications and ADRs to select features that are good predictors for ADRs in a logistic regression setting. PhLeGrA is an analytic method implementing Hidden Conditional Random Fields to allow the calculation of the probability of drug reactions given a input drug and a knowledge graph of drugs, proteins, pathways and phenotypes [18]. Similarly, Muñoz et al. [19] propose a specific way to extract features from knowledge graphs for ADR prediction. They show that several multi-label learning models perform well for this task. Our work is similar to some extent, however it uses simpler but explanatory classifiers, and goes a step further by identifying features, subsequently proposed as explanatory elements for ADRs. Indeed, we hypothesize that within the large set of considered features, those that are both highly predictive and associated with a good level of interpretability may suggest to experts plausible elements of explanation. In Dalleau et al. [20], knowledge graph mining serves a completion perspective and aims at inferring links between drugs and genes. All PhLeGrA, Muñoz et al., Dalleau et al., and the present work illustrate the interest of aggregating several LOD (linked open data) sets for knowledge discovery and data mining tasks, as discussed by Ristoski and Paulheim [21].

Shi and Weninger [22] use a similar approach to ours, but from a fact checking perspective. Indeed, for each relation type p, a set $D^k_{(o_u,o_v)}$ of discriminative paths is learned. This set contains anchored predicate paths $o_u \xrightarrow []{r_1} \xrightarrow []{r_2} \dots \xrightarrow []{r_{k}} o_v$ of length k that describe a statement $o_u \xrightarrow []{p} o_v$, where $o_u$ and $o_v$ are respectively the ontology classes associated with nodes u and v. To check whether a triple $s \xrightarrow []{p} t$ is true, they use the learned set of discriminative paths for the relation $D^k_{(o_s,o_t)}$. In such sets, only paths with the most discriminative power are kept. Similarly to our approach, they use ontology class generalization but apply it only to start and end nodes s and t of the fact to be checked. This differs from our approach as we apply generalization on each intermediate node (see Methods Section for details). Additionally, their path modelling allows reverse traversal, i.e., $\xrightarrow []{p^{-1}}$, and constrains both source and target nodes, while we only constrain source nodes. Previous works also use knowledge graph mining to provide explanations. Those include Explain-a-LOD that enriches statistical data sets with features from DBpedia. It uses correlation between attributes and rule learning to provide hypothesis explaining statistics [23]. Explain-a-LOD relies on FeGeLOD to extract features from the DBPedia knowledge graph [24]. In particular, FeGeLOD extracts two types of features similar to ours: paths of size 1 ($\xrightarrow []{r} e$) starting at the entities of interest; paths of size 1 ($\xrightarrow []{r} t$) where the original entity (e) is replaced with ontology classes (t) it instantiates. In the same vein, KGPTree [25] extracts paths of the form root $\rightarrow$ predicate $\rightarrow$ entity $\rightarrow$ predicate ...$\rightarrow$ entity, while allowing generalizations on both predicates and entities, whereas we offer generalization on entities only. However, they only allow a generalization to a unique and broad type denoted with a wild card (*), while we allow a granular generalization following ontology class hierarchies. FeGeLOD and KGPTree extract only paths and path patterns, whereas one may want to extract other common structures such as subtrees. Mustard Python library offers such functionalities applying Graph Kernels on RDF graphs, plus additional facilities such as detecting hubs or low frequency patterns [26, 27]. However, it does not allow generalization operations.

The contribution of our work is twofold: first, we show that knowledge graphs provide sufficiently diverse features to enable simple and explainable models to distinguish between drugs that are causative or not, for two types of ADRs commonly monitored; second, we manually evaluate in a mini-study that in this setting, predictive features constitute good candidates for investigating ADR mechanisms further. The following sections present materials and methods, obtained results and a discussion about our methodological choices and results.

Materials and methods

Data sources

PGxLOD

PGxLOD is a linked open data (LOD) knowledge graph built for pharmacogenomics (PGx) and encoded in RDF [7]. It aggregates data mainly about drugs, genetic factors, phenotypes and their interactions from six data sources: PharmGKB, ClinVar, DrugBank, SIDER, DisGeNET and CTD; but also includes references to Gene Atlas, UniProt, GOA and KEGG. In particular, it includes pharmacogenomic relationships, i.e., n-ary relations that represent how a genomic factor may impact a drug response phenotype. These relations are compiled from PharmGKB and the literature. We used PGxLOD version 4 that encompasses 88,132,097 triples. Table 1 presents its main statistics. PGxLOD is available at https://pgxlod.loria.fr.

Table 1 Types and numbers of entities available in the PGxLOD knowledge graph

In our study, PGxLOD enables to associate phenotypic and molecular features to drugs, by the exploration of their neighborhood in the graph.

Drug expert classifications and their preprocessing

We experiment with two expert classifications of sets of drugs labeled as either causative or not for ADRs. The first concerns drugs causative for drug induced liver injury (DILI), and the second drugs causative for severe cutaneous adverse reactions (SCAR).

DILI classification We built our DILI classification from DILIRank, a list of 1036 FDA-approved drugs classified by their risk of causing DILI [9]. DILIRank distinguishes between four classes listed in Table 2. This classification was obtained by (1) the curation of information gathered from FDA-approved drug labels, setting an initial list of 287 drugs; (2) a semi-automatic approach that completes the list up to 1036 drugs, by combining information from hepatotoxicity studies and the literature. We sub-sampled from DILIRank 370 drugs (146 DILI$\oplus$, and 224 DILI$\ominus$) that fulfill criteria required for our subsequent analysis: being either in the Most—or No—DILI concern classes, being associated with a SMILES (simplified molecular-input line-entry system) description, and being mapped to PGxLOD. The latter mapping was obtained with PubChemIDs, which are available both in DILIRank and PGxLOD (coming from PharmGKB, DrugBank, and/or KEGG). Drugs satisfying these criteria and classified as Most-DILI concern constitute the DILI$\oplus$ subset, and those classified as No-DILI concern constitute the DILI$\ominus$.

Table 2 Classes and size of the original DILIRank expert classification

SCAR classification Our SCAR classification relies on a manually built classification called “Drug notoriety list”, shared by members of the RegiSCAR project. This list was originally assembled for the evaluation of the ALDEN algorithm, which assesses the chance for a drug to cause Stevens-Johnson Syndrome and Toxic Epidermal Necrolysis, a specific type of SCAR [10, 28]. This classification lists 874 drugs and assigns them to five classes representing various levels of association with SCAR. These classes are listed in Table 3. We sub-sampled from the RegiSCAR drug notoriety list 392 drugs (102 SCAR$\oplus$, and 290 SCAR$\ominus$) fulfilling two criteria: being mapped to PGxLOD and having a SMILES description. The mapping starts with 874 drug names, which is the only description available in the RegiSCAR list. Drug names are matched with lists of synonyms associated with drugs in PharmGKB and DrugBank. Drugs satisfying these criteria and classified as Very probable, Probable or Possible constitute the SCAR$\oplus$ subset, and those classified as Unlikely or Very unlikely, the SCAR$\ominus$.

Table 3 Classes and size of the original RegiSCAR drug notoriety list

Methods

An overview of the steps of the proposed method is provided in Fig. 1.

Fig. 1

Method overview. Main steps are in orange, data in blue, parameters in green and models in pink

Graph-based feature construction

Graph canonicalization In knowledge graphs, several nodes may co-exist, while representing the same entity. For example, a unique drug in PGxLOD can be represented by two nodes: one from PharmGKB and another from DrugBank. Each of these two nodes have their own connections to other nodes in the knowledge graph. Therefore, the union of their edges constitutes all the available knowledge about the drug they both represent. As they represent the same drug, they are connected through an owl:sameAs edge. In order to avoid traversing such edges, we canonicalize the knowledge graph, that is to say, nodes linked by owl:sameAs edges are grouped into a unique node as a pre-processing step, before building graph-based features. Such a procedure corresponds, in graph theory, to the contraction of owl:sameAs edges.

Paths, path patterns, and neighbors as drug features

From an arbitrary set of drugs D (such as DILI$\oplus$ $\cup$ DILI$\ominus$) and a canonicalized knowledge graph (such as PGxLOD), a set of graph-based features is built on D. We distinguish three kinds of features: paths, path patterns and neighbor nodes. First, we build paths rooted by drugs from D. The neighborhood of each drug $d \in D$ is explored in the graph with a max distance of k, generating sequences of predicates and nodes of max length k. Accordingly, if a path $d \xrightarrow {p_1} n_1 \xrightarrow {p_2} n_2$ is found in the knowledge graph (with $k=2$), the drug d is associated with the feature $\xrightarrow {p_1} n_1 \xrightarrow {p_2} n_2$ in the output matrix.

Second, we build path patterns that generalize paths by considering ontology classes instantiated by nodes. The aim of path patterns is to offer more general descriptions, which have more chances to be shared by several entities. For example, if $n_1$ instantiates $C_1$ and $n_2$ instantiates $C_2$, we add the following path patterns: $\xrightarrow {p_1} C_1 \xrightarrow {p_2} n_2 ,\ \xrightarrow {p_1} n_1 \xrightarrow {p_2} C_2,\ \xrightarrow {p_1} C_1 \xrightarrow {p_2} C_2 ,\ \xrightarrow {p_1} \top \xrightarrow {p_2} n_2 ,\ \xrightarrow {p_1} n_1 \xrightarrow {p_2} \top ,$ and $\xrightarrow {p_1} C_1 \xrightarrow {p_2} \top , \xrightarrow {p_1} \top \xrightarrow {p_2} C_2 \text{, and} \xrightarrow {p_1} \top \xrightarrow {p_2} \top$. It is noteworthy that, to allow a high level of generalization, we always consider the top-level class $\top$ in the generalization procedure. To leverage hierarchies organizing ontology classes, a node n is only generalized by ontology classes at a distance of at most t from itself, following instantiation and subsumption edges.

Third, we list neighboring nodes, i.e., any node that can be reached from $d \in D$ within a distance of max k. We do not keep track of the distance between d and its neighbors: if a node n can be reached from a drug $d_1$ after 2 hops and from a drug $d_2$ after 3 hops, these two drugs will be associated with the neighbor n in the output matrix. Interestingly, a neighbor n can be represented by the general path pattern $(\xrightarrow {*}*)^{\{0,k-1\}}\xrightarrow {*}n$.

Features can potentially be noisy and numerous, leading to a combinatorial explosion of their number. That is why, we add several constraints to only keep the most interesting features. First, we only keep the most specific paths and path patterns, considering that a node is more specific than the classes it instantiates and a class is more specific than its superclasses. A path or path pattern $P_1$ is more specific than another path or path pattern $P_2$ if each node/class in $P_1$ is more specific than the node/class at the same position in $P_2$. Thus, if $\xrightarrow {p_1} n_1 \xrightarrow {p_2} n_2$ is associated with the exact same drugs as $\xrightarrow {p_1} C_1 \xrightarrow {p_2} C_2$, then the path pattern is removed and only the path constitutes a feature. Similarly, if $\xrightarrow {p_1} C_1 \xrightarrow {p_2} C_2$ is associated with the exact same drugs as $\xrightarrow {p_1} \top \xrightarrow {p_2} \top$, then we only keep the first one. Second, the exploration is stopped at nodes whose degree is greater than a parameter deg. Such nodes are usually called hubs [27], and expanding a path ending at a hub may generate a large number of new paths associated with the same drugs. Third, we consider minimal and maximal supports of features (denoted $s_{\text {min}}$ and $s_{\text {max}}$) as additional filters. The support of a feature consists in the number of drugs associated with a feature, i.e., neighbors of a same entity or rooting a path or a path pattern. Only features associated with more than $s_{\text {min}}$ drugs and less than $s_{\text {max}}$ drugs are output in the final matrix. Fourth, two blacklists avoid traversing noisy or unwanted edges. Edges are not traversed if their predicate is blacklisted (in $b_{\text {predicates}}$) or if the adjacent node instantiates (directly or indirectly) a blacklisted class (in $b_{\text {exp-types}}$). For example, we blacklist predicates of PROV-O that are used to describe provenance metadata. By blacklisting rdf:type in $b_{\text {predicates}}$, we make sure the exploration is performed through entities and not classes of PGxLOD. Relations between classes are only considered when generalizing paths. Aside from noisy features and combinatorial explosion, we use these blacklists to prevent considering features that “obviously” carry the signal we are trying to predict. For example, we blacklist all “side-effects” links from SIDER, which may directly link drugs in D with the side effect we aim at predicting. We also blacklist all classes from MeSH related to SCAR or DILI to avoid taking into account nodes instantiating them. A third blacklist ($b_{\text {gen-types}}$) avoids generalizing nodes in paths by blacklisted classes. This is particularly useful to withdraw general classes (e.g., Drug) that increase the number of generated path patterns while not adding useful information.

Besides previous topological constraints, we also perform a domain-driven filtering configured by parameter m. Indeed, in our objective of explaining ADRs, experts highlighted that features that may be explanatory to them mention pathways, genes, GO, and MeSH terms. For this reason, we propose three post-processing atomic filters, only keeping neighbors, paths or path patterns containing at least a pathway ($m={\texttt {p}}$), a gene or a GO class ($m={\texttt {g}}$), or a MeSH class ($m={\texttt {m}}$). Such atomic filters can be combined to form disjunctive filters. For example, the $m={\texttt {pg}}$ filter keeps neighbors, paths or path patterns containing at least a pathway or a gene or a GO class.

Table 4 summarizes the parameters that limit the number of features in the output matrix. Additional details about our method of feature construction are provided in [29].

Table 4 Parameters used to limit the number of features

Cross-validation strategy

Once drug features are extracted from the knowledge graph, they are given to a machine learning algorithm that learns a model, which mimics the expert classification. To quantitatively evaluate our approach, we adopt a 10-fold cross-validation strategy, meaning that we repeated the following steps of our learning pipeline (i.e., feature selection, and training) 10 times, holding out each time one tenth of our labeled data for testing. The split in 10 folds is performed once, randomly.

Feature selection with boosting

The concept of boosting relies on the sequential learning of classifiers that successively focus on getting correct the examples that were wrongly classified by previous classifiers, with a weight system: correctly classified examples loose weight, whereas falsely classified ones gain weight [30]. The final classification is built by a weighting of the results obtained from the various classifiers. Ensembles of decision trees can be used in a boosting strategy for the estimation of feature importance in order to select the most important features for a subsequent learning model [31]. This approach is frequently named “wrapper-based feature selection”. We adopt this approach, using the AdaBoost algorithm, with up to 10 decisions trees learned from the whole train set [32]. In case of perfect fit, the learning procedure is stopped early. For the Decision Tree algorithm we used CART from scikit-learn, Gini impurity, all features considered at each split and a $minimal\ number\ of\ examples\ per\ leaf = 5$. Every train example is used, but classes are artificially balanced: examples are associated with a different weight depending on whether they belong to the $\oplus$ or $\ominus$ class. Features that appear in at least one of the decision trees are selected for the subsequent step of our learning pipeline.

Performance and robustness evaluation

For evaluating our capacity of distinguishing between drugs associated or not with an ADR, we train a last decision tree, using selected features only. Note that because the selection step is repeated at each iteration of the cross validation strategy, selected features may vary from one iteration to another. Algorithm and parameters are the same as those of the selection step, i.e., CART from scikit-learn with Gini, all features considered at each split, minimal number of examples per leaf = 5, and weighted instances for class balancing.

Performance results are reported in term of Precision, Recall, F1-score (reference class $\oplus$), accuracy and AUC-ROC. Metrics are averaged over the 10 iterations of the cross validation.

A robustness evaluation is also performed to assess the impact of the train set (i.e., expert classifications) on the final classifier. To this aim, first, we reproduced the experiment, but with a shuffled class assignment in the train set. This leads to a train, denoted by $\cdot ^{shuffled}$ (where $\cdot$ is either DILI or SCAR) where drugs are associated randomly to either the class $\oplus$ or $\ominus$. This first sanity evaluation mainly checks the presence of a nonrandom signal in the train. Second, we replace the set of drugs from $\ominus$ by a set of drugs randomly picked out of our knowledge graph. We repeat this random draw five times for each expert classification resulting in 10 train sets denoted $\cdot ^{random \ominus _{i}}$ (where $\cdot$ is either DILI or SCAR and i is an index taking values from 1 to 5). In each case the draw is made from nodes of the knowledge graph that instantiate (directly or indirectly) pgxo:Drug; are identified with a URI from PharmGKB or DrugBank namespaces; are linked by a x-pubchem predicate to a PubChem URI (in order to have a SMILES associated with the drug); and are not drugs of the original DILI or SCAR expert classifications. In the case of DILI, this is a draw of 224 nodes out of 5893 in the canonical graph. In the case of SCAR, this is a draw of 290 nodes out of 5921. This second evaluation checks the impact of the selection of negative examples on classifier performances. In addition, we count the number of features that are present in 5, 4, 3 or 2 of the five $\cdot ^{random \ominus _{i}}$ experiments.

Qualitative evaluation by human experts

To go beyond performance evaluation, we produce a classification model made of rules using the RIPPER algorithm, and its Java implementation named JRip [33]. JRip has the advantage over decision trees to provide classification rules that are more concise and by consequence easier to interpret for humans. JRip actually implements a propositional rule learner, produces relatively non redundant rules in comparison to rules that could be learned following branches of a decision tree. However, JRip and CART decision tree usually perform very similarly, since they implement similar pruning strategies and stopping criteria. To evaluate this assessment, we performed a 10-fold cross validation of the JRip approach and compared performances with CART. JRip generates rules of the form:

$$\begin{aligned} \displaystyle \left( \bigwedge _{\forall a \in A} a \right) \wedge \left( \bigwedge _{\forall b \in B} \lnot b \right) \Rightarrow c \end{aligned}$$

where A is a set of attribute-value pairs affirmed, B is a set of attribute-value pairs negated, and c is the minority class of the classification problem. Accordingly, $c=\oplus$ in our study. Note that A or B can be an empty set, but not both at the same time. In other words, JRip rules consist in the conjunction of affirmed and negated features.

As for CART Trees, we start with an initial step of feature selection with AdaBoost, with the same parameters but this time considering all examples of the train set. Following, JRip rules are built considering also all examples of the train set, and features that appear in at least one of the trees built by AdaBoost. To be consistent, with the CART Tree setting, we set to 5 the minimal number of instances per rule, which can be compared to the minimal number of examples per leaf.

JRip rules are post-processed to facilitate their interpretation by our experts. First, among features that are path patterns, we discard those only involving generic classes such as Resource or Drug. Indeed, such path patterns turned out to be impossible to interpret. Second, features are translated into a readable format, by resolving ids with associated labels (e.g., drugbank:BE0003543 is resolved as “Cytochrome P450 1A1”) and by interpreting and rewriting paths and path patterns in an understandable form (e.g., $\xrightarrow {\texttt {drugbank\_vocabulary:enzyme}}\texttt { \small geneatlas\_vocabulary:Resource}$ is turned in $\xrightarrow {\text {interactsWith}} \text{Enzyme}$). Because of the limited number of features in rules, this translation is made manually, on the basis of descriptions of predicates, classes, and entities found in their original knowledge graph or database.

We asked three experts in pharmacy and pharmacology to review independently each attribute (i.e., feature) of the rules, to evaluate if they may be explanatory for ADRs. Each expert has to answer a voluntarily simple three-way question: “according to your own knowledge or the state of the art, do you think that the feature is explanatory for DILI?” (SCAR, respectively). Possible answers are “yes”, “maybe (possible, but not obvious)” and “no (probably not explanatory)”. We allow expert to check the literature or any state-of-the-art resource, but up to 15 min, since we consider that more time causes to fall in the “no (probably not explanatory)” option. To guide their decision on each feature, experts are provided with two Web links: one pointing to the list of drugs from the $\oplus$ train set that supports the feature; one to the page of the main entity mentioned in the feature (i.e., the neighbor node, or the final node of the path or path pattern) in an expert database: DrugBank, ChEBI, KEGG, QuickGO, or MeSH browser, depending on the namespace of the node. After expert reviews, answers were normalized, under their supervision, to guarantee all experts interpret the negation of features the same way. For each feature we check if at least one, two or three of the experts think it is or may be explanatory, if the three think it is explanatory and if the three think it is not. In addition, we compute Cohen’s kappa coefficient to evaluate the average agreement between experts with two different settings: considering the three different answers as distinct, or considering answers “yes” and “maybe” as a unique positive answer.

Results

Graph-based feature construction

We experimented graph exploration with combinations of the following parameters values $k \in \{1,2,3,4\}, t \in \{1,2,3\}, deg = 500, undirected = {\texttt {false}}, s_\text {min} = 5, s_\text {max} = +\infty$ and $m \in \{{\texttt {p}}, {\texttt {g}}, {\texttt {m}}, {\texttt {pg}}, {\texttt {pgm}}\}$. $k = 4$ was only tested with $t = 1$ because of memory issues caused by the high number of generated features with greater values of t. However, we report only the best results, which were obtained with $k = 3, t = 3, m = {\texttt {pgm}}$ for DILI and $m = {\texttt {pg}}$ for SCAR. These values of m enable to conserve only features that includes an entity that is either a pathway, a gene, or a GO term (for ${\texttt {pg}}$); or a pathway, a gene, a GO, or MeSH term (for ${\texttt {pgm}}$).

The construction of graph-based features is obviously limited by the amount of available memory. We enforce $s_\text {min}$ to be set to allow the construction of paths and path patterns, in order to avoid combinatorial explosion. Accordingly their number is reported only once we reduced the number of all possible combinations, which we were not able to compute. We used a server with a Xeon E5-2680 v4@2.40GHz CPU, 28 cores/56 threads and 768GB of memory. As an illustration, we obtained the features associated with the DILI expert classification under $k = 3, t = 3$ in approximately 1 h using 62 GB of RAM, and under $k= 4, t= 1$ in 4 days using 380 GB of RAM.

To provide with an idea of the size of the considered neighborhood with regards to all reachable nodes, Table 5 reports statistics about numbers of neighbors, paths and path patterns reachable with different level of filtering. In particular, we report sizes of the full neighborhood of drugs, and of 3 levels of filtering. The first level of filtering consists in prohibiting the expansion of the neighborhood through nodes with a degree higher than 500 ($\textit{deg}=500$). In the full neighborhood and first level of filtering, k and t are not constrained since no path or path pattern is computed, but we report max k and t reached in the neighborhood. We note that k is surprisingly lower in the larger neighborhood, i.e., 19 and 20 versus 23 and 23 with the first level of filtering. This can be explained by the fact that with “hubs” (nodes with $\textit{deg} > 500$), the full neighborhood can be reached through smaller paths. We also observe that because of this first filtering, certain nodes are not accessible anymore (when every possible path to them pass through a hub), which results in a smaller number of neighbors. The second level of filtering comes on top of the first, and constrains neighbors, paths, and path patterns to have a minimal support set of 5 ($s_\text {min} = 5$), a max length of 3 ($k=3$), and a max depth 3 for generalization of paths into patterns ($t=3$). The filtering level 3 comes on top of the second, and constrains neighbors, paths, and path patterns to contain a pathway, a gene, a GO term, or a MeSH term for DILI ($m = \texttt {pgm}$), or to contain a pathway, a gene, or a GO term for SCAR ($m = \texttt {pg}$). The filtering level 3 is the one used in the following experiments, because it is computable in our setting, while providing the best performances in our set of experiments.

Table 5 Numbers of drug features extractable from the knowledge graph, with different levels of filtering

Quantitative and robustness evaluation

Performances of our CART decision trees to distinguish between drugs associated or not with ADRs are reported in Table 6. With both types of ADRs, we obtained accuracy and AUC higher than 0.70, illustrating the fact that learned classifiers reproduced a large part of expert classifications, on the basis of features from the knowledge graph.

Robustness of classifiers is illustrated by the results provided in Table 7. $DILI^{shuffled}$ and $SCAR^{shuffled}$ are associated with AUC of 0.49 and 0.51, respectively. This illustrates that a random assignment of class labels in the train set, leads to a classifier that randomly assign labels to test examples. This is expected, but illustrates that expert classifications encompass a signal that our classifiers learn and reproduce, to some extent.

Classifiers trained with a random pick of negative examples ($\cdot ^{random \ominus _{i}}$) instead of negatives picked by experts are significantly better for the three performance metrics (t test, $t > 36, p <1.7 \times 10^{-6}$). This reveals it is harder for our classifier to discriminate between positives and negatives of expert classifications, than it is between positives and randomly picked drugs. This lets us assume that negatives from expert classifications are somehow similar to positives (they may share some properties) and harder to distinguish for the classifier, even if not associated with the studied ADR. When comparing features used in the five random pick of $DILI^{random \ominus _{i}}$, we observed that, out of a mean of 122 features (sd=8), 6, 12, 27 and 99 were common to respectively 5, 4, 3 and 2 picks. With $SCAR^{random \ominus _{i}}$, out of 108 features (sd=8), 2, 6, 23 and 81 were common to respectively 5, 4, 3 and 2 picks.

Table 6 Quantitative evaluation of our classifiers of drugs associated with ADRs or not (DILI or SCAR)

Table 7 Robustness evaluation of our classifiers

Expert evaluation

JRip produced 6 and 5 rules for DILI and SCAR, respectively. The translation of these rules is available in Additional file 1: Tables S1 and S2. After removing uninformative features, we obtained 11 and 13 distinct features, respectively. Those are provided in Additional file 1: Tables S3 and S4. These features are those reviewed by our three experts in pharmacology (CB, CNC, and NP). Quantitative performances of the JRip algorithm are presented in Table 6 for comparison with CART. On both datasets (DILI and SCAR) differences in performances (Precision, Recall, Accuracy, F-measure) are not statistically significant (t test, $p<0.05$).

The ratio (and number) of features for which experts reached an agreement, or for which at least one, two or three of the experts think they are or may be explanatory (answers “yes” or “maybe”) are provided in Table 8. Detailed results of the manual evaluation are provided in Additional file 2. We observe that no feature generated by JRip is considered as unexplanatory by all three experts. In other words, every feature is thought as possibly explanatory by at least one expert. The ratio of features having a full agreement between experts on the possibility of being explanatory is reduced compared to those having a partial agreement but stays relatively high (0.73 features) for DILI and moderate (0.38 features) for SCAR. Full agreement for features being explanatory (all three experts answer “yes”) remains minor, but exists. Kappa’s Cohen agreement score is $\kappa _{n=3} = 0.26$ when considering answers “yes”, “maybe”, and “no” independently, but reaches $\kappa _{n=2} = 0.70$ when the problem is reduced to two classes by aggregating “yes” and “maybe” answers. Note that Table 8 reports in its fifth column the ratio of features that reach full agreement for our three experts when “yes” and “maybe” answers are aggregated. Additional file 2 contains results of the manual evaluation of the features.

Examples of features and elements of interpretation

Three features reach an agreement for being explanatory (i.e., three answers “yes” per feature). Those can be interpreted as elements that are well known for being explanatory, or at least associated, with DILI or SCAR mechanisms. As a first illustration, reached an agreement for DILI. This is explained by the fact that endoplasmic reticulum is known, in particular in liver tissues, to host primarily cytochrome P450 enzymes, well known for being involved in drug metabolism [34]. As a second illustration, Cytochrome P450 2B6 reached an agreement for SCAR, whereas genomic variations in the gene coding for this enzyme are known for being associated with SCAR [35]. One might consider fairly that this feature does not bring new explanatory elements, although it can be considered a minimum that our method highlights well established explanatory elements.

All other features did not reach an agreement, or reach one for “maybe”. Each of those is interesting to explore for further interpretation, but for the sack of briefness, we will only discuss two of them here. First, reaches an agreement for potentially being explanatory (i.e., three answers “maybe”) for DILI. This path pattern is relatively complex to interpret since it is long ($k=3$) and negated. Experts searched for literature reporting associations between DILI and Calcium signaling pathway. They found that a relatively old literature (old is seen as lacking confirmation by some experts) were reporting such an association [36]. A more recent bioinformatics article by Chen et al., also reported such an association, but with a finer grain of information, since they report an association with hepatomegaly (a secondary example of DILI), and a negative association with hepatitis (a primary example of DILI) [37]. Accordingly, negative results from this study are consistent with our finding of this latter negated feature. We note that Chen et al. study is computational, as is ours, and that we may also be impacted by similar bias. Second, for SCAR, obtained very diverse opinions, with one “no”, one “maybe” and one “yes”. Even if this disagreement could be perceived as inconclusive, it may also point to a promising candidate for explanation. In this very one case, Electron transport is known for being perturbed in mitochondrion in many types of ADRs, including hepatotoxicity [38]. However, we did not find any study reporting a link with skin toxicities.

Table 8 Ratio of features that either reach a full agreement for being unexplanatory, explanatory or are considered as possibly explanatory to various extents

Discussion

In our work, simple, but explainable, classifiers (CART Decision Trees and JRip) were preferred to more advanced machine learning methods, even if we are convinced that methods based on deep neural networks, such as graph embedding with Graph Convolutional Networks (GCN), should lead to better performances than those obtained [39, 40]. However, acquiring explanatory elements about decisions made by such models necessitates an additional step of neural network analysis, such as saliency maps [41], which provides information such as the layer or neurons activated by some instances. We consider this information of high interest for data scientists, but such methods require high level interpretation before being understandable by typical domain experts, unfamiliar with neural networks [42]. Consequently, such direction did not seem mature enough to reach our objectives. For instance, heatmaps that are to some extent explanatory for image classification, are still hard to transpose to knowledge graphs [43]. However, it would be of interest to evaluate performances of a GCN on the classification task to measure the gap caused by our choice of simple classifiers. We also hope that our work will motivate studies on explainable subsymbolic approaches.

Our approach is reproducible for other applications. For instance, the same rational could be applied for the investigation of the mechanism of drug beneficial effects. This would necessitate to change our expert classifications for lists of drugs with a same indication ($\oplus$) and drugs without effect for this indication ($\ominus$), which could be obtained from SIDER or DailyMed.

An objective of our work is to illustrate various advantages of mining knowledge graphs, and particularly Semantic Web ones. First, they provide human-readable features, that may subsequently be interpreted by experts. Second, predictive features may come from various connected data sets and be jointly used in a single rule, which would have not been found if data sets were considered isolated. In addition, using Semantic Web standards eases the addition to our graph of novel data, following other owl:sameAs links. Third, Semantic Web knowledge graphs are associated with a formal semantics we benefit from at two steps: at the initial canonicalization, and at the generalization of path patterns. In this regards, one may ask if we could benefit from additional reasoning mechanisms, such as generalization over predicates. In our very specific case, predicates are not associated with any hierarchy, so it would not have changed our results, but from a general point of view, path patterns would benefit from this mechanism. Similarly, we could think of a canonicalization, not only with owl:sameAs links, but also following properties carrying similar semantics (e.g., skos:exactMatch) or by applying matching approaches such as PARIS [44].

When testing with values between 1 and 3, we observed that larger t and k are associated with better performances. To achieve this, we adopted a rational approach for scaling the mining of RDF knowledge graphs, which is presented in [29]. This approach reaches its limits with $k>3$ and $t>3$, but we think that additional optimization in the graph mining algorithm is still possible and would enable going further.

We used only binary features as they are easy to consider as explanations. Other strategies (e.g., counting, relative counting [45]) could also have been considered while maybe hindering the descriptive power of such features. To maximize the descriptive power of candidate features and avoid redundancy, one could use specific metrics (e.g., approaches relying on hierarchies [46] and/or extent of classes [47] of ontologies). Such metrics could also be considered within the decision tree algorithm, to propose to the algorithm additional features (more or less aggregated according to generalization) that may be associated with best split with regards to the Gini index (or others). This would lead to the consideration of formal knowledge directly in the mining algorithm [21]. Also, we proposed, with our parameter m, an hardcoded way of selecting features of interest. However, one may think of an interactive selection by user, following the possibilities offered by an ontology.

A usual difficulty in human annotation or human evaluation is the normalization of expert answers. Despite a 1-h training about the task, the interpretation of negated features has been heterogeneous among experts. One pitfall was to think that if the affirmation of a feature is true, its negation is wrong. This is misleading because it is possible that a feature is explanatory for some examples and that its negation is also explanatory for other examples or in another context. To ensure normalization of negated features, we considered a feature as explanatory if its affirmation or its negation is explanatory to the expert. This change has been considered in a normalisation batch of reviews done in cooperation with experts.

Our review by human experts evaluates how many features highlighted by our approach are relevant (similarly to what Precision measures), but does not evaluate how many relevant features we may miss (similarly to what Recall measures). It would be of interest to ask experts what are features such as pathways, drug targets, cellular functions that are known to be associated with DILI and SCAR ADRs to enable a final comparison. However, establishing an exhaustive list from the state of the art would be complex and time consuming for experts. Text mining approaches could be of interest to guide them in this matter.

Conclusion

We illustrate in this work that life science knowledge graphs provide sufficiently diverse features to enable simple and explainable models to distinguish between drugs that are causative, or not, for two severe ADRs. These features take the form of paths, path patterns or simple neighboring nodes from the graph, which have the advantage, when adequately selected, of being human-readable and interpretable by experts. We quantify through a small-sized human evaluation that such features are not only discriminative, thus predictive for the classification, but also appear to be good candidates for providing explanatory elements of ADR mechanisms. In conclusion, this work illustrates that simple models, fed with diverse and explicit knowledge sources such as those connected in the form of linked open data constitute an alternative to complex models, efficient but hard to interpret. A natural perspective is to combine such rich sources of background knowledge with models that are both highly performing (such as GCN) and interpretable.