Since an RBM can be efficiently applied to learn the distribution of multidimensional networks and reconstruct their inputs, we developed an RBM-based model, DDIT, to predict different phenotypic types of DDAs. Figure 1 shows the flowchart of DDIT.

Drug indications are gold standard dataset from PREDICT (Gottlieb et al., 2011). The data for drug side effects is obtained from SIDER (Kuhn et al., 2016). The data for drug contraindications are from MED-RT (https://ncit.nci.nih.gov/ncitbrowser/pages/vocabulary.jsf?dictionary=MED-RT), produced by The Veterans Health Administration (VHA). For this, we downloaded the archive content Core_MEDRT_2019.11.04_XML.zip (https://evs.nci.nih.gov/ftp1/MED-RT/Archive/) and then extracted the relationship of “CI-with,” which describes co-morbid contraindication of a drug (see Supplementary Figure S1). In total, we collected 2,816 drug–indication pairs, 132,150 drug–SE pairs, and 10,443 drug–contraindication pairs.

As the names of drugs and diseases in different datasets often vary in their vocabulary, this required consideration and adjustment for standardization. For example, the drug names in the indication dataset, side effect dataset, and contraindication dataset were from DrugBank (Wishart et al., 2007), ATC (Miller and Britt, 1995), and RxNorm (Nelson et al., 2011), respectively, while the corresponding disease names in these three datasets were from OMIM (Hamosh et al., 2005), UMLS (Bodenreider, 2004), and MeSH (Leydesdorff et al., 2016), respectively. To unify the drug and disease names, we mapped their names to Unified Medical Language System (UMLS) ontology. We accessed UMLS Knowledge Sources Metathesaurus 2019AB using Rest API for Java (https://github.com/HHS/uts-rest-api). All the data were completely mapped to UMLS ontology, and no other tools were used to unify these drug and disease terms.

During the data collection phase, we collected 2,816 drug–indication pairs, 132,150 drug–SE pairs, and 10,443 drug–contraindication pairs. As these data sets were unbalanced, to balance the data, we first selected a disease subset S, then we randomly select 2,816 (2,816 is the lowest count among three types of DDAs) associations from drug–SEs and 2,816 associations from drug–contraindications with disease in the subset S. To obtain S, we have selected the diseases meeting to one of following two criteria: (1) diseases included in the known drug indications and (2) diseases shared by two different types of clinical phenotypes. Finally, we obtained the 2,816 data points for each type of associations.

An RBM is an undirected graphic model (Salakhutdinov et al., 2007) that can be used to learn probability distributions over input data using a layer of binary hidden units. As shown in Figure 2, an RBM consists of a layer of visible units (v) and a layer of hidden units (h). Each visible unit is connected to all hidden units and has no intralayer connections between any pairs of visible units or any pairs of hidden units. The state of each unit possesses a binary value.

In our study, we built an RBM model for each drug. In other words, for a drug, we adopted a two-layer RBM with diseases as visible layer and 400 hidden units as hidden layer. The hidden layer represented the hidden factor, and it cannot be observed. Each RBM model for a drug only had diseases related to the drug as visible units. Thus, different drugs had different RBM models. However, different RBMs of drugs shared the connection weight between each visible disease unit and hidden unit pairs. We assumed that for each drug, the RBM model had n visible units, m hidden units, and l association types encoded in a visible unit. In our context, each visible unit represented a disease. Therefore, we let binary vector v i = ( v i 1 , ⋯ , v i k , ⋯ , v i l ) denote the state of the ith visible unit, where visible variables v i k = 1 if the kth type of DDA is observed in the input data, and v i k = 0 otherwise. For example, for indication, the binary vector is v _i = (1,0,0), and for both side effect and contraindication, the binary vector is v _i = (0,1,1). With a 3-bit vector, it will be able to distinguish the three types of DDA at the same time. For each hidden unit, the state of jth hidden unit is expressed as h j ,   j = 1,2 , ⋯ , m . Let W i j k denote the weight of the connection between visible variable v i k and hidden variable h j , and it is shared by different RBMs of drugs. The vector v = ( v 1 , v 2 , ⋯ , v n ) denotes the input layer, while h = ( h 1 , h 2 , ⋯ , h m ) denotes the hidden layer. Figure 3 shows the modeling of four drugs and two diseases. Through the CD algorithm, the model can be effectively trained (Hinton, 2002). An RBM can learn the distribution of multidimensional networks well and reconstruct the input. This will predict the DDAs that are not yet known in the input. RBM details are supplied as Supplementary Material (Supplementary Text S1).

As the verified drug–disease associations provide more reliable information than those that are as yet unknown, we further introduced a conditional RBM to incorporate this additional information to affect the states of hidden units (Salakhutdinov et al., 2007). In this, we let r i = ( r 1 , ⋯ , r j , ⋯ , r m ) be a binary vector, in which r i = 1 if disease j has association with the current drug i, and r i = 0 otherwise. Details about the conditional RBMs are supplied as Supplementary Material (Supplementary Text S2).

The algorithm is implemented in Python. Through grid search, we determined the best parameters: (“m”: 400, “learning_rate”: 0.5, “epochs”: 300).

According to the computational method of similarity, KNN (Guo et al., 2003) was divided into drug- and disease-based KNN. For drug-based KNN, the hypothesis was that similar drugs should have similar effects on the same disease. For any given drug, we identified the other top k drugs that were most similar to it and then calculated its phenotypes by averaging the phenotypes of its neighbor drugs. The drug similarity is calculated by Jaccard Index based on the drug-related disease profile with the formula defined in Eq. 1 as follows: J = M 11 M 01 + M 10 + M 11 (1)

The approach for disease- and drug-based KNN was similar where the similarity between two diseases was calculated by disease-related drug profile similarity.

For KNN, we further optimized for k systematically. For drug-based KNN, we searched k ranging from (1, 10) and choose a k = 6 for the best AUC and AUPR with 0.736 and 0.824, respectively, in the indication prediction. For disease-based KNN, we searched k ranging from (1, 10) and choose k = 5 for the best AUC and AUPR with 0.777 and 0.916, respectively, in indication prediction.

We adopted Random Forest as a classifier for comparative analysis. For each drug–disease pair, we combined the drug-related disease profile and disease-related drug profile together as the feature vector to train the Random Forest prediction model. For three association types, we constructed three random forest models, respectively. Each model was a problem of binary classification. We implemented this algorithm by using the “RandomForestClassifier” function in the sklearn package with default parameters.

We further carried out XGBoost on our dataset. We modeled the same data as the input of random forest. We adopted the Scikit–Learn Wrapper interface for XGBoost to create the XGBoost model. We optimized the hyper parameters using “GridSearchCV” function and obtained the best parameters (“gamma”: 0.25, “learning_rate”: 0.1, “max_depth”: 5, “reg_lamba”: 0, “scale_pos_weight”: 1). We then trained the model with the best parameters.

We used the 10-fold cross-validation method to evaluate the model. In this method, the data set was randomly divided into 10 sub-parts with nine of them used as the training set in turn and the remaining one being the test set.

To access how trained models can be generalized into groups of drugs that the models have never trained on before, we further make a leave-one-drug-class-out cross-validation (Yao et al., 2019). We first mapped the UMLS concept to the ATC code, then divided the drugs into 15 classes by the ATC code (See Supplementary Table S1). The data sets were divided into 15 parts according to drug classes. Fourteen of them were used as the training set in turn, the remaining one being the test set.

To access how trained models can be generalized into groups of diseases that the models have never trained on before, we further made a leave-one-disease-class-out cross-validation. We divided the diseases into 23 MSH classes (see Supplementary Table S2). Twenty-two of them were used as the training set in turn, the remaining one being the test set.

The web server of DDIT was built using modern frontend–backend architecture with three main components: front end, backend server, and a relational database containing the information of drugs, diseases, and their phenotypic associations. The database was built using Mysql 5.6. The backend was implemented in Java using SSM (Spring + SpringMVC + Mybatis) as a framework and provided the REST API (Sohan et al., 2017). The front-end was built using React (Gackenheimer, 2015) and several other libraries. The backend was deployed in Apache Tomcat (Vukotic and Goodwill, 2011), while the front-end was deployed in Nginx (Nedelcu, 2013). This architecture provided for the easy maintenance of each module.

Receiver operator characteristic (ROC) and precision–recall (PR) curves were used as evaluation metric for predictive performance. We compared DDIT integrating three phenotypic types of DDAs with one single phenotypic type. Both in terms of ROC (Figure 4A) and PR curves (Figure 4D), DDIT with integrated three phenotypic types performed >0.079 better than single indication data in indication prediction. Similarly, in prediction of side effect and contraindication, area under ROC curve (AUC) and area under the PR curve (AUPR) had improved by >0.086 (Figures 4B, E) and >0.102 (Figures 4C,F) respectively. This suggests that data integrating multiple clinical phenotypic types provide more information than single analysis and simultaneously improve prediction performance.

We then evaluated the performance by comparing DDIT with the drug-based KNN, disease-based KNN, Random Forest, and XGBoost (see Materials and Methods). As shown in Figure 5, DDIT represented improvement by at least 0.217 in AUC and 0.072 in AUPR compared with the other four methods. The AUC and AUPR for leave-one-drug-class-out is also shown in Supplementary Figure S2. Here, DDIT represented improvement of at least 0.135 in AUC and 0.075 in AUPR, compared to the other four methods.

We then searched for an external validation dataset from CTD (Davis et al., 2021), DrugBank (Davis et al., 2021), and DynaMed (https://www.dynamed.com/) to evaluate the prediction results. We collected the novel DDAs from the CTD database. These had not been used for building DDIT, but the drugs and diseases of these novel DDAs were contained in our modeling datasets. Table 1 shows the top 10 predicted indications. Seven of these 10 predictions could be found in CTD, DrugBank, or DynaMed databases. The remaining three predictions may represent candidate drugs for new indications. For example, DDIT predicted that Bleomycin is indicated for small cell lung cancer. This is conceivable as a true positive, since Bleomycin, as recorded in DrugBank, is a drug for the treatment of malignant neoplasms and operates by inhibiting DNA synthesis. That Buspirone is a candidate drug for obsessive–compulsive disorder may also be true positive, as it is labeled as indicated for anti-anxiety in DynaMed, which is a symptom of obsessive–compulsive disorder.

For the prediction of side effect, Table 2 shows the top 10 predicted side effects. Nine of 10 can be found in DynaMed. DDIT inferred that dermatitis was a side effect of vincristine. Vincristine is a chemotherapy medication used to treat various types of cancer. The prevalent cutaneous side effects in patients affected by tumors undergoing chemotherapy are skin rash, xerosis, pruritus, paronychia, hair abnormality, and mucositis (Fabbrocini et al., 2012). This may suggest that our inference is again a possible true positive.

As for contraindications, Table 3 shows the top 10 predicted contraindications. Nine of ten could be found in DynaMed. The prediction atrioventricular block as a contraindication of atenolol may also be true positive as DynaMed notes that atenolol can cause atrioventricular blocks in cases of severe positioning.

Altogether, these results suggest that DDIT is a powerful computational tool that integrates multiple clinical features for the facilitation of drug repurposing.

Figure 6 represents the web interface of DDIT. Three core functions are implemented in DDIT:

– Drug/Disease: The page allows users to search for either (i) predicted Is, SEs, and CIs by inputting a drug name or (ii) predicted related that can cure, cause the disease, or as contraindicated in people with the disease, by giving the disease name.

The exact data for the RBM is a 3D array (2141*780*3) named A in .npy format in Numpy (provided in Supp. Data S1 file). The 0^th dimension represents drugs, A [i,,] means the i-th drug; the 1^th dimension represents diseases, A [,j,] means the j-th disease; the 2^th dimension represents the phenotypic types. A [i,j,k] means the k-th type between drug i and disease j. For example, A [i,j,0] denotes indication between drug i and disease j, A [i,j,1] denotes side effect between drug i and disease j, while A [i,j,2] denotes contraindication between drug i and disease j. The index i and j is calculated by the drug_id -1 and disease_id -1 respectively because the array in numpy starts from index 0 rather than 1. The mapping between drug id and drug name is provided in Supp. Data S2 file, while the mapping between disease id and disease name is provided in Supp. Data S3 file.

The ranked predictions of three types are provided in supplementary files of Supp. Data S4, Supp. Data S5, Supp. Data S6 respectively. The first column represents drug id, the second column is disease id, the third column is the prediction score, and the last column represents the status, status = 0 if the association type is not included in our dataset, status = 1 if the association type is known in our dataset.