SemaTyP: a knowledge graph based literature mining method for drug discovery

Background Drug discovery is the process through which potential new medicines are identified. High-throughput screening and computer-aided drug discovery/design are the two main drug discovery methods for now, which have successfully discovered a series of drugs. However, development of new drugs is still an extremely time-consuming and expensive process. Biomedical literature contains important clues for the identification of potential treatments. It could support experts in biomedicine on their way towards new discoveries. Methods Here, we propose a biomedical knowledge graph-based drug discovery method called SemaTyP, which discovers candidate drugs for diseases by mining published biomedical literature. We first construct a biomedical knowledge graph with the relations extracted from biomedical abstracts, then a logistic regression model is trained by learning the semantic types of paths of known drug therapies’ existing in the biomedical knowledge graph, finally the learned model is used to discover drug therapies for new diseases. Results The experimental results show that our method could not only effectively discover new drug therapies for new diseases, but also could provide the potential mechanism of action of the candidate drugs. Conclusions In this paper we propose a novel knowledge graph based literature mining method for drug discovery. It could be a supplementary method for current drug discovery methods. Electronic supplementary material The online version of this article (10.1186/s12859-018-2167-5) contains supplementary material, which is available to authorized users.


Background
Drug discovery is the process through which potential new medicines are identified. High-throughput screening (HTS) and computer-aided drug discovery/design (CADD) are the two main drug discovery methods for now [1]. Despite advances in technology and understanding of biological systems, drug discovery is still a lengthy and expensive process with low rate of new therapeutic discovery [2,3]. Developing a new drug is estimated to take 14 years and cost approximately $1.8 billion [4]. In contrast, Literature-Based Discovery (LBD) is a safe and low-cost approach to identify new drugs for indications. LBD seeks to discover new relationships in existing knowledge from unrelated literatures [5]. Drugs are often discovered on the serendipitous observation that a drug effect may be therapeutically useful if it induces a desired *Correspondence: yangzh@dlut.edu.cn 1 College of Computer Science and Technology, Dalian University of Technology, Hongling Road, 116023 Dalian, China Full list of author information is available at the end of the article effect or counters a disease phenotype [6]. For instance, Don R. Swanson (1924Swanson ( -2012 proposed fish oil as a new treatment for Raynaud's disease in 1986 after noting the association "high blood viscosity is observed among Raynaud's Syndrome sufferers" in some biomedical articles and another association "dietary fish oil lowers blood viscosity" in other articles [7]. This hypothesis was verified in medical experiments two years later. Basic LBD techniques search for a set of intermediate terms that frequently co-occur with a source term and a target term [5]. As shown in the above example, "blood viscosity" is the intermediate term in associating the "dietary fish oil" with the "Raynaud's Syndrome". In addition, more sophisticated LBD methods first employ natural language processing (NLP) techniques to extract relations between entities from biomedical literature. Then novel discoveries could be analyzed from the extracted relations [8]. For example, Hristovski et al. used SemRep to extract rela-tions among entities from biomedical literature [9]. These extracted relations could then be used for inferring novel relationships in literatures [8]. More recently, a number of recent LBD methods have explored methods that utilize certain graph data structures. For example, Cameron et al. introduced a graph-based method that automatically finds clusters of contextually similar paths in a semantic graph [10,11]. These clusters are used to elucidate the latent associations between disjoint concepts in the literatures. These existing LBD methods have several limitations. The main issue of of term co-occurrence approach is that the extracted relationships lack logical explanations [12]. NLP-based methods strongly depends on the availability of domain-specific NLP tools [13]. Graph-based methods don't consider the different semantic types of nodes in the graph. Most importantly, all existing methods have not exploited all available published biomedical literature for drug discovery. They only focus on part of the abstracts related to disease of interest. This could lead to missing the valuable informations existing in the filtered literature.
In this paper, we propose a biomedical knowledge graph based inference method to discover drug therapies from literature. Knowledge graphs (KGs) are collections of relational facts, which have proven to be sources of valuable information that have become important for various applications [14]. The famous knowledge graphs include Freebase [15], DBpedia [16], Nell [17] and YAGO [18], etc. Here, we first construct a biomedical knowledge graph called SemKG with relations extracted from PubMed abstracts. Then based on SemKG, a drug discovery method called SemaTyP (Semantic Type Path) is introduced to exploit the semantic types of paths to discover drug therapies. The experimental results show that our method could not only discover new candidate drugs for new diseases, but also could provide the mechanism of action of the candidate drugs. To summarize, the contributions of the paper is: First, we introduced a biomedical knowledge graph -SemKG -which is constructed by integrating information extracted from PubMed abstracts. Second, this is the first method that discovers candidate drugs by using biomedical knowledge graph. Our method could be a supplementary method for current drug discovery methods, which could improve the successfulness in discovering new medicine for recently incurable diseases.

Materials and tools
The biomedical knowledge graph used in this study is constructed based on the predications (subject-relationobject triples) extracted from PubMed abstracts by Sem-Rep. In this section, the datasets and tools used in this study are briefly introduced.

PubMed
PubMed is a free search engine accessing primarily the MEDLINE database of references and abstracts on life sciences and biomedical topics. It provides now access to more than 26 million citations, adding thousands of records daily [19].

UMLS semantic network
The Unified Medical Language System (UMLS) semantic network consists of 133 semantic types and 54 relationships that exist between the semantic types. In this paper, the abbreviations are adopted to represent the semantic types. For example, 'podg' represents 'Patient or Disabled Group' and 'topp' is 'Therapeutic or Preventive Procedure' .

Metamap
MetaMap is a widely available program providing access from biomedical text to the concepts in the unified medical language system (UMLS) Metathesaurus [20]. It could be applied for biomedical name entity recognition, word sense disambiguation (WSD) and other natural language processing tasks [21].

SemRep
SemRep is a relation extraction tool which first uses MetaMap to map noun phrases to UMLS concepts [22] then extracts semantic predications from biomedical free text [23]. For example, from the sentence "We used hemofiltration to treat a patient with digoxin overdose that was complicated by refractory hyperkalemia", Sem-Rep extracts four predications: On the right of symbol '|' is the abbreviation of entity's semantic type (black bold).

Construction of SemKG
Knowledge graph is a multi-relational graph composed of entities as nodes and relations as different types of edges. In this work, we constructed a biomedical knowledge graph, called SemKG, with the predications which are extracted from PubMed abstracts by SemRep. In the SemKG, let E = {e 1 , e 2 , . . . , e N } denote the set of n entities, R = {r 1 , r 2 , . . . , r M } denote the set of relations between entities and T = {t 1 , t 2 , . . . , t K } denote semantic type of entities. The elements of R and T are all from the UMLS semantic network. The edge between entities e i and e j is weighted by the number of predications that have been extracted. Besides, the attribute of edge includes the abstracts' PubMed ID (pmid) from where the predications are extracted. A prototype example of the SemKG is illustrated in Fig. 1. Figure 2 is an illustration of an edge of the SemKG, it shows that there are three different relations between "hydrocortisone" and "sleep, slow wave" which are extracted from four abstracts (pmid 15714228, 3657191, 3725299 and 4495256). The relation "AFFECTS" is extracted from two abstracts (pmid 15714228 and 3657191) simultaneously. Figure 2 shows the same entity could be assigned with different semantic types. For example, the "hydrocortisone" is a kind of "hormone" (horm) in the predications extracted from the two abstracts (pmid 15714228 and 3657191) and it also could be "Pharmacologic Substance" (phsu) in other predications (pmid 4495256).

SemaTyP method Path exploration
Given a knowledge graph KG, a path π is defined as a sequence of predications e 0 r 0 e 1 r 1 . . . r −1 e , where is the length of path π. For a gold standard drug i − target i − disease i case, which provides information about targeted disease i and the corresponding drug i directed at the target i . SemaTyP first constructs training data by obtaining all paths π = ρ(drug i → disease i ; target i , ), which encodes a path of length reaching node disease i from source node drug i and crossing node target i . Then p = π 1 , π 2 , π 3 , π 4 . . . is the set of all length paths. All paths in ¶ = {p 2 , p 4 , p 5 , . . . , p } are considered as positive training data. The minimum length of path in ¶ is 2, which represents the path drug i − target i − disease i . Similarly, the corresponding negative training data is obtained from a set of false cases drug j − target j − disease j .

SemaTyP feature selection
For each path π i , a training data (x i , y i ) is constructed, where x i is a vector of semantic types and y i is a boolean variable indicating whether π i is a positive case. The process of constructing x i for π i is as follows: The symbol c denotes component of path π i . (c) constructs an occurrence number vector of semantic types for c. T_E =[ te 1 , te 2 , . . . , te K ] is a vector of semantic type of entities, the entry of vector is the number of occurrence of corresponding semantic type. Similarly, T_R = [ tr 1 , tr 2 , . . . , tr M ] denotes a vector of relations and the entry is the number occurrence of corresponding relation. The symbol is concatenation of two vectors. For π i , a length of K * ( + 1) + M * training vector is constructed, where K is the length of T_E and M is the length of T_R. Figure 3 shows an prototype example of constructing one training data. As shown in Fig. 3, the T_E collects the number of occurrence of all semantic types of corresponding entity, and the T_R collects the number of occurrence of all relations between its two entities. For the drug − entity 1 − target − entity 2 − disease case, a length of (K * 5 + M * 4) vector is constructed.
For other path π m i (m < ), it is extended to length by reduplicating entity target. For example π m i = e 0 r 0 t . . . r m−1 e m is converted to e 0 r 0 tr 0 tr 0 t . . . r −1 e , where t denotes target in this example.

Training model
Given a set of training vectors, a logistic regression model is trained to predict conditional probability P(y|x; θ). We treat the number of semantic types as features for the logistic regression model.
Where the θ i are appropriate weights for the number of semantic types. The parameter vector θ is estimated by Fig. 1 The prototype example of SemKG. The symbol e, r and t represent entity, relation and the type of the entity, respectively. no is the number of occurrences and pmid is PubMed ID Where λ 2 controls L 2 -regularization to prevent overfitting. o i (θ) is the per-instance weighted conditional loglikelihood given by Where p i is the predicted probability The trained logistic regression model is used for discovering candidate drugs for each disease.

Implementation of SemaTyP
To evaluate a potential treatment case drug candidate − target candidate − disease, first a set of paths ¶ candidate = {ρ(drug candidate → disease; target candidate , 2... )} are obtained by aforementioned method. Then the score of the drug candidate for disease is: where χ(π i ) is the feature selection process for π i and n is the number of paths in ¶ candidate .
Since the treatment of the interested disease is unknown, all drugs or chemicals could be one of candidate drugs for the disease. Then all combinations of the drugs and targets are constructed to be hypothetical treatments. Finally, the candidate drugs are ranked by their score.

Baseline method
Random walk algorithm (RWA) generates finite Markov chains, which can be viewed as random walk on a directed graph [24]. RWA has been employed to resolve a series of problems due to the wide applicability of the algorithm [25]. Here, we compare our method with RWA and other two RWA-based methods, which are considered as the baseline methods.

Basic notions of RWA
Let G = (V , E) be a directed graph with n nodes and m edges. A random walk on G is considered as follows: RWA starts at a node υ 0 ; if t-th step is node υ t , RWA moves to the neighbor of υ t with probability 1/deg(υ t ). The output of a random walk is a Markov chain (υ t : t = 0, 1, ...). We denote by P t the distribution of υ t : We denote by M = (p i,j ) i,j∈ϒ the matrix of transition probabilities of this Markov chain. So Let A G be the adjacency matrix of G and let D denote the diagonal matrix with (D) ii =1/deg(i), then M = DA G . The rule of the walk can be expressed by the equation the distribution of the t-th point is viewed as a vector in R V , and hence It follows that the probability p t ij that, starting at i, the algorithm reaches j in t steps is given by the ij-entry of matrix M t .

Two RWA-based competing methods
In addition to RWA method, we compared our method with two state-of-the-art drug repositioning methods which are NRWRH [26] and TP-NRWRH [27]. NRWRH is a network-based random walk algorithm with restart on heterogeneous network. TP-NRWRH is a two-pass random walk with restart on the drug-disease heterogeneous network. Both of these two methods focus on predicting new targets for a drug of interest.

Implementation for drug discovery
To evaluate a potential drug candidate for treating disease i , the starting node υ 0 of RWA-based methods is set to drug candidate . Figure 4 illustrates an example of evaluating "chlorpromazine" to be the treatment of "cardiachypertrophy". Figure 4a is a weighted semantic graph with 7 nodes and 9 edges. Figure 4b shows the results of RWA with starting node "chlorpromazine". It shows that when the step of RWA is 1, "chlorpromazine" can't reach "cardiachypertrophy", then the score of "chlorpromazine" of step_1 RWA is 0. Similarly, the score of "chlorpromazine" for treating "cardiachypertrophy" is 0.697 when the step is 4. For each disease i , RWA scores all candidate drugs of the disease. After that the candidate drugs can be ranked by their scores.

Results
In this section, we firstly introduce the details of the SemKG and the training data constructed in our experiment. Then, several metrics are introduced to measure the performance of SemaTyP. After that, case studies are conducted to confirm the ability of SemaTyP to find potential drugs for indications.

The SemKG and training data The SemKG
The predications extracted from all abstracts in PubMed (before June 1, 2013) are used to construct the SemKG. Since the performance of SemRep is not perfect: its precision, recall, and F-score are 0.73, 0.55, and 0.63, respectively [28],and the low precision (73%) means many false semantic associations will be returned [12]. We filtered out all the predications that are only extracted once in order to ensure the quality and accuracy of the extracted predications. Table 1 shows the details about the SemKG. Figure 5 is the distribution of top 20 types of entities in the SemKG. For example, the first five types in SemKG are dysn (Disease or Syndrome), podg (Patient or Disabled Group), bpoc (Body Part, Organ, or Organ Component), aapp (Amino Acid, Peptide, or Protein) and topp (Therapeutic or Preventive Procedure).

Training set
In this work, 7144 drug − target − disease are extracted from Therapeutic Target Database (TTD) as true cases (Additional file 1). The is set to 4, K is 133 and M is 52. Based on the aforementioned construction of training data, 19,230 positive data are obtained. Each data is a length of 873 (133*5+52*4) vector. On the b a Fig. 4 Random Walk Algorithm for drug discovery other side, for each drug − target − disease, we random replaced the drug, target and disease with other drug, target and disease. If the new triplet doesn't exist in TTD, then it is considered as a false example, which is denoted as drug − target − disease . Similarly, 19,230 negative training data is obtained from false cases.

Evaluation metrics
To systematically evaluate the performance of our method, we conduct ten-fold cross validation and drug rediscovery test.
In the ten-fold cross validation, all training data are randomly divided into ten subsets with equal size. In each cross validation trial, one subset is taken in turn as the test set, while the remaining nine subsets constitute the training set. After performing prediction, each test case is given a predicted score. According to the final predicted scores, the case is assigned a boolean label indicating whether it is a positive case. In this study, the Precision, Recall and F-score are adopted to measure the performance of SemaTyP method. In our study, drug rediscovery test is performed to evaluate the effectiveness of the SemaTyP when predicting potential drugs for new diseases. For each disease of interest, a list of candidate drugs are constructed to be scored by SemaTyP. Considering the fact that the predicted topranked results are more important in practice, we measure the performance of our method in terms of the top-ranked results, i.e. the mean ranking of true therapies and the proportion of correct therapies ranked in the top 10. Usually, it is regarded as more effective if the method can rank more true therapies in top portions.

Ten-fold cross validation
We explored a range of values for the L 2 -regularization parameters λ 2 using cross validation on the training data. Figure 6 shows that parameter λ 2 ranging from 0.0001 to 100 has little effect on the prediction performance and a small amount of L 2 -regularization can slightly improve performance of SemaTyP. In this study, we set the parameter λ 2 to 1.0. The precision, recall and F-score are 0.907, 0.879 and 0.892, respectively. In addition, we also compared the L 2 penalty with Lasso (L 1 ) regularization [29]. As same to L 2 regularization, the parameter λ 1 of Lasso regularization ranges from 0.0001 to 100. Table 2 shows the comparison results of L 1 and L 2 regularization. The results show that the model achieves higher performance with L 2 regularization. This is because L 1 regularization is often used for feature selection [30] when the number of potentially relevant features is very large. However, in this work the number of features we selected is not large (873).
We vary the number of training data to see how training data size affects the quality of the model. Figure 7 shows that our method benefits from more training data, and it is especially evident when more than half of all the data are used. Figure 7 shows that the increase in training data significantly improves the performance of SemaTyP when less than 50% training data are used. After that, the increase in training data slightly improves the performance of the method.
Additionally, we vary the settings of to see how pathway length affects the results. The was set to 2, 3 and 4, respectively. Table 3 shows the results of our model with different . It shows that when the is 2, 32 training data was obtained by aforementioned method. It means there are only 32 drugs connect to their indications by directly crossing corresponding targets. We didn't train the model with the training data, since 32 training data is not enough for training a machine learning model. As shown in Table 3, 1742 data was obtained when is 3. The performance of our model trained by the 1742 data is shown in Table 3. Table 3 shows that the performance of our model with equals 4 is better than equals 3 as expected. As Fig. 7 shows that the increase in training data could significantly improve the performance of our model. When is 3, the size of training data is 9.06% of the training data obtained by equals 4.
In this work, the is set to a value less than 5, it's because: 1) Although more training data could be obtained when exceeds 4, Fig. 7 shows that when the training data exceeds certain size, the performance of our method is relatively stable. 2) As increases, longer paths starting from a drug to a disease are obtained. However, Fig. 6 The performance of SemaTyP more entities in a drug-disease path might reduce the quality of training data. Therefore, in this work, we set the to 4.

Drug rediscovery test
To evaluate the capability of our method in discovering potential drugs for new diseases, we conduct the drug rediscovery test. In this test, 360 drug − disease relationships (Additional file 2) are selected from TTD as gold standard to form test set. Each disease i in test set has one known associated drug i , but the drug mechanism of action is not clear. For each disease i we randomly selected other 99 drugs or chemicals from TTD as candidate drugs for this disease. We report the mean of those predicted ranks of drug i and the hits@10, i.e. the proportion of known drugs ranked in the top 10. If the known drug of a disease is not rediscovered, then the score for the drug is set to -1 and the ranking number is 101. Specifically, for disease i and candidate drug j , 5,785 drug j −target candidate −disease i are constructed. This is due to that the targets of disease i are unknown, then each target (protein) in TTD could be the target candidate of disease i . For disease i , the comparison methods also scores and ranks all 100 candidate drugs. The step of RWA is set from 1 to 10. The NRWRH and TP-NRWRH methods are configured to their recommended settings in their papers. Table 4 shows the results and the "Not found" column is the number of known drugs which are not found by the method. As we can see from Table 4, there are 262 gold standard drugs are not discovered by RWA_1 (random walk algorithm and the step is set to 1). It means that only 98 (360-262) drugs directly connect to the disease in the SemKG. The "Not found" number decreases when the step number of RWA increases. Table 4 shows that all drugs could be found by RWA when step length exceeds 3. It's because all drugs could be connected to the disease in the SemKG through a semantic path whose length is greater than 3. Table 4 shows that there are 19 and 17 drugs are not found by NRWRH and TP-NRWRH, respectively. Although the step of the two RWA-based methods is 3, NRWRH and TP-NRWRH are both random walk algorithm with restart. This could result in the diseases fail to reach the appropriate drugs within 3 steps.
For the "Mean ranking" column, the worst result is obtained by RWA_1 (72.28), it is due to there are 262 known drugs are not found by RWA_1. As the step length of RWA increases to 2 the meaning ranking decreases to 26.59, it's because more drugs could be discovered by RWA_2 than RWA_1. But when the step of RWA continues to grow, the mean ranking improves. It's because  although all known drugs could be discovered when the step of RWA exceeds 3, more other candidate drugs also could be found. The more discovered candidate drugs could lead the ranking of true drugs decreasing. Table 4 shows that NRWRH and TP-NRWRH achieve better performance than RWA method, it's because: 1) The best performance of RWA on "Mean ranking" is achieved when the step is 3, and the step of NRWRH and TP-NRWRH is 3. 2) NRWRH and TP-NRWRH methods integrate biomedical background knowledge to choose next step rather than randomly step to next node. For "Hits@10", the value of "Hits@10" decreases when the step of RWA increases. For RWA method, Table 4 shows that RWA_3 and RWA_4 achieve the best performance: 1) almost all drugs could be discovered and 2) the "Mean ranking" value is relatively small and the "Hits@10" is relatively large. In addition, Table 4 shows NRWRH and TP-NRWRH achieve better performance than RWA method. We could see from Table 4, our method achieves the best performance in both tests. The "Mean ranking" of our method is 26.31 and the "Hits@10" is 48.61%. The reasons of our method outperform others are: 1) we could know from Table 4 that when the step of RWA is 3 or 4 the RWA achieves the best performance. Our method could cover all the paths whose length is 2 to 4. 2) Our method Bold values denote the best scores corresponding to specific metric scores the semantic path based on the distribution of their semantic types other than only based on the structure of the SemKG.

Case study
We conduct 12 case studies to demonstrate the efficacy of our methods (Table 5). For each disease, SemaTyP can predict the potential drugs and the corresponding targets simultaneously. For example, TTD has reported that testosterone and ap22408 are known drugs for osteoporosis. These two drugs are ranked 1st and 3rd as potential drugs for osteoporosis by our method. What's more, SemaTyP also provides corresponding targets for the drugs, which have not been discovered for now. For instance, terikalant is predicated to treat cardiac arrhythmia by acting on actin. Aspirin, is predicted to treat cardiovascular disease by acting lymphoid cell, etc. These prediction instances further confirm that SemaTyP not only has the potential to predict novel drugs for disease, but also could provide potential mechanism of action for the drugs.

Discussion
To the best of our knowledge, this is the first method that employs knowledge graph for solving LBD tasks. This paper showed that use of implicit semantic types to find drugs from literature can be effective for LBD. Our overall approach however, has several limitations. The first limitation is the construction of knowledge graph -SemKG -relies heavily on effective NLP tools. On one hand, the accuracy of MetaMap reduces in the presence of ambiguity, which leads its inability to resolve word sense disambiguation [20]. On the other hand, although the isolated predications are filtered out in order to improve the quality of the SemKG, there are still considerable number of false predications existing in the knowledge graph, which could lead to our method inferring lowerquality results. In addition, in the process of constructing SemKG, more than half the initial predications are filtered out, which might lead to possible selection biases in the step. The second limitation is SemaTyP relies on the semantic types of nodes and edges to infer associations, hence our method is effective only when the required ontology are easily available. Another limitation is SemaTyP needs to obtain all paths between candidate drug and disease. When the scale of knowledge graph is large, it's difficult for our method to obtain long paths. These and other limitations suggest the next steps in this research. In future, high-quality NLP tools need to be developed to improve the quality of SemKG. Additionally, another representation of nodes and edges in SemKGgraph embedding -could be useful for our method to obtain long paths.

Conclusion
In this work, we have presented a novel method named SemaTyP uncovering the potential associations between drugs (chemicals) and diseases from literature. We first constructed a biomedical knowledge graph by integrating informations extracted from PubMed biomedical literature. Then based on the knowledge graph, we devised a novel model to discover potential drugs and corresponding targets. Finally, we test our method on two different tests. The experimental results show that our method can effectively discover drugs for diseases from literature. Our method has potential to accelerate drug development and benefit the field of target identification.

Additional files
Additional file 1: Supplementary Data 1. The gold standard drug-target-disease cases used in this work. The 7144 drug-target-disease cases which are extracted from Therapeutic Target Database (TTD) as true cases for constructing training data. (TXT 466 kb) Additional file 2: Supplementary Data 2. The gold standard drug-disease cases extracted from TTD. There are 360 drug-disease relationships are selected from TTD as gold standard to form test data for drug rediscovery test. Each disease i in test set has one known associated drug i , but the drug mechanism of action is unclear. (TXT 10 kb) Abbreviations CADD: Computer-aided drug discovery/design; HTS: High-throughput screening; LBD: Literature-based discovery; NLP: Natural language processing; TTD: Therapeutic target database