Causal reasoning over knowledge graphs leveraging drug-perturbed and disease-specific transcriptomic signatures for drug discovery.

Network-based approaches are becoming increasingly popular for drug discovery as they provide a systems-level overview of the mechanisms underlying disease pathophysiology. They have demonstrated significant early promise over other methods of biological data representation, such as in target discovery, side effect prediction and drug repurposing. In parallel, an explosion of -omics data for the deep characterization of biological systems routinely uncovers molecular signatures of disease for similar applications. Here, we present RPath, a novel algorithm that prioritizes drugs for a given disease by reasoning over causal paths in a knowledge graph (KG), guided by both drug-perturbed as well as disease-specific transcriptomic signatures. First, our approach identifies the causal paths that connect a drug to a particular disease. Next, it reasons over these paths to identify those that correlate with the transcriptional signatures observed in a drug-perturbation experiment, and anti-correlate to signatures observed in the disease of interest. The paths which match this signature profile are then proposed to represent the mechanism of action of the drug. We demonstrate how RPath consistently prioritizes clinically investigated drug-disease pairs on multiple datasets and KGs, achieving better performance over other similar methodologies. Furthermore, we present two case studies showing how one can deconvolute the predictions made by RPath as well as predict novel targets.

This is a PLOS Computational Biology Methods paper.

Introduction

The representation of biomolecular interactions occurring within cells is often intuitively organized in the form of biological networks. These networks can be used to inherently model biological processes through the use of nodes denoting biological entities and edges representing their relationships. While homogeneous networks, such as protein-protein interaction networks, can represent relationships between a single entity type, knowledge graphs (KGs) can incorporate a broad range of biological scales, from the genetic and molecular level (e.g., proteins, drugs, and biochemicals), to biological concepts (e.g., phenotypes and diseases). These KGs can then be utilized for several applications in drug discovery, such as providing insights into molecular mechanisms and therapeutic targets [1-2], side effect prediction in the early stages of drug development [3], target prioritization [4], and drug repositioning [5].

Given the flexibility of KGs, multiple heterogeneous relation types can be modeled to represent biological processes that are governed by interactions occurring between component entities [6]. Even though a variety of relation types (e.g., literature co-occurrence, associations, etc.) can be leveraged by network-topology algorithms for various applications, causal relations are particularly useful as they can be used to infer the effect of any given node on another by reasoning over the KG [7]. Nonetheless, not all interactions included in a given KG are necessarily biologically relevant as they may be context-specific, such as to a particular cell type, tissue or disease. Furthermore, as the complete human interactome remains unknown, KGs modeling PPIs are also incomplete and the interactions which are modeled tend to be biased towards well-studied proteins and their relationships [8]. One approach to address these challenges is to jointly leverage prior knowledge in KGs with data-driven -omics experiments [9-13].

Experimental datasets have been widely employed by recent drug repurposing approaches to identify drug candidates for a given disease using the anti-correlation in biological processes or pathways at the transcriptomic- or proteomic- level between drugs and diseases as a proxy [13-16]) (see [17] for a recent review and for a survey of such methods). While these approaches use prior knowledge in the form of pathways (gene sets), this concept has yet to be applied on KGs for drug discovery. However, by mapping the signatures of an -omics experiment to a KG, we can not only verify which causal interactions are observed within a specific context, but also prioritize and identify the mechanism of action of a drug for a given disease with high precision.

Currently, there exist numerous algorithms that leverage causal relations for the interpretation of -omics data. In general, these algorithms operate by assessing the concordance between transcriptomic or proteomic signatures and the predicted causal effects encoded in these relations [18-19]. For instance, the Reverse Causal Reasoning (RCR) [20] and Network Perturbation Amplitude (NPA) algorithms [21-22] assess and score this concordance employing causal graphs consisting of up-stream and down-stream proteins (nodes) representing regulations occurring in biological pathways. Subsequently, the scores obtained from these algorithms can be used for the interpretation of -omics data commonly derived from contrast experiments. Although the interpretations obtained from these algorithms may be relevant for several downstream applications, such as drug target prediction, disease characterization, and side effect prediction, the algorithms themselves cannot be directly used for these applications. Additionally, these algorithms have been specifically designed for bipartite graphs, thus, simplifying biological pathways to a single relation between two proteins.

While traditionally, these algorithms were applied on small causal networks, they have recently begun to be applied on large-scale KGs, given the increasing availability of causal information, including proteins, drugs and phenotypes. For instance, a recent algorithm we published, drug2ways, reasons over all paths between a drug and a disease in a KG to predict the effect of the drug as the cumulative effect of all directed interactions between these two nodes [23]. Reasoning over all paths overcomes the limitation of earlier algorithms that exclusively account for shortest paths on protein-protein interaction networks, oversimplifying the effect exerted by one node on another, as all other paths between the two nodes are ignored [24-25]. Nonetheless, paths in large-scale KGs can grow exponentially, many of which may not be relevant in a true biological context. Thus, incorporating signatures from context-specific experimental datasets along with prior knowledge in a KG can enable us to reason over the entire network-structure and ensure only paths which can be observed in a biologically meaningful context are retained. In doing so, we can address several of the limitations of the above-mentioned methods for drug discovery.

Here, we present RPath, a novel algorithm that prioritizes drugs for a particular disease by reasoning over causal paths in a KG, guided by both drug-perturbed and disease-specific transcriptomic signatures (. We demonstrate how RPath is able to recover a large proportion of clinically investigated drug-disease pairs on multiple transcriptomic datasets and KGs, performing better than other network-based methods. Furthermore, we show two additional applications where we illustrate how our approach can also assist in hypothesis generation and target prioritization.

Schematic representation of the RPath algorithm.

Step 1) All acyclic paths of a given length between a drug and a disease in the KG are calculated. If there exist causal acyclic paths connecting the drug and the disease, a subgraph involving all these paths is inferred. This subgraph represents the proposed mechanism of action by which the drug may be a therapeutic target of the given disease. Step 2) Transcriptomic signatures observed from a drug-perturbed experiment are overlaid onto each corresponding node present in these paths. Then, RPath traverses through each path and evaluates whether the inferred direction of regulation (i.e., activation or inhibition) at every step is concordant with the up- and down- regulations (i.e., red and green nodes, respectively) observed in the transcriptomic signatures. Step 3) In a similar manner, transcriptomic signatures observed within a specific disease context are overlaid onto each corresponding node in the concordant paths from the previous step (if any). Next, RPath evaluates whether the disease transcriptomic signatures contradict the paths that were concordant with the drug signatures. If this is the case, the specific drug-disease pair is prioritized.

Results

This section is divided into three subsections that outline the different applications of RPath presented in this manuscript. First, we demonstrate how RPath can be used to identify potential drug candidates for various diseases using a variety of KGs and datasets, outperforming numerous link prediction methods. Next, we leverage the inherent interpretability of KGs to generate hypotheses for the predictions made by RPath. Finally, we outline how RPath can be reversed-engineered and alternatively used to predict targets for a given disease.

Identification of drug candidates

To demonstrate the ability of our algorithm to accurately identify drug candidate for a given disease, we evaluated its performance to recover clinically investigated drug-disease pairs using two distinct KGs and four transcriptomic datasets (i.e., two each containing numerous drug-perturbed and disease transcriptomic experiments). In this task, RPath consistently prioritized a significantly larger number of clinically investigated drug-disease pairs across all datasets and in both KGs compared with the precision expected by chance (i.e., probability of randomly picking a positive label among drug-disease combinations that are connected through a path) (.

The highest precision values were found for the L1000-GEO datasets with 80% and 66.67% for the OpenBioLink and custom KGs, respectively. In the remaining datasets, the precision was approximately 50%, except for the CREEDS-Open Targets datasets in the custom KG that exclusively yielded a single drug-disease pair which was not in clinical trials. While the precision expected by chance approximately varied between 10% and 42%, RPath consistently achieved higher precision values across nearly all datasets and KGs, ranging between 50% and 80% (e.g., more than five times higher for the L1000-GEO dataset in OpenBioLink running RPath (80%) vs. chance (17.42%)). Notably, the number of prioritized drug-disease pairs were constrained for two reasons: i) RPath requires transcriptomic information for a given drug and disease and, ii) the pair must also be present in the KG (see . Furthermore, apart from the low number of drug-disease pairs that fulfilled these criteria, RPath filters the majority of pairs with paths between them after overlaying the transcriptomic signatures in Steps 1 and 2 (see of the algorithm (. For example, in the case of the CREEDS-GEO datasets and the OpenBioLink KG, the total number of diseases was 10, resulting in only a couple of drug-disease pairs being prioritized. Nonetheless, we were still able to validate our methodology across multiple datasets and KGs, observing that RPath performed significantly better than chance at identifying clinically investigated drug-disease pairs.

Finally, we benchmarked RPath against 11 alternative methods [30-31] that have been used to predict drug-disease links in a KG with the same characteristics as the ones used in this work. The precision of these methods varied between 5% and 43% (. Furthermore, since the majority of these methods prioritize a drug and a disease based on their network proximity (e.g., shortest paths and number of shared nodes), these methods recurrently prioritized the same set of drug-disease pairs. Thus, these methods could not be used to prioritize drugs outside the vicinity of disease-associated proteins since only a minority of drug-disease pairs are connected by a single protein, but the majority of them contain longer paths that are not considered by these methods ( Lastly, we also conducted permutation experiments where we permuted both the binarized gene expression values (i.e., +1, -1, 0) observed in the transcriptomic datasets and the edges of the KGs, while maintaining their underlying structure. The results of our experiments showed how the number of prioritized drug-disease pairs significantly decreases when permuted datasets and KGs are employed and that none of these few prioritized pairs were clinically investigated ().

Interpretation of the mechanisms of action of the proposed drug candidates

In a case study, we sought to explore the results obtained by running RPath on the custom KG. Of the prioritized drug-disease pairs (see , we studied the paths between two of the pairs to demonstrate how our approach can potentially be used to deconvolute the mechanism of action of some drugs ( We selected bicalutamide and ponatinib as these two anti-cancer drugs were the top-ranked prioritized drugs and already approved for prostate cancer and acute myeloid leukemia, respectively. Furthermore, since the mechanisms of action of these drugs have been widely studied, we can compare the mechanistic paths identified by RPath against known interactions and pathways reported in scientific literature.

Devoncoluting the mechanism of action of a drug through RPath.

By investigating all the paths of a given length between a drug and a disease in a KG, we can analyze the different mechanisms that are proposed by RPath. a) Visualization of the custom KG. Proteins are colored in blue, diseases in red and drugs in green. Sankey diagram illustrating a sample of the paths between ponatinib and AML (b) and bicalutamide and prostate cancer (c) for the custom KG. Activatory relations in the Sankey diagrams are colored in red and inhibitory relations in blue.

First, we investigated ponatinib, a multi-targeted tyrosine-kinase inhibitor, which is used to treat acute myeloid leukaemia (AML) (. Among the targets of this drug present in the concordant paths for this pair, we were able to identify fms-like tyrosine kinase 3 (FLT3), which is mutated in approximately 20% of AML patients [32] and several members of the FGFR family proteins. Furthermore, we observed other proteins including KDR, LYN, and SRC, all of which are kinase-associated targets in AML. As a downstream target of these proteins, we found JAK2, a well-studied player in myeloproliferative diseases, with known mutations and hypermethylation events. We further identified the transcription factor, CEBPA, that is critical for normal development of granulocytes and is also implicated in AML [33] and the SPI1 gene, from which circSPI1, a circular RNA derived from the gene, has recently been shown to be highly expressed in AML patients [34]. Other proteins that are inhibited as a result of the signaling cascade triggered by ponatinib include KIT, which is implicated in cell death in AML [35]. RAS family members NRAS and KRAS, both of which are associated with the prognosis of solid tumors and hematological malignancies, including AML [36] were also implicated.

The second studied drug-disease pair is bicalutamide, used for the treatment of prostate cancer. Bicalutamide is an anti-androgen medication that binds to the androgen receptor (AR), as illustrated in . The paths between bicalutamide and prostate cancer point to several downstream targets of this drug, including the epigenetic regulator KMT2D, which is known to sustain prostate carcinogenesis by epigenetic mechanisms [37], and NECAB3, known to enhance the activity of HIF1A, thus promoting glycolysis under normoxic conditions and enhancing tumorigenicity in cancer cells [38]. Furthermore, we were able to identify CTNNB1, which plays a role in the development of numerous prostate cancers [39]. Interestingly, we also observed novel players that have not yet been reported in the literature, such as GNAI1, SYMPK, UBR5, and MEF2C that may provide new insights on the mechanism of action of this drug.

Target prioritization

Prior to the identification of a therapeutic drug candidate for any given disease, a crucial first step is often to identify biologically relevant protein targets. Ideally, the perturbation of a particular protein target in a disease state should result in the reversal of the observed phenotype. In a similar manner to the above-mentioned applications, by reasoning over the KG guided by disease signatures, RPath can be used for target prioritization. Since, as per our knowledge, there are no large datasets that contain information about known targets for a wide variety of indications, we could not conduct a validation strategy similar to the analyses presented in the subsection Identification of drug candidates. Instead, we focused on evaluating the top prioritized protein targets across all diseases using literature evidence (.

Among the top protein target-disease pairs proposed by RPath, two have already been associated with AML, including PRKCA, for which several drugs already exist [40-41] and CXCL8/IL-8 [42-44]. Furthermore, CDC42, which has been proposed as a candidate target for medulloblastoma, plays a role in several cancers. Specifically, CDC42 has been shown to act as a regulator of medulloblastoma-associated genes [45] and compounds for its inhibition have also been proposed [46].

Discussion

In this work, we present a novel methodology that leverages prior knowledge from causal relations across multiple biological modalities in KGs and assesses their concordance with transcriptomic signatures for drug discovery. In the past, several algorithms have been primarily introduced for the interpretation of transcriptomic signatures by reasoning over shortest paths [24-25] or bipartite graphs [20-22]. Though these algorithms could also be indirectly applied for drug discovery, they present some shortcomings: i) they operate on homogeneous causal graphs with a single entity type (e.g., protein nodes), ii) they are solely conducted on single contrast experiments (e.g., drug-treated vs. control), and iii) they do not fully exploit all possible paths in these causal graphs. RPath addresses these shortcomings by reasoning over all possible causal paths in a multimodal KG and leveraging both drug and disease transcriptomic signatures. First, our algorithm reasons over the ensemble of paths between a given drug and a disease in a KG. Second, it evaluates the concordance of these paths against the transcriptomic changes experimentally observed for that drug. Third, it assesses whether the effect of these paths is opposite to the transcriptomic signatures observed within the disease context. In a final step, the algorithm identifies potential drug candidates as those whose paths correlate with drug-perturbed transcriptomic signatures and are anti-correlated to the disease transcriptomic signatures. We have validated our methodology in eight independent analyses, finding that RPath consistently identifies a large proportion of clinically investigated drug-disease pairs over multiple datasets and KGs. Additionally, we conducted several robustness experiments and benchmarked the algorithm against 11 network-based methodologies. Finally, we also showed how our approach can be used to deconvolute the mechanism of action of a drug as well as to prioritize protein targets for a given disease.

We acknowledge a few shortcomings in our work that are worth discussion. Firstly, we were limited by the availability of high-quality annotated transcriptomic datasets for drugs and diseases, as only four of the approximately 30 datasets that we identified met our requirements. Furthermore, the coverage of measured genes varied largely across experiments. For instance, while the average number of genes measured in the Open Targets dataset was approximately 900, that number dropped to 500 in the CREEDS dataset (. In contrast, the total number of proteins in the KGs were in the range of several thousands. As RPath requires that signatures from both these drug and disease datasets be mapped to the KG, most of the proteins in the KGs could not be quantified. Thus, we allowed for up to one error when calculating the concordance in the path between a drug and a disease. Furthermore, two other reasons justified an error within the path. Firstly, introducing an error limits the impact of an arbitrary fold change cut-off, which ultimately determines the up-/down-regulation of each protein. Secondly, some paths might contain causal relations that do not reflect a change at the transcription level of the affected protein (e.g., phosphorylation of a protein kinase) [18-19]. We expect that this challenge we faced of quantifying proteins in our KGs will be overcome by high-quality, consistent datasets such as those generated in large pharmaceutical enterprises and emerging data-driven biotech companies looking to leverage large-scale computational technologies. Another characteristic of our approach is that the identification of a potential drug for a given disease requires knowledge of the protein target and the effect of the drug on it. However, this information is not always available or must be inferred using computational approaches. Finally, the interpretation of the mechanism of action of a proposed drug with the help of scientific literature comes with the caveat that the individual interactions were themselves derived from the literature. Nonetheless, it is still possible to interpret the mechanism of action of a drug irrespective of the aforementioned limitation as the paths of the proposed drug-disease pairs include only those which are concordant with observed data-driven transcriptomic signatures.

While we have demonstrated our novel algorithm across multiple datasets and KGs, we envision multiple other applications. Firstly, by incorporating time series data into the analysis, we can determine how the paths between the drug and the disease are altered over time following the concept outlined by [47]. Secondly, although we have demonstrated our methodology using transcriptomic data, other modalities can be used if the KG contains causal relations for these entities (e.g., metabolomics). Additionally, although we have employed transcriptomic signatures in this work, we acknowledge that RNA levels may not directly reflect the functional activity of proteins. However, given the growth in the availability of proteomic data, we envisage the application of our approach on proteomic experiments from databases such as PRIDE [48], ProteomicsDB [49], and L1000 [26] in the future. Furthermore, although a multimodal KG may lack the context within which each relation occurs, RPath inherently takes this into account by removing the paths which do not match the observed transcriptomic signatures. However, the algorithm could also be applied on a disease-specific KG in order to model the pathophysiological mechanisms characteristic of a given phenotype [50-51].

Methods

Theoretical background

We denote a KG as a set of nodes and edges, where nodes correspond to three distinct biological entities (i.e., chemicals, proteins, and diseases) connected through causal relations, representing activatory or inhibitory effects. Causal relations within the KG connect drug-protein, protein–protein, and protein–disease nodes. A (directed) path in a KG is defined as a sequence of two or more biological entities connected through causal relations. Paths in the KG can be either cyclic or simple. A cyclic path refers to paths in which one or more nodes repeat, whereas a simple path corresponds to a path in which no nodes appear more than once. The length of a path is defined by the number of edges that connect the nodes within the path.

RPath algorithm

The algorithm used in our framework, RPath, reasons over the paths in a KG to identify all possible effects a given drug can have on a disease (. Each of these paths can be divided into three main sequential parts that attempt to represent the mechanism of action of a drug: i) the drug activates/inhibits a protein target (drug-protein edge), ii) the protein target triggers a signaling cascade (a set of protein-protein edges), and iii) the signaling cascade reverts the disease condition (protein-disease edge). Furthermore, since every causal edge contains information on the effect each node exerts on another (i.e., activation or inhibition), we can infer the direction of regulation (i.e., up-/down-regulated) for each node at each step of a path [24-25].

Once the causal acyclic paths between a particular drug and disease in the KG have been calculated (, the next step of RPath is to overlay transcriptomic signatures from a drug-perturbed experiment (. We hypothesize that because a number of paths might represent the biologically relevant mechanism of action of this drug, the observed transcriptomic signatures for proteins in the KG should be concordant with the inferred up- or down-regulations at every step of the path. For example, if in a given path, a drug inhibits a protein target and that target activates a signaling cascade, we expect the inhibition of the protein target as well as the inhibition of the proteins downstream of the target. We would like to note that a gene is considered to be differentially expressed if its expression is significantly altered with respect to a reference sample (i.e., control). Keeping this in mind, a cut-off is applied to each measured gene in the experimental dataset based on the fold change; this measurement is used to define whether the gene is up-/down-regulated or unchanged.

Similarly, the final step of RPath involves overlaying disease-specific transcriptomic signatures to the nodes in the paths of the KG (. We hypothesize that, in contrast to the overlaying of drug-perturbed signatures, transcriptomic signatures in a disease context should be anti-correlated to both the drug-perturbed signatures as well as the inferred up- or down-regulations for every node in the path. This final step is inspired by previous work that exploited the anti-correlation between drug and disease signatures at the pathway level for drug repurposing [15-16]. In summary, RPath aims at prioritizing a specific drug for a given disease if i) there exist causal paths between the drug and disease in the KG, ii) the causal effects on these paths are aligned with the transcriptomic changes observed in the drug-perturbed experiment, and iii) both the drug signatures and the paths are anti-correlated with the transcriptomic dysregulations observed in the disease. outlines the pseudocode of the described logic of the algorithm.

Pseudocode of the RPath algorithm.

Given a KG, drug, disease and a defined path length (i.e., lmax), the core function of the algorithm, is_drug_prioritized, returns whether a drug should be prioritized or not. For this, the function calculates all acyclic paths between a drug-disease pair in the KG. For each path found, drug-perturbed (i.e., drug_tr) and disease-specific (disease_tr) transcriptomic signatures are overlaid onto their corresponding protein nodes. The function then prioritizes the drug if at least one path is concordant with the observed drug-perturbed transcriptomic signatures (evaluated via Function 1, is_concordant) and the same path is anti-correlated with the observed disease-specific transcriptomic signatures (evaluated via Function 2, is_anti_correlated). Paths which match both the drug-perturbed signatures and contradict disease-specific signatures are then returned by RPath as promising drug candidates.

As an additional application, the algorithm can be modified following the same logic for target prioritization (see This variant of the algorithm begins from a disease of interest and calculates all paths from the disease to all proteins for a given path length (e.g., a path length of 6). Next, it calculates the concordance between the paths for each potential protein target and the transcriptomic signatures of the given disease to assess whether there are proteins that could be key up-stream regulators of the observed phenotype. We would like to note that this application exponentially increases the running time of the algorithm as it requires querying paths from a disease to several thousands of proteins in the KG, as opposed to only a handful of chemicals.

Datasets and validation

In this subsection, we present drug-perturbed and disease-specific transcriptomic datasets as well as the KGs used to demonstrate our methodology. We then introduce the strategy we follow to validate our methodology.

Drug-perturbed and disease transcriptomic datasets

We identified four databases that were suitable for our approach (; drug-perturbed transcriptomic data were obtained from CREEDS [29] and L1000 [26] while disease transcriptomic data were collected from Open Targets [28] and GEO [27]. All experimental datasets from these resources (downloaded on 15.02.2021) contained gene expression changes measured in humans. Drugs and diseases from datasets obtained from these databases were then mapped to PubChem compound identifiers and the Mondo Disease Ontology (MONDO), respectively, for consistency with the entities of the knowledge graphs presented in the next subsection. Similarly, gene identifiers in all datasets were harmonized to ENTREZ. Of the four databases, datasets from L1000 contained a binarized value for the direction of dysregulation for every gene (i.e., up-regulation and down-regulation), while for the remaining databases, fold changes were binarized for significantly dysregulated genes using |log2 fold change| = 1 as a cutoff (. As fold change thresholds tend to be arbitrary selected [52], we opted to select a threshold of 1 as opposed to a more stringent one (e.g., |log2 fold change| > 2) to ensure a larger number of dysregulated genes would be retained. Finally, we conducted a systematic search for databases that contained either a large number of drug-perturbed or disease-specific transcriptomic datasets. While this search initially resulted in 27 candidate databases (see , the majority of them were not suitable for our study as they either contained too few transcriptomic datasets or the drugs/diseases in these datasets were not in the KGs used to demonstrate our methodology.

Knowledge graphs

We demonstrate our methodology using two established publicly available KGs that contain causal relations across drugs, proteins, and diseases: OpenBioLink KG [53] and a custom KG [23]. Both KGs are originally generated from a compedia of independent databases; thus, containing unique causal interactions depending on the source databases they include. As outlined in the algorithm, the KGs are required to encompass three types of causal edges: drug-protein (i.e., drug activates/inhibits protein), protein–protein (i.e., protein activates/inhibits protein), and protein–disease (i.e., protein activates/inhibits disease). Furthermore, the original node identifiers for drugs and diseases in both KGs were respectively mapped to PubChem compound identifiers and MONDO to be consistent with the transcriptomic datasets. Next, we removed drugs and diseases that were not present in any of the four transcriptomic datasets presented in the previous subsection as the paths between these drug-disease pairs cannot be validated. shows the final statistics of both KGs after the previously outlined filtering steps. summarize the overlap between the genes measured in each of the four transcriptomic datasets and their corresponding protein nodes in the KGs.

Distribution of node and edge types in the custom and OpenBioLink KGs.

The properties of each of the two networks are detailed in .

Validation

In line with other network-based approaches designed for drug discovery [54-55], we used drug-disease pairs that have been clinically investigated as positive labels, extracting this information from ClinicalTrials.gov (accessed on 28.09.2020). Clinical trials are commonly used as a proxy for highly validated, medically relevant biological interactions, independently of whether the clinical trial was successful or not, as multiple in vitro and in vivo studies must first validate the interaction in order for the drug to proceed to a clinical trial. Notably, this assumption may result in a worse performance being reported than the actual performance of the algorithm as some of the drug-disease pairs that are considered as negative labels may in fact be positive ones.

Since drugs and diseases in ClinicalTrials.gov are formalized using MeSH identifiers, we harmonized these identifiers to the ontologies used in the KG (i.e., PubChem compounds for drugs and MONDO for diseases) (mappings are available at the GitHub repository). After the harmonization, any drug or disease in the KGs that was not present in any clinical trial or did not have any path to any disease in the KG was subsequently removed, as its corresponding node could not be used in the presented validation. Details on the harmonization procedure are provided in

To test the robustness of RPath in identifying these clinically investigated drug-disease pairs, we conducted eight independent analyses, one for each of the combinations of the two drug datasets, the two disease datasets, and the two KGs (e.g., CREEDS-GEO-OpenBioLink, L1000-GEO-OpenBioLink, etc.). For each of these eight analyses, we ran RPath over a given KG to prioritize drug-disease pairs among all possible drug-disease combinations. We would like to note that these pairs prioritized by the algorithm are those whose paths are both correlated with drug-perturbed transcriptomic signatures and anti-correlated with disease transcriptomic signatures. Furthermore, we make two assumptions in the design of the algorithm. Firstly, paths with cycles or a length greater than 7 edges between a given drug and disease are not considered, assuming that the effects exerted by paths beyond this length are less biologically relevant [23]. Secondly, we allow for at most one error between the transcriptomic data and a given path (see pseudocode of the algorithm in . We refer to an error in the path as the disagreement between the type of causal interaction (i.e., activation or inhibition) and the direction of dysregulation of genes in the transcriptomic datasets (i.e., up or down -regulation), or the absence of the drug-perturbed and/or disease-specific transcriptomic signature. We restricted the number of allowed errors to at most one as, without any errors, running the algorithm over the two KGs with most dataset combinations will not yield any prioritized pairs, and permitting more than one error will result in an exponential increase in the number of prioritized pairs. For example, in the latter case, for a path of length 5 (i.e., a sequence of 3 proteins), an allowance of two errors (e.g., missing expression values for 2 of the 3 proteins) would still result in the prioritization of the drug-disease pair if the remaining protein both correlated with the drug and anti-correlated with the disease, obfuscating results.

From this set of prioritized drug-disease pairs, we expect to retrieve a larger proportion of clinically investigated drug-disease pairs (i.e., positive labels) than expected by chance (i.e., proportion of positive labels in the dataset that also have a path between the drug and disease). Here, it is important to note that, as in any drug discovery task, there is a class label imbalance where the vast majority of the drug-disease pairs are negative labels while the proportion of positive labels ranges from anywhere between 9% and 41% for each of the eight analyses (. Furthermore, this type of validation falls into the so-called early retrieval problem. In other words, from the thousands of drug-disease pairs that are tested, we are exclusively prioritizing the top-ranked pairs that have been equally prioritized by the algorithm. This small subset represents the interesting drug-disease pairs that would be further investigated in the drug discovery process. In such cases, it is inadequate to apply metrics such as receiver operating characteristic (ROC) curves as they operate on a full ranked list. Therefore, it does not necessarily evaluate the ability of a model to prioritize the most promising drug-disease pairs candidates [56]. Additionally, considering that not all drug-disease pairs have been clinically studied, a number of the negative labels might be falsely classified as positive. To address these issues, we evaluated the performance of RPath based on the ratio of true positives that appear in the prioritized drug-disease pairs (i.e., precision was used as the performance metric). As a baseline, we assessed whether the prioritized drug-disease pairs found through the algorithm contain a larger proportion of positive labels (i.e., drug-disease pairs investigated in clinical trials) than expected on average by chance.

As a benchmark, we compared RPath against 11 equivalent approaches that can be used to prioritize drug-disease pairs based solely on network structure, as outlined by [26] and [27] (. The choice of these approaches is motivated by the fact that, as per our knowledge, there are no network-based methods that operate on multimodal KGs using transcriptomic signatures for the prioritization of drug-disease pairs. Additionally, we conducted a validation experiment where we simultaneously randomly permuted the directionality of the genes measured in the transcriptomic datasets and the KGs using the XSwap algorithm [57] while both preserving network structure and the original gene expression distributions. Using these, we then rerun the eight analyses to compare the significance of our results [57].

Implementation details

The RPath algorithm and the benchmarked methods are implemented in Python leveraging NetworkX (v2.5) (https://networkx.github.io). Network visualizations were done using WebGL, D3.js, Three.js, Matplotlib and igraph. Source code, documentation, and data are available at https://github.com/enveda/RPath. The validation presented in the paper can be reproduced by running the Jupyter notebooks available at https://github.com/enveda/RPath/tree/master/notebooks.

Processing of transcriptomic datasets.

(DOCX)

Click here for additional data file.

Benchmarked methods.

(DOCX)

Click here for additional data file.

Prioritized pairs.

(DOCX)

Click here for additional data file.

Clinical trial information mapped to the OpenBioLink and custom KGs.

For each possible pairing of a drug-disease database (i.e., column 2) with entities that could be mapped to either the OpenBioLink or custom KG, we report the proportion of drug-disease pairs contained in ClinicalTrials.gov (i.e., column 3). We consider these clinically investigated drug-disease pairs as positive labels for the validation of our approach. The number of unique drugs (PubChem compound identifiers) and diseases (MONDO identifiers) are reported in columns 4 and 5, respectively, while the total number of possible combinations of these unique drugs and diseases are presented in column 6.

(XLSX)

Click here for additional data file.

Percentage of drug-disease pairs at different steps for every dataset combination and KG.

(XLSX)

Click here for additional data file.

Evaluation of the 11 benchmark methods using precision as a metric.

None of the benchmarked methods achieve a precision greater than 50%. Furthermore, we would like to note that most of the drug-disease pairs prioritized by each of these methods are the same since they are based on network proximity. Thus, if a drug and a disease share a large number of nodes, they will consistently be prioritized by most of these methods.

(XLSX)

Click here for additional data file.

Statistics on the genes measured in the four transcriptomic datasets used.

(XLSX)

Click here for additional data file.

Investigated datasets.

(XLSX)

Click here for additional data file.

Overlap of the transcriptomic dataset and the KGs.

The total number of drugs and diseases in the datasets which can be mapped to the KG as well as the proportions of them that are present in the KG are given in columns 3 and 4, respectively. Similarly, column 5 displays the total number of mapped proteins as well as the proportions of proteins that are present in the KG. Details about each individual drug/disease are available at https://github.com/enveda/RPath/blob/master/data/drug_disease_overview.tsv.

(XLSX)

Click here for additional data file.

Properties of the OpenBioLink and custom KGs.

(XLSX)

Click here for additional data file.

Drug repurposing approaches exploiting anticorrelation of transcriptomic signatures.

(XLSX)

Click here for additional data file.

Permutation experiments across the four dataset combinations using permuted KGs and gene expression datasets.

(XLSX)

Click here for additional data file.

Pseudocode of the RPath algorithm designed for target prioritization.

(TIF)

Click here for additional data file.

29 Nov 2021

Dear Mr Domingo-Fernández,

Thank you very much for submitting your manuscript "Causal reasoning over knowledge graphs leveraging drug-perturbed and disease-specific transcriptomic signatures for drug discovery" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

As you will see from the reviews, both reviewers, but especially Reviewer 2, have considerable concerns regarding the quality of the validation of your methods. They make several suggestions how this could be improved, and we would recommend carefully considering these options. Also, several points of the procedure appear to be somewhat arbitrary (cutoffs, etc...). These should be extensively explained and justified.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Carl Herrmann, Ph.D.

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: This paper describes a novel algorithm that overlays knowledge graphs with drug-perturbed signatures and disease-specific signatures as constructed from 4 data sources. The description of the four data sources can be extended: it is important to know how many drugs, diseases and genes are covered and how much overlap between them (and the KG).

The validation section leaves space for a number of questions and unclarities that need to be addressed. In the description of the algorithm two assumptions/restrictions are described. Paths with cycles or length > 7 are removed. How many are there found? Furthermore, one error in the path is permitted. What does this mean? What is an error? An inconsistent edge? Again, how often does this occur?

For the validation drug-disease pairs are obtained from ClinicalTrials.gov. They are considered as positive labels, even if the ClinicalTrial was aborted? What about negative labels? Did you consider to use a set as a proxy to negative cases? How many drug-disease pairs from CT remained after the harmonisation step? How many of those pairs were found in the prioritisation step (correlation with drug-perturbed signatures and disease signatures)?

Do I read well from table 1 (looking at the TP+FP) that the number of drug-disease pairs validated was resp. 5, 11, 2 and 2? If that is the case this validation is quite limited. If this is not the case more explanation is needed to understand the figures. Why not running all drug-disease pairs from CT against the different sets?

For the interpretation task tow drugs are selected: ponatinib and bicalutamide. How and why have these two drugs selected?

For the target prioritization task it would be better to show the results for all protein target-disease pairs for both Custom and OpenBiolink KG. It is not clear what determines which KG is selected.

Figure 4 is unreadable for me. Too small to see what texts are hidden in the Sankey figure. In the text I think erroneously 4a was referenced for the potatinib case (should be 4b) and 4b for the bicalutamide case (should be 4c).

Baseline: the approach is nice and interesting, the validation needs more attention and requires also a better understanding of the selections/choices made.

Reviewer #2: Review PCB 2021

Domingo-Fernandez et al present a method to prioritise drugs by combining knowledge graphs (KGs), and transcriptomics profiles. The motivation and background are very reasonable, and authors outline and the paper does a good job at pointing directly at the shortcomings we have to deal with at the state of our current knowledge. However, we found that the work dies not find good solutions to actually address them and both the implementation as well as the paper in itself are quite shallow, as we outline below. It might well be that the exiting knowledge and data we have is not ready to be used in this way.

Major comments:

1- the method is based on a combination of 3 ‘step’s in analysis and feels like it i mixing three very different things that might make sense by themselves and ‘see what happens’.

First, finding pathways in KGs seems sensible for a drug, but then overlying transcirptomics of that drug on those pathways is much less clear - the authors argue:

“the observed transcriptomic signatures for proteins in the KG should be concordant with the inferred up- or down-regulations at every step of the path”.

What is the justification for this? e.g. they elaborate: “we expect the inhibition of the protein target” …. Can authors prove this on the datasets they use?

Indeed, the effect of a drug does not need to affect the genes on the same pathway, and that is in fact the fundament of the methods based on ‘causal reasoning’ .

And the comes yet another overlay of transcriptomics in step 3, now from the disease, but only ’on the concordant paths from the previous step (if any)’

2. Further, a fundamental limitation of mapping gene expression to KGs is that the level where authors have “causal” information (proteins) is largely detached from the levels at which they have measurements (RNA). Can authors run their methods using proteomic data to be closer to the actual KGs?

Of note, even if they had adequate protein-level measurements, their method would still only be applicable to the drugs whose effects remain on the protein level exclusively, which not even their selected examples do.

3. Authors perform a number of seemingly arbitrarily methodological choices that need justification (there is some vague general statements in the discussion):

3.1. - Keeping this in mind, a cut-off is applied to each measured gene in the experimental dataset based on the fold change

Which threshold was used and how chosen?

3.2. - Paths with cycles or a length greater than 7 edges between a given drug and disease are not considered,

Why exactly 7?

3.3… - Secondly, we allow for at most one error between the transcriptomic data and a given path

Why 1?

3.4.why 'non-concordant paths are removed.’ In step 2?

4. Benchmarking. Authors benchmarked RPath against 11 equivalent approaches that can be used to prioritize drug-disease pairs based solely on network structure,

Why not methods that use gene expression? It is not fair to compare against methods that use different input data.

5. Validation and overstating results.

5.1. Calling finding ‘ a larger proportion of positive pairs than expected on average by chance.’ a validation is over-selling the results.

5.2. Evaluating based on literature evidence”, seems that they took the output of their prioritisation run and looked up the proteins in PubMed? If so, they are validating their KG by the very data that the KG very likely was built from, and so finding papers in favour of a connection of their prioritised drugs with the disease is not surprising at all. They don’t even attempt to address any confounding of this kind or research bias in general. The discussion also does not mention it.

6. Lack of enough description, in part relating to the specific points above:

6.1. Section 2.3.1 lacks fundamental description of which diseases they are able to look at in the model, and the remainder of the paper also does not mention the pool of diseases and drugs studied.

6.2. Section 2.3.3 argues at length why they can’t apply standard metrics (class label imbalance etc), but they don’t describe their chosen method. BTW, AUC-PR is not an ROC curve, as their explanation suggests.

6.3. Same can be said for the results, presenting 4/5 as “80%” without a measure of confidence.

6.4. Pharmacologically, not explaining the “hits” is underwhelming. Which drug, which disease, does the mechanism of action make sense? They even mention mechanism of action in their abstract…

6.5. The chosen examples look cherry-picked and are not representative of the general notion of “we predict drug mechanisms”, because both are cancer drugs. Both also exemplify the detachment of inference and measurement, as both have extensive transcriptional effects that are not accounted for in the KG (androgen receptor…).

Based on the points above, authors should tone down their claims.

Minor points

1. - The link to the Jupyter notebooks https://github.com/enveda/RPath/tree/master/src/notebooks is broken

2 - “we harmonized these identifiers to the ontologies used in the KG […] using proprietary harmonization scripts”

Disappointing that this part of the work can not be checked as proprietary and lacks reproducibly and transparency

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No: links to code is broken. harmonization is based on 'proprietary software'

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Erik M. van Mulligen

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

2 Dec 2021

Submitted filename: RPath - rebuttal.docx

Click here for additional data file.

17 Jan 2022

Dear Mr Domingo-Fernández,

Thank you very much for submitting your manuscript "Causal reasoning over knowledge graphs leveraging drug-perturbed and disease-specific transcriptomic signatures for drug discovery" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

While Reviewer 1 is satisfied with your response, Reviewer 2 still has strong reservations, and feels that you failed to adequately address some of their previous comments. We would encourage you to carefully read the comments and address them in a revised version. In particular, you should consider the following points and make appropriate changes:

While I understand that the inclusion of proteomic data goes beyond the scope of your methods, you should at least acknowledge this resource and the fact that corresponding proteomic data is increasingly becoming available, for example through the database mentioned.

Given the number of parameters used in your method, it is indeed questionable how they impact the results. While a full exploration of the parameter space is not realistic, I would ask you to provide evidence how changing the number of allowed errors impacts the results, e.g. as a supplementary table.

Regarding the comparison with other methods: No alternative method will exactly fit the one developed by the authors. The authors have shown in their comparison that including transcriptomic data in the context of KG improves compared to using the network topology alone. They should however also address the complementary question, namely how using KG together with expression signatures improves over transcriptomic-only approaches for the prediction of drug-disease pairs. In particular, numerous studies have tackled the problem of drug repurposing using transcriptomic signatures.  We would urge the authors to at least discuss in depth those methods, and if possible, provide a comparison with some of these transcriptomic-only methods.

Regarding the precision of the method and the comparison with the random null hypothesis: the small number of cases makes the significance of the precision indeed questionable. An alternative null hypothesis should take into account the number of drug diseases-pairs that are already connected in the KG, and not all drug x disease pairs, as the algorithm already selects for connected pairs.

I would urge the authors to provide the full list of 30 drug-disease pairs, and discuss these identified hits, to support the claim of generality beyond the two examples already discussed.

Please consider uploading the supplementary files to zenodo and providing the corresponding doi as an alternative to github, which should be restricted to providing the code.

As previously mentioned, I would strongly recommend considering the other issues raised by reviewer 2 beyond those listed above.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Carl Herrmann, Ph.D.

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Thanks for the response on my previous remarks. Nice manuscript that deserves to be read and cited a lot.

Reviewer #2: We acknowledge the efforts of the authors to address our points, but we feel that they could and should have put more effort in doing so. Some comments are simply dismissed, as showed below. The fact that the revised version was sent back in only a few days suggests that indeed not too much effort was put in revising the manuscript.

As bottom line, the paper tries to introduce KG-based learning as drug prioritisation tool, but fails to convince us that this can be actually useful, given the current state of our knowledge and the algorithm at hand. The paper for us is too detached from clinical reality and does not try to get into the details of how a drug works and why this could be mirrored by their algorithm. If the “validations” are not significantly extended, we do not think that one should believe that these results can explain drug mechanisms of action and prioritise drugs for new diseases.

Below some comments on the answers to our points:

# Some questions have been misunderstood it seems and thus the answer does not relate to the question. For example to the question:

Indeed, the effect of a drug does not need to affect the genes on the same pathway, and that is in fact the fundament of the methods based on ‘causal reasoning’. And then comes yet another overlay of transcriptomics in step 3, now from the disease, but only ’on the concordant paths from the previous step (if any)’

Authors answered:

* Although the effect of a drug does not always have to affect the genes on the targeted path, it is reasonable to think that for the majority of the cases, an effective drug for a given disease will likely alter the signatures of the targeted path. The idea of overlaying transcriptomic signatures from diseases and drugs is widely established and it is the basis of well-known drug discovery libraries such as CMap and L1000.

But The paradigm is to match indeed signature, but that was not the question, rather why matching the expression to the path where the drug targets.

# An example of how authors are dismissive about our comments is the question about proteomics vs transcriptomics, which they answer with factually wrong answers:

1 Proteomics data does not yet have the same coverage as RNA-seq or other gene expression techniques that we have used in our study (thousands of different transcripts are measured vs. a few proteins).

Mass-spectrometry proteomics provides coverage of thousands of proteins. It might not be genome wide, but there is enough coverage to at least try.

1. As of yet, there are no major repositories of proteomic datasets as there are for gene expression (e.g., GEO or ArrayExpress).

This is incorrect. See e.g. PRIDE, https://www.ebi.ac.uk/pride/archive/

# Methodological choices of e.g. parameters should have bene further justified or at least explored the implications of changing them. For example,

3.3... - Secondly, we allow for at most one error between the transcriptomic data and a given path. Why 1?

Using zero as a threshold, there were no drug-disease pairs prioritized because the chance of finding a perfect concordant path is extremely low (see reasons in the second paragraph in the discussion section). Thus, we decided to allow for at most one error, precluding the possibility also of too high a leniency.

This seems an arbitrary decision, why not exploring e.g. 2, 3,…?

# Comparisons to other methods should be further elaborated:

As per our knowledge, there are no network-based methods that operate on multimodal KGs using both transcriptomics modalities and that also have software available that can be used to prioritize drug-disease pairs.

Can authors not use any network base method for gene expression that was perhaps not developed for KGs but for molecular networks and extend to KGs?

# Some answers are not clear, e.g.

> Nonetheless, because the majority of the interactions in the KG are not relevant to any particular disease condition (except for the last edge in the path connecting a protein to a disease), the interactions encoded in a KG do not directly reflect literature knowledge, nor favour connections of prioritized drugs with the disease.

This is a confusing statement “KG do not directly reflect literature knowledge,” we thought they do reflect literature knowledge, that is the nature of a knowledge graph?

#validations are still unconvincing in their scope:

> As mentioned in the previous point, we used precision as a metric. In this case, 5 drug-disease pairs were prioritized and among them 4 have been clinically investigated, resulting in the reported precision of 80%. Furthermore, we demonstrate in several datasets and two different KGs that, nearly across the board, the precision the algorithm achieves is better than the expected precision by chance (Table 1).

In a small set of 5 “validators”, the chance for spuriously high precision is quite significant. Considering that the algorithm uses gene-disease relationships (final connection) and allows only one mismatch, it effectively filters for a set of paths that have higher likelihood of being disease-related already. That the precision here is higher than average should be self-explanatory. Combined with the multitude of possibilities given by omics measurements, the chance of spurious association seems too high (see Ioannidis 2005, https://doi.org/10.1371/journal.pmed.0020124).

> There are 30 prioritized pairs whose investigation would require a considerable manual effort. Thus, we focused extensively on the proposed mechanism of actions of only two of these drugs in the Interpretation section 2.2.

Since the authors claim generality and also acknowledge that there is no baseline to compare their method against (see above), they should supply proof of this generality. Otherwise, at least tone down the claims.

For the reader, it is not apparent whether the chosen two results are representative or cherry-picked. This is of particular relevance because the abstract states that "The paths which match this signature profile are then proposed to represent the mechanism of action of the drug.”

For someone knowledgeable in pharmacology, it is not a lot of effort to sanity-check 30 predictions. Also, why can the results not be given as a supplementary table, such that the knowledgeable reader can form their own opinion of how meaningful are these predictions?

For instance, given the now-available pool of drugs and diseases, we would be interested on whether glucocorticoids or 4-ASA show up in the treatment of inflammatory bowel disease or ulcerative colitis. Antidiabetics in diabetes? It is also intriguing that psychiatric diseases are in the list, while psychiatric medications are quite underrepresented in the pool. Instead, there are several anti-infectives, which do not target the human organism at all; and cancers are severely overrepresented as well.

Minor:

# Provide all information in a long-term manner:

> We have now added a file which displays the drugs and disease with respect to their source data annotation on GitHub (https://github.com/enveda/RPath/blob/master/data/drug_disease_overview.tsv). This file is now referenced in the Supplementary Information.

This is not adequate to point to a repository; a manuscript should contain all information that is required for the reader to understand capabilities and implications of the algorithm/method. Additionally, GitHub is not a permanent storage facility, authors should deposit everything relevant for the paper in e.g. zenodo.

* For the sake of reproducibility and transparency of the paper, we released the mappings corresponding to the diseases analyzed in the paper. The mappings are available at https://github.com/enveda/RPath/blob/master/data/manual_drug_disease_mapping.tsv and the notebooks have been accordingly updated. The corresponding sentence has been modified accordingly.

Are the ‘proprietary harmonisation scripts’ also provided or only the results mapping? Otherwise, can the tool be used by others in other contexts?

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No: is unclear of code is available or only resulting mappings. also as noted some data is in a github repository so could get lost or hidden again.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Erik M. van Mulligen

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

24 Jan 2022

Submitted filename: RPath - response to reviewers.docx

Click here for additional data file.

9 Feb 2022

Dear Mr Domingo-Fernández,

We are pleased to inform you that your manuscript 'Causal reasoning over knowledge graphs leveraging drug-perturbed and disease-specific transcriptomic signatures for drug discovery' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Carl Herrmann, Ph.D.

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************************************************

22 Feb 2022

PCOMPBIOL-D-21-01914R2

Causal reasoning over knowledge graphs leveraging drug-perturbed and disease-specific transcriptomic signatures for drug discovery

Dear Dr Domingo-Fernández,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsanett Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol

Table 1

Evaluation of RPath in multiple datasets across the two KGs using precision.

Each row corresponds to the results of running RPath on a specific drug-disease dataset combination. The second and fourth columns show the performance that is expected to be achieved by chance.

-	OpenBioLink KG		Custom KG
Dataset combination	Precision (TP/TP+FP)	Expected precision by chance	Precision (TP/TP+FP)	Expected precision by chance
L1000 [26]–GEO [27]	80% (4/5)	17.42%%	66.67% (2/3)	13.74%
L1000 –Open Targets [28]	54.55% (6/11)	15.01%	50% (2/4)	9.62%
CREEDS [29]–Open Targets	50% (1/2)	32.66%	0% (0/1)	24.40%
CREEDS–GEO	50% (1/2)	41.15%	50% (1/2)	34.08%

Table 2

Top 5 prioritized protein target-disease pairs.

These results were obtained by running RPath over both KGs with the GEO and Open Targets datasets using the same path length as the drug discovery task (see . Pairs were prioritized based on the number of concordant paths. The vast majority of pairs were prioritized using the disease transcriptomic signatures from the GEO dataset given its larger coverage of measured genes compared to Open Targets (.

Protein target	Disease	Concordant paths	Nodes in the concordant paths	KG	Transcriptomic dataset
NOG	AML	18,456	1,008	Custom KG	GEO
PRKCA	AML	12,861	669	Custom KG	GEO
CXCL8 / IL-8	AML	7,234	465	Custom KG	GEO
NOG	Plasma cell myeloma	5743	616	Custom KG	GEO
CDC42	Medulloblastoma	5,651	91	OpenBioLink	GEO