We adapted two benchmarking datasets in this article, which are the PREDICT dataset and the LRSSL dataset. The PREDICT dataset was obtained from the literature (Gottlieb et al., 2011), which was originally collected from the DrugBank database (Wishart et al., 2018) and the OMIM (Online Mendelian Inheritance in Man) database (Hamosh et al., 2005). It contains 1933 known drug–disease associations between 593 drugs and 313 diseases. The LRSSL dataset was extracted from the literature (Liang et al., 2017). It contains 3,051 drug–disease associations between 763 drugs and 681 diseases.

Without losing generality, we take the PREDICT dataset as an example to describe our method. Let p _i (i = 1, 2, … , 313) be the i-th disease and d _j (j = 1, 2, … , 593) the j-th drug. We defined the relationship between p _i and d _j as: a i j = { 0 p i  and  d j  have no association, 1 o t h e r w i s e . (1)

An adjacency matrix was then created to describe the drug–disease associations in the whole dataset which can be noted as A ∈ R ^313×593.

The basic assumption of this work is that drug-related diseases are more likely to be similar when the drugs are more similar. Three different similarity kernels of drugs and two different similarity kernels of diseases were applied. By fusing these similarity kernels, our method can identify potential associations between drugs and diseases.

The similarity kernel fusion (SKF) algorithm was applied to integrate three drug similarity kernels into a drug comprehensive similarity kernel and two disease similarity kernels into a disease comprehensive similarity kernel. Without losing generality, we take the drug similarity kernels as an example to explain this fusing process.

First, we normalize the aforementioned three drug similarity kernels (M _d,u , u = 1, 2, 3) by using the following formula: θ d , u ( i , j ) = M d , u ( i , j ) ∑ v = 1 593 M d , u ( v , j ) , (6) where θ _d,u (i, j) represents a normalized similarity corresponding to M _d,u (i, j). The matrix composed by the normalized kernel is noted as: Θ d , u = { θ d , u ( i , j ) } 593 × 593 , (7)

Second, we established a neighbor-constrained normalization kernel for each drug similarity kernel. Given the drug d _i and M _d,u , we collected the k most similar drugs as a set N _d,u (i, k). The neighborhood-constrained normalization of the M _d,u can be defined as follows: φ d , u ( i , j ) = M d , u ( i , j ) I d , u , k ( i , j ) ∑ v = 1 593 M d , u ( i , v ) I d , u , k ( i , v ) , (8) where I d , u , k ( i , j ) = { 1    d j ∈ N d , u ( i , k ) 0   d j ∉ N d , u ( i , k ) , (9)

The matrix composed by the neighborhood-constrained normalization is noted as: Φ d , u = { φ d , u ( i , j ) } 593 × 593 , (10)

Finally, the normalized kernel Θ _d,u and the neighbor-constrained normalization kernel Φ _d,u were fused by using the following iterative process: Θ d , u ( t + 1 ) = 1 2 α ( Φ d , u   ∑ r ≠ u Θ d , r ( t )   Φ d , u T ) + 1 2 ( 1 − α ) ∑ r ≠ u Θ d , r ( 0 ) , (11) where t is the number of iterations, α, a weight parameter between 0 and 1, T, the transpose operator in matrix algebra, and Θ d , r ( 0 ) = Θ d , r . (12)

After z rounds of iterative computations, we obtained the integration kernel as follows: Θ d = 1 3 ( Θ d , 1 ( z ) + Θ d , 2 ( z ) + Θ d , 3 ( z ) ) , (13)

Although useful information was extracted in the fusion process, noise is inevitable simultaneously. The information of k most similar drugs for each drug was further extracted to suppress the noise influence by defining a weight matrix; we defined an indicator function as follows: w d , k ( i , j ) = {   1      I d , 1 , k ( i , j ) = I d , 2 , k ( i , j ) = I d , 3 , k ( i , j ) = 1   0     I d , 1 , k ( i , j ) = I d , 2 , k ( i , j ) = I d , 3 , k ( i , j ) = 0 0.5                         o t h e r w i s e , (14)

Finally, the drug comprehensive similarity kernel was defined as the following formula: S d , k = { θ d ( i , j ) w d , k ( i , j ) } 593 × 593 , (15) where θ _d (i, j) is the element in the i-th row and the j-th column of the matrix Θ _d .

By applying Eqs 6–15 and employing two disease similarities, we obtained the disease comprehensive similarity kernel S _p,k . The value of k in the computing disease comprehensive similarity kernel is not necessarily the same as that of the drugs.

In this work, Laplacian regularized least squares (LapRLS) were employed to build the prediction model to uncover the potential drug–disease associations. Based on the characteristics of the model, we could predict the drug–disease associations from either the drug subspace or disease subspace.

In terms of the drug subspace, the Laplacian similarity matrix L _d was represented as follows: L d = D d − 1 / 2 ( D d − S d , k ) D d − 1 / 2 , (16) where D _d is a diagonal matrix whose diagonal elements are the sum of the corresponding row elements of the matrix S _d,k .

The drug–disease association prediction matrix F _d was calculated by minimizing the following loss function: min F d ( ‖ A − F d ‖ F 2 + β d ‖ F d T L d F d ‖ F 2 ) , (17) where A is the original drug–disease association matrix, F _d , the predicted association matrix in the drug subspace, β _d , the weighting coefficient, and ||.|| _F , the F-norm operator. The first term of the loss function aims to reduce the difference between the original matrix and prediction matrix. The second term is used to avoid the over-fitting problem.

The derivation of the optimization algorithm of LapRLS was introduced in Xia et al. (2010). We can calculate the predicted association matrix in the drug subspace by using the following formula: F d = S d , k ( S d , k + β d L d S d , k ) − 1 A . (18)

Similarly, by using Eqs 16–18 on the disease subspace, the predicted association matrix F _p in the disease subspace was calculated as follows: F p = S p , k ( S p , k + β p L p S p , k ) − 1 A , (19) where L _p is the normalized similarity matrix in the disease subspace. L _p was defined as follows: L p = D p − 1 / 2 ( D p − S p , k ) D p − 1 / 2 , (20) where D _p is the diagonal matrix of the matrix S _p,k .

With the two predicted association matrices F _d and F _p , we defined the final predicted association matrix as follows: F = λ F p + ( 1 − λ ) F d , (21) where λ ∈ (0, 1) is a weighting parameter.

We used novel association prediction and novel drug prediction to estimate the prediction performance of our method. In the novel association prediction, we applied both 5-fold and 10-fold cross-validation schemes. For the 10-fold cross validation, all drug–disease associations in the benchmarking dataset were divided into 10 nonoverlapping subsets randomly with almost the same size. While for the 5-fold cross validation, the number of the nonoverlapping subsets is five. The novel drug prediction was used to evaluate the prediction performance on new drugs. All drugs, not the drug–disease associations, were randomly divided into 10 subsets of approximately equal size. In each trial, one set was used in turn to act as the testing set, and other sets were used as the training set. In all aforementioned cross-validations, when one subset was used as the testing set, all prior knowledge of the testing set was removed before computing the association similarity kernels for drugs and diseases. The known associations corresponding to the testing set were reset to unknown. This guarantees that there is no information leak in the testing process.

All pairs of drug–disease associations were scored. Given a threshold, the drug–disease pair with a score larger than the threshold is predicted to be associated, otherwise nonassociated.

Due to the requirement of performance comparison and the performance values that are reported by the existing methods, we took a different set of performance measures in different contexts. We used the AUROC (area under receiver operating characteristics) curve and AUPR (area under precision-recall) curve to evaluate the prediction performance of our method in the context of 10-fold cross-validations. AUROC and AUPR are threshold-free metrics that are capable of measuring the overall performance of prediction models (Jiao and Du, 2016). While AUROC is powerful and popular as a main performance evaluation index, the value of AUROC may be misleading when the dataset is highly imbalanced (Saito and Rehmsmeier, 2015). It should be noticed that the value of AUPR can be used as an alternative metric to evaluate the prediction performance when the data is imbalanced (Davis and Goadrich, 2006). The value of AUPR tends to be smaller relative to the value of AUROC because of the highly imbalanced dataset (Saito and Rehmsmeier, 2015). In the context of 5-fold cross-validations, we applied several additional performance measures other than the AUROC and the AUPR. Six statistics, including sensitivity (SEN), specificity (SPE), precision (PRE), accuracy (ACC), F1-score (F1), and Matthew’s correlation coefficient (MCC), were applied in measuring the prediction performance of DDA-SKF in the context of 5-fold cross-validations. The threshold was optimized to maximize the F1-score. These performance measures were defined as follows: S E N = T P T P + F N , (22) S P E = T N T N + F P , (23) P R E = T P T P + F P , (24) A C C = T N + T P T N + F N + T P + F P , (25) F 1 = 2 P R E ⋅ S E N P R E + S E N , a n d (26) M C C = T P ⋅ T N − F P ⋅ F N ( T P + F P ) ( T P + F N ) ( T N + F P ) ( T N + F N ) , (27) where TP, TN, FP, and FN represent the number of true positives, true negatives, false positives, and false negatives in the 5-fold cross-validation, respectively.

There are three parameters in the SKF, which are the number of neighbors k, the weighting coefficient α, and the number of iterations t. There are also two parameters in the LapRLS, which are the weighting coefficients β and λ. To find the best parameter combination, we performed many trials manually with arbitrary values of parameters to optimize the AUROC value. We also performed a systematic exploration using a floating forward grid search strategy to further analyze the effects of different parameters. We finally fix the parameter values α and k as 0.2 and 15 in the drug subspace, while 0.7 and 10 in the disease subspace, respectively. The parameters β and λ were fixed as 2⁻¹⁶ and 0.4, respectively. The number of iteration t was determined using a stopping criterion. We defined the relative error E _d,u in the drug space as follows: E d , u ( t ) = ∥ Θ d , u ( t + 1 ) − Θ d , u ( t ) ∥ ∥ Θ d , u ( t ) ∥ , (22a) where u = 1, 2, 3. The relative error E _p,u in the disease subspace is defined similarly. In the disease subspace, if the E _d,u < 10⁻⁷ and t ≥ 10, the iteration will be terminated, while in the drug subspace, the iteration will be terminated if E _p,u <10⁻¹⁰ and t ≥ 10. The number of iterations was finally fixed as 10.

In this section, to evaluate the prediction performance of our model on the benchmarking datasets using both novel association prediction and the novel drug prediction. We compared our model against other state-of-the-art methods, including BNNR (Yang et al., 2019a), DisDrugPred (Xuan et al., 2019), DRRS (Luo et al., 2018), SCMFDD (Zhang et al., 2018b), MBiRW (Luo et al., 2016), and DRIMC (Zhang et al., 2020). We also compared our method against the PREDICT (Gottlieb et al., 2011), LRSSL (Liang et al., 2017), and NTSIM methods (Zhang et al., 2018a).

To eliminate the randomness in the cross-validation, we repeated each cross-validation five times with different random data partition schemes. The average value and the standard deviations of every performance measures were reported. However, it should be noted that the standard deviations were not reported by other methods in comparison.

For the novel association prediction, the evaluation results of all the methods on the benchmarking dataset are collected in Table 1. As in Table 1, DDA-SKF obtained good values in AUROC and AUPR. DDA-SKF achieved an AUROC of 0.937, which is 0.536, 5.281, 0.861, 31.601, and 2.854%, respectively, higher than that of BNNR’s 0.932, DisDrugPred’s 0.890, DRRS’s 0.929, SCMFDD’s 0.712, and MBiRW’s 0.911. DDA-SKF also has an AUPR of 0.533, which is only slightly lower than that of BNNR. Therefore, we believe that DDA-SKF has a better, or at least comparable, performance to all state-of-the art methods in novel association prediction. The integration of disease and drug similarity kernels is effective.

For the novel drug prediction, we evaluated the performance of all models in the benchmarking dataset. The results are recorded in Table 1. DDA-SKF achieved an AUROC of 0.845, which is 8.892, 1.198, 10.458, 15.280, and 5.890%, respectively, higher than that of BNNR’s 0.776, DisDrugPred’s 0.835, DRRS’s 0.765, SCMFDD’s 0.733, and MBiRW’s 0.798. DDA-SKF also obtained an AUPR of 0.270, which is only slightly lower than that of DRIMC. Therefore, we conclude that DDA-SKF outperformed most state-of-the-art methods in novel drug predictions.

In the comparison with PREDICT, LRSSL, and NTSIM methods, 5-fold cross-validation was applied rather than 10-fold cross-validations, as all three existing methods reported their performances in 5-fold cross validations. According to the results in Table 2, DDA-SKF achieved better performances in almost every comparison. DDA-SKF achieved an AUROC of 0.929 and 0.931 for the PREDICT dataset and LRSSL dataset, which are 0.869 and 3.215% higher than that of the second model NTSIM, respectively. As for AUPR, DDA-SKF achieves an AUPR of 0.497 and 0.382 on the PREDICT dataset and LRSSL dataset, which are the best values among the four. For the F1-score, DDA-SKF achieves 0.504 and 0.427 for the PREDICT dataset and LRSSL dataset, which are 25.373 and 26.331% higher than that of NTSIM, respectively. In summary, DDA-SKF achieved better performances in this comparison.

Although the prediction performances of DDA-SKF in terms of AUROC and AUPR are slightly lower than those of the best method for the novel association prediction and novel drug prediction, our method can work without enough disease similarities. This is useful when the orphan drugs, whose indication is still very limited, are considered.

In this section, the capability of our model in predicting novel drug–disease associations is tested. We performed three tests based on the PREDICT dataset in this part: one is drug-repositioning prediction, the second is orphan drug indications prediction, and the last is complex disease prediction.

To predict novel indications for all drugs, all known drug–disease associations in the benchmarking dataset were used as the training set, and the unknown drug–disease associations were regarded as the candidate set. DDA-SKF was applied to obtain scores for all candidate drug–disease associations. All candidate associations were ranked according to the prediction scores. Top-ranked candidate associations were identified as novel drug–disease associations.

We verified newly predicted associations by comparing them against the Comparative Toxicogenomics Database (CTD) (Davis et al., 2019). CTD contains curated and inferred chemical–disease relationships, which are divided into marker/mechanism and therapeutic. We only compared the therapeutic relationships for verification. For each of the 593 drugs, we collected the top-5 and the top-20 prediction results. The predictions for all drugs were listed in Supplementary Table S1 and Supplementary Table S2 as OMIM IDs in supplementary materials. In the DDA-SKF prediction results, 156 of top-5 and 377 of top-20 predictions have been confirmed by CTD.

We listed several top-5 prediction results in Table 3. Although some top-ranked predictions had not been verified, we believe that these may provide new indications for approved drugs. For example, vincristine for B chronic lymphocytic leukemia (CLL), cisplatin for acute lymphoblastic leukemia (ALL), and methotrexate for neuroblastoma. These articles provide some possible hints that consist of our predictions (Lj et al., 1997; Vilpo et al., 2000; Lau et al., 2015).

In particular, DDA-SKF can perform drug–disease association prediction without disease similarities. This may be important in studying some orphan drugs. Since the indications of orphan drugs are usually very limited, the similarities between diseases may be misleading. Therefore, by excluding the disease similarities from the model, DDA-SKF should have more potential in finding new indications for orphan drugs, which may eventually decrease the cost of orphan drug utilization.

We selected several orphan drugs from Orphanet (http://www.orpha.net), including celecoxib, methotrexate, and doxorubicin. For each orphan drug, all known drug–disease associations were removed. The orphan drug would have no association information during the process of prediction. We fed the drug similarity kernels with orphan drug information and the aforementioned association matrix to our model. The results of top-20 predictions for each orphan drug were recorded in Supplementary Table S3 in supplementary materials. The top-5 predicted results are summarized in Table 4. These successful prediction instances further confirm that DDA-SKF has the potential to predict novel indications for orphan drugs without disease similarity information.

Last, several commonly studied complex diseases were chosen to evaluate DDA-SKF. For instance, Alzheimer’s disease (AD), Parkinson’s disease (PD), and amyotrophic lateral sclerosis (ALS) were considered. We listed the top-5 predicted results in Table 5 along with the literature evidences. For each disease, the uncovered associations will be ranked based on the prediction scores, and the top-ranked unknown associations were identified as the prediction results. The AD and PD results were obtained on the PREDICT dataset, while the ALS results were obtained using the LRSSL dataset. These results indicate that DDA-SKF has the potential to uncover new drugs for these complex diseases.

In this work, different drug similarity kernels were integrated as a drug comprehensive similarity kernel in the drug subspace. Also, different disease similarity kernels were integrated as a disease comprehensive similarity kernel. To demonstrate the effectiveness in integrating multiple similarity kernels, we calculated the values of AUROC and AUPR with every single similarity kernel on the PREDICT dataset. The results are shown in Figure 2 and Figure 3. We also compared the effectiveness of the multiple similarity kernels fusion and all single similarity kernels. The results are recorded in Table 6. We can see that the comprehensive kernel produced higher performances than every single kernel.