Rare diseases are low-prevalent disorders caused by pathogenic mutations or harmful environmental factors that can have chronic, debilitating, or life-threatening effects [1]. Currently, there are over 7000 rare diseases that affect approximately 10% of the European and American populations, yet fewer than 6% have an approved treatment option [2]. This highlights the pressing need for developing therapies targeting rare diseases. However, the development of a new drug can be a time-consuming and costly process, taking up to 15 years and costing as much as US$2.5 billion [2]. Consequently, the development of novel drugs for rare diseases, which affect only a small number of individuals, is not pursued frequently, as it is less likely to provide a return on investment for pharmaceutical companies [1]. A cost-efficient and faster way to provide drugs for rare diseases is via computational drug repurposing.
Drug repurposing is the process of identifying for an already approved or investigational drug a new use outside the scope of the original medical indication. For example, a drug could be repurposed for a different disease, based on the knowledge that drugs target particular pathways and disease mechanisms that may be shared by multiple diseases [3]. Computational drug repurposing aims to predict novel drug-disease associations, which can be achieved by predicting drug-target interactions (DTIs) [4]. The computational prediction of new DTIs can provide insights into potential pathological and drug mechanisms, as well as drug repurposing and design, helping researchers to generate testable hypotheses in the lab [5]. Network-based data integration and machine learning-based methods for DTIs prediction can mitigate costly and time consuming experimental verifications and are the current state of the art approaches in computational drug repurposing [6,7,8,9,10].
The landscape of biomedical information resources is heterogeneous and broad, yet most current methods for predicting DTIs are limited to homogeneous networks or bipartite models, failing to account for the intricate relationships among diverse data sources [6]. To fully exploit the potential of computational drug repurposing, we propose an automated workflow for predicting DTIs that does take complex relationships among diverse data sources into account. For our DTI prediction, we use the biophysical drug repurposing approach, which is based on the hypothesis that structurally similar drug molecules share similar targets. We extract biological data from multiple online databases using the BioKnowledge reviewer library (https://github.com/SuLab/bioknowledge-reviewer), a tool developed by Queralt-Rosinach et al. [11]. This library integrates heterogenous knowledge and data into a knowledge graph, which is a machine-readable semantic representation of relational information, where concepts are encoded as nodes and relationships between them as edges. Now the prediction of DTIs can be framed as a link prediction problem, where the goal is to identify missing edges in the knowledge graph between drugs and targets that represent potential DTIs. To address this challenge, our proposed automated drug repurposing workflow leverages both network-based analysis and machine learning methods. Network-based methods help to identify potential interactions based on network topology and structural features, while machine learning methods can use more complex data features to make predictions. By combining these approaches, our automated workflow is able to discover new DTIs that can be further used by biologists to generate drug repurposing hypotheses that can be tested in the lab.
The drug repurposing workflow generalises to any rare disease, but as a proof-of-concept, we focus on Huntington's disease (HD). HD is a rare neurodegenerative disorder of the central nervous system characterized by dementia, involuntary movements due to the movement disorder chorea and behavioural and psychiatric disturbances [12]. There are some symptomatic treatments available but because their effects are limited, there is a constant need for better, modifying drugs to treat symptoms of the disease [12]. We present the automated data mining workflow as a web application, using Flask, which makes the workflow accessible for researchers with no technical expertise. The Python code for running this app is accessible as a Docker container, available at https://github.com/carmenreep/ DrugRepurposing, which runs on a laptop with minimal specifications (at least CPU 2.80 GHz). Besides finding missing edges (potential DTIs) in the knowledge graph, the machine learning model also predicts the interaction types of these missing edges. Good embeddings for this link prediction task were achieved through enriched information on the drugs and target sites.
To obtain both human and animal biological data and metadata, the workflow uses the Monarch Biolink API version 1.1.14 (https://api.monarchinitiative.org/api). The Monarch Initiative, a collaborative, open science project, seeks to semantically integrate genotype-phenotype information from numerous sources and species [13]. To integrate drugs into the knowledge graph, the workflow utilizes the Drug-Gene Interaction Database (DGIdb), a web resource that aggregates information on drug-gene interactions and druggable genes from various sources, including publications, databases, and web-based resources [14]. We obtain drug-gene information from DGIdb using its API (version v2) available at https://dgidb.org/api.
The workflow leverages the BioKnowledge reviewer library to extract and integrate data from online sources into a knowledge graph. Starting with a list of seed nodes, a Monarch network is created by including the first layer of neighbours and relations from Monarch for each seed node, along with their ortholog-phenotype nodes. A seed node can take the form of a disease phenotype MIM number, such as '143100' representing Huntington's disease. The edges are formatted as triples, where each triple includes additional information such as the reference Uniform Resource Identifier (URI), the date when the information was obtained, and more information about the semantics of the relation. Nodes in the graph are identified using different biomedical ontologies in the OBO Foundry [https://doi.org/10.1093/database/baab069] maintained or used by Monarch and contain other attributes such as semantic group, URI, label, name, synonyms, and description. The URI serves as a link to a web page that provides a more detailed description of the ontology term representing the node [11].
To obtain drug-target information, we use DGIdb and take all genes (targets) in the Monarch graph as seeds. First, we need to map the Monarch genes to Entrez Gene identifiers (Entrez ID), which are used as a standard gene identifier system [15]. We accomplish this using the BioThings MyGene.info API, accessed with the Python wrapper biothings_client version v0.2.6 [16]. For each gene, we obtain a list of drugs (ID, name) that interact with the gene, along with the type of interaction and interaction source. The drug identifiers are from either the ChEMBL Database [17] or the Wikidata knowledge base [18]. There are various interaction types, such as 'activator', 'blocker', but to improve our predictions, we used the interaction direction (inhibits or activates) instead of the interaction types themselves [14]. Because some relations lack direction, we introduced a third category called 'regulates'. Table 1 shows each interaction direction category along with the interaction types belonging to that category. We mapped these three interaction groups to URIs using the OBO Relations Ontology (RO) https://www.ebi.ac.uk/ols/ontologies/ro, version 2022-05-23.
To enable the biophysical drug repurposing approach, it is necessary to identify structurally similar drugs. In our workflow, we use the Tanimoto coefficient [19] to measure the similarity between drugs. To achieve this, we first retrieve the SMILES chemical structure notation for each drug in our graph using the BioThings MyGene.info API accessed with the biothings_client version v0.2.6 [16]. Subsequently, we convert the SMILES structures into RDKit molecule objects using the RDKit Python package version 2022.3.2 with the Chem module [20]. The RDKit molecule objects are then transformed into Morgan fingerprints using the GetMorganFin-gerprintAsBitVect() function of the AllChem RDKit module. Using the BulkTanimotoSimilarity() function of the DataStruct module from RDKit, we can calculate Tanimoto coefficients between every possible pair of drugs in the graph. This results in a large number of weighted edges, which can lead to a complex network. To mitigate this, we adopt a method by Thafar et al. [21], where all similarity scores are ranked in descending order and only the top-10 most similar drugs are retained, similar to the k-nearest neighbours algorithm. Finally, we label all similarity edges with the 'CHEMINF:000481' ID, 'http://semanticscience.org/resource/CHEMINF_000481' URI, and the human readable string description 'similar to'.
The final graph is transformed into a Resource Description Framework (RDF) graph [22], where each entity, relationship, and entity class (gene, drug, etc.) is represented as an ontology term by its URI. As Monarch did not provide identifiers for the entity classes, we manually mapped the entity class labels to terms in ontologies and used their URIs. To achieve this, we utilized the Semanticscience Integrated Ontology (SIO) [23] (http://semanticscience.org/ontology/sio.owl) version 1.53 and obtained URIs using the URI resolution service identifiers.org (https://registry.identifiers.org/registry/sio, accessed June 2022). Table 2 presents the specific URIs we used for the entity classes. To perform this transformation, we extended the BioKnowledge reviewer by using the RDFLib Python package version 6.1.1 [24] and ensured that the graph was stored in Turtle format [25].
To prepare the knowledge graph for embedding, we first remove all known drug-gene interactions from the graph. This is important to prevent bias in the prediction task, as keeping these edges would make the embedding vectors of the drugs very similar to the embedding vectors of the genes they interact with. We therefore split the graph into two separate graphs, one for drugs and one for genes. Each graph is then embedded separately using a graph embedding algorithm. After the embedding, the drug vectors are fused with the gene vectors to obtain drug-gene edges, which are used for training the machine learning model, as explained in more detail in the XGBoost prediction model section below. Celebi et al. [26] compared different knowledge graph embedding methods for drug-drug interaction prediction, and found that RDF2Vec with Skip-Gram generally outperforms other methods. Therefore, this workflow employs RDF2Vec for graph embedding. RDF2Vec adapts the language modelling approach of Word2Vec to RDF graph embeddings [26]. First, random walks are performed over the graph to generate sequences of entities and relations. Then, the Skip-Gram model is used to learn one embedding for each entity/relationship in the graph. After training, semantically and syntactically similar entities/relationships have similar embeddings [26]. For the prediction task in this study, only the drug and gene vectors are of interest, and therefore, only these vectors were selected for further computation.
For RDF2Vec, the Python function RDF2VecTransformer from the rdf2vec module of the package pyrdf2vec version 0.2.3 is used [27]. The maximum depth of one walk is set to 4 and for each entity in the graph, the maximum number of walks is set to 10.
To train a supervised machine learning model, it is essential to have both positive and negative samples of data [26]. The positive samples are all known interactions (regulates, inhibits, or activates). The negative samples can be obtained from unknown interactions between drugs and genes. Edge embeddings for positive and negative samples are generated by adopting a node embedding fusion approach. For every possible drug-gene combination, we obtain one embedding by fusing the drug embedding and the gene embedding with the Hadamard operator, which is a strong operator for learning edge features in link prediction tasks [28]. We then add the class of the interaction (inhibits, activates, or regulates) to the resulting embedding. For edges that do not exist in the graph, we assign the label "unknown" to represent the unknown interaction class.
The prediction data for our machine learning model includes all unknown drug-gene interactions involving genes that contribute to the disease phenotype of interest. The negative samples of the training data are all unknown interactions that are not the prediction data. However, the number of negative samples significantly outweighs the number of positive samples. Including all of these negative samples could result in data imbalance and affect the performance of our model [26]. To address this issue, we decided to downsample the negative cases by randomly selecting negative samples with a sample size equal to the class in the positive set with the largest number of interactions (regulates, inhibits, or activates).
For our machine learning model, we utilized the XGBoost classifier proposed by Thafar et al. in 2021 [21]. To implement the model, we used the XGBClassifier() function from the Python package xgboost (version 1.3.3) [29]. We set the learning objective to 'multi:softmax', which allows XGBoost to optimize the likelihood of each class label and assign a probability to each possible class.
To address minor class imbalance in our positive sample, we computed sample weights using the compute_sample_weight() function from sklearn [30] version 1.1.1. These weights are then used for training the model, which provides some bias towards the minority classes during training.
To optimize the hyperparameters of our model, we conducted a randomized search on the search space presented in Table 3 using the RandomizedSearchCV() function from the model_selection module of the sklearn Python package (version 1.1.1) [30]. We set the number of parameter setting combinations to be tested (n_iter) to 20.
Given the challenge of identifying negative examples of drug-target pairs, as unlinked drugs and targets may simply represent drug-target pairs that have not been identified yet, we opted against conducting an error analysis. Instead, the model performance is assessed using the repeated stratified k-fold cross-validation technique alongside the F1-score metric.. The number of subsets for the k-fold cross validation is set to 10 and the number of repeats is set to 5. During each iteration of the process, the F1-score is calculated and averaged for each class X. Finally, the average F1-score over all k iterations and number of repeats is computed to obtain the final evaluation metric.
The best hyperparameters are used to build the final model. The confidence of each prediction is obtained using the predict_proba() function of xgboost version 1.3.3, which returns the probability of an interaction belonging to its predicted class.
For every gene in the graph that is associated with the symptom of interest, an interaction type and score is predicted for every drug in the dataset, provided that this interaction does not exist in the graph. To prioritize the most promising drug candidates for further investigation, we perform a ranking step based on the predicted positive interactions. First, we remove all predictions with a confidence score lower than 0.9 to focus only on the most confident predictions. Next, we rank the drugs based on the number of positive interactions they have with the genes associated with the symptom. This ranking approach is based on the hypothesis that drugs with more positive interactions with genes that cause a symptom are more likely to be effective in alleviating that symptom. In the case of drugs with the same number of interactions, we use the sum of prediction confidence scores as a secondary ranking criterion, with drugs having higher sums being ranked higher.
Our workflow was initially run with the terms "huntington's disease" (OMIM number '143100') and "chorea" ('HP:0002072') as seeds representing the disease and symptom fields respectively, for constructing the knowledge graph. The graph was created on 2022-06-27 and has in total 2189 nodes and 17467 edges. Figure 2 provides an overview of all entities and relationships in the graph and Table 4 shows the identifiers and descriptions of each relation between nodes in this graph. The graph includes 1352 drugs and 284 genes, resulting in a total of 383,968 edge representations, of which 1753 are known (1301 regulates, 391 inhibits, and 61 activates). This graph has 200 genes that are associated with the symptom chorea, which are the genes of interest, and there are 1077 known drug-gene edges with these 200 genes, indicating that the prediction data consists of 269,323 unknown interactions of potential interest. To deal with the imbalance between the larger number of negative samples and the comparatively smaller number of positive samples, 1301 negative samples were randomly selected, to balance the number of the largest interactions in the positive class (regulates). The best XGBoost hyperparameters can be found in Table 5, and the F1 score with these hyperparameters is 0.867. The trained XGBoost model was used to predict the classes of the unknown drug-gene interactions of interest.
Table 6 shows the predicted top ten ranked drugs that interact with genes that are associated with the phenotype chorea. We manually explored the top ranked predictions for HD that are associated with chorea. Below we present the top three candidates. Table 7 presents the two highest ranked drugs and the genes that these drugs have a positive predicted interaction with. The top predicted ranked drug is CHEMBL29097. CHEMBL29097 (synonym MK-886) is an inhibitor of 5-lipoxygenase-activating protein activity, currently in preclinical phase. It has been found that 10 microM MK-886 can abolish the biosynthetic production of cysteinyl leukotrienes (CysLTs), which is suggested to be involved in brain inflammation and neurological diseases [33].In addition to its anti-inflammatory activity, MK-886 has been shown to activate the proteasome which is known to have a causative role in HD [34]. Impaired function of the proteasome leads to the formation of intracellular aggregates in the nucleus as the proteasome cannot clear efficiently misfolded huntingtin proteins [35].
The second highest ranked drug is baicalein. Baicalein (CHEMBL8260) is a flavonoid isolated from the traditional Chinese medicinal herbal Scutellaria baicalensis Georgi, currently on Phase 2. Baicalein has known anti-inflammatory and neuroprotective efficacy in neurodegenerative disease models [36]. Rui et al. [36] studied the effects of baicalein on inflammasome-induced neuroinflammation in Parkinson's disease(PD) and found that baicalein can suppresses MPTPinduced nigral dopaminergic neuron death, glial activation, and motor dysfunction in mice by suppressing the NLRP3/caspase-1/GSDMD pathway. In addition, several studies have have demonstrated that baicalein protects neurons in animal models of Alzheimer's disease (AD) and PD by inhibiting neuroinflammation [37].
Amphotericin b (CHEMBL267345) was another drug on our list that ranked very high. Ampotericin b is an approved antifungal drug used to treat serious fungal infections. Experimental evidence shows that some antibiotics and antifungal medication have neuroprotective action through anti-aggregating activity on disease-associated proteins [38]. Although this drug has been shown to cause a delay in the formation of amyloid-𝛽, it was also found to induce toxicity [39]. However, Soler et al., [40] developed a derivative of amphotericin that has anti-aggregating action but lacks toxicity and antimicrobial activity [38].
To demonstrate the reusability of our approach, we also applied our methodology to another rare disease that currently lacks treatment; Spinocerebellar ataxia type 1 (SCA1). We used as seeds the terms "SCA1" (OMIM number '164400') and the symptom "hyporeflexia" (HP:0001265) to run our drug repurposing workflow and below we describe few of the top hits.
The first prioritized drug by our workflow was Dovitinib (CHEMBL522892) currently in phase 3. Dovitinib is a pan receptor tyrosine kinase (RTK) inhibitor that has anti-tumor activity in pre clinical models of several cancers [41]. It has been recently suggested as a candidate treatment for AD because it normalizes 𝛽 amyloid mediated transcriptional responses by targeting the CREB3L2-ATF4 heterodimerization which is responsible for the majority of the transcriptional changes occuring in AD neurons [42]. Its well tolerated safety profile and the ability to cross the blood brain barrier [42] makes it an interesting candidate for AD but also potentially for other neurodegenerative disorders that exhibit similar disease mediated changes like AD.
The second predicted drug on the list was broquinaldol (CHEMBL1394319), a small molecule that has antifungal and antibacterial activity. This is an investigational drug that was found to have activity against thyroid cancer in a high throughput screening experiment [43]. However, there is currently no evidence for being associated with neurodegenerative diseases.
Number three on the candidate drug list for SCA1 was an interesting compound, astemizole (CHEMBL296419). Astemizole is an approved second generation antihistamine drug [44] that has been found to rescue motor phenotype in a Drosophila model of PD [45].
This work presents a novel disease-drug profiling approach to identifying potential candidate compounds that could alleviate the symptoms of a rare disease. It combines two established, and widely accepted approaches, of mining rare disease-specific data from multiple public databases into a knowledge graph [11], and graph-based machine learning approaches to identifying drugtarget interactions [21]. The result is an automated workflow which makes disease-drug profile predictions targeted to specific rare diseases. We demonstrated its utility using predictions from Huntington's disease and SCA1.
The advance that this work provides to the field is the use of rare disease specific knowledge graphs. Using BioKnowlegde reviewer in a drug repurposing automated workflow enables to learn from a comprehensive view of the underlying druggable rare disease biology and pathogenesis of interest. This is advantageous over current integration methods used in rare disease research, which use information about thousands of complex disorders [2], because it leverages knowledge for precision medicine. Another advantage is that by comparing to existing solutions [8,9,10], our method harnesses heterogeneous and expressive semantic graphs for DTI prediction beyond bipartite networks. Integrating new types of entities with Semantic Web technologies enables us to represent more complex relations around drugs and targets, and it opens the possibility of learning from them and exploiting the semantics by means of methods such as RDF2Vec graph embedding methods.
Through sophisticated graph-based algorithms, we can traverse the knowledge graph to identify patterns in the data and predict potential new relationships between drugs and diseases. We demonstrated that knowledge graphs and graph machine learning used streamlined in an automated workflow gives testable DTI hypotheses for drug repurposing in the rare disease area. This can support researchers to systematically generate compound prioritization coupled with well-designed validation experiments to discover treatments for rare diseases in a timely and cost-effective manner. One limitation is that we did not integrate domain expert knowledge on disease pathobiology with patient data, which can be the basis for highly innovative drugs. While [11] gave a solution to include expert knowledge in graphs, access to patient data is a serious problem in health research. However, projects such as the EJP-RD1 are providing Semantic Web based solutions for patient data sharing.
Our results provide some interesting candidates that could potentially be of great value for the rare disease community. Some of our prioritized drugs are already associated with other neurodegenerative disorders (AD and PD) targetting neuroinflammation, which is a hallmark of the HD pathology. Other candidates (Broquinaldol and Amphotericin b) belong to the class of antifungal and antibacterial medication. These types of drugs have drawn a lot of attention and although they are mainly used to treat infections new applications are being discovered. It has been reported that antibiotics, for example Doxycycline and minocycline, have neuroprotective effects due to their anti-inflammatory properties [46].
The workflow is presented as a web application and yields promising results in a reusable and reproducible way for the rare diseases community. In the future, it could be extended and improved by the addition of experimentally validated negative interactions from reliable databases. Our current approach uses unknown drug-gene pairs as negative samples for classification. It is therefore possible that this set includes currently unknown positive interactions, which may adversely affect model training [47]. Additionally, the knowledge graph could be extended to include further information about each drug, such as side-effects and drug interactions, and additional input seeds could be obtained from the Monarch database, such as genes and metabolites associated with the particular disease, or related diseases. Moreover, integrating the predicted embeddings into the knowledge graph would enable the evaluation of performance across diverse prediction methods, offering valuable insights into model efficacy and versatility. Lastly, it is important to note that the RDF graphs are stored in turtle format at the location where the code is executed, to enable further use in code execution and analysis. The RDF graph is currently not served as a live RDF/SPARQL endpoint, which presents an opportunity for improvement in our approach. This aligns with broader efforts in the field to enhance the transparency and accessibility of machine learning outcomes.
All relations in the HD chorea KG with their identifiers and description (sources: OBO Relations Ontology [31]; GENO ontology [32]). Column 'count' shows the number of edges in each relation group.
Integrating data into knowledge graphs with state-of-the-art graph-based machine learning methods results in a novel automated drug repurposing workflow, specifically suited for rare diseases, where data tends to be sparse and distributed. The workflow produces a ranked list of candidate compounds, which serve as new hypotheses for drug treatments. The drug repurposing workflow generalises to any rare disease, but as a proof-of-concept, we focused on Huntington's disease, and a related condition, SCA1. We identified several promising candidate drugs for Huntington's Disease for the symptom chorea, demonstrating the potential of our approach. With further testing and validation, these candidates could be explored as potential treatments for the disease. The workflow is provided as a web application, in a publicly available Docker container. It is therefore accessible for researchers with no technical expertise and is a reusable and reproducible application for the rare disease community.
Table 6 The top ten ranked drugs for HD chorea. URI name https://identifiers.org/chembl:CHEMBL29097 CHEMBL29097 https://identifiers.org/chembl:CHEMBL8260 BAICALEIN https://identifiers.org/chembl:CHEMBL221137 EMBELIN https://identifiers.org/chembl:CHEMBL267345 AMPHOTERICIN B https://identifiers.org/chembl:CHEMBL308688 5,7-DIMETHOXYISOFLAVONE https://identifiers.org/chembl:CHEMBL2110660 IGMESINE https://identifiers.org/chembl:CHEMBL275809 FR-122047 https://identifiers.org/chembl:CHEMBL161343 ARACHIDONOYL GLYCINE https://identifiers.org/chembl:CHEMBL585 TRIAMTERENE https://identifiers.org/chembl:CHEMBL1269845 CHEMBL1269845