In this study, we took advantage of miRNA–disease association data from the HMDD v2.0 database (Yang et al., 2013), which contained 5,430 verified associations between 495 miRNAs and 383 diseases. For convenient calculation, we constructed an adjacency matrix A ∈ R^nd×nm to represent the miRNA–disease relationship. We set nd and nm to represent the number of diseases and miRNAs, respectively. Specifically, the element A(i, j) is equal to 1 when disease d_i is shown to be connected with miRNA m_j; otherwise, it is set to 0. In order to clearly demonstrate the detailed information of matrix A, we visualized it in Figure 1. We used white points and black points to denote known associations and unknown associations, respectively. According to the distribution of points in Figure 1, number of known associations is far less than the number of unknown associations, which means the matrix A can be considered a sparse matrix.

Based on the theory proposed by Wang et al. (2010), disease similarity can be calculated using semantic information. We used SD₁ ∈ R^nd×nd to denote disease semantic similarity, which could be obtained by utilizing the disease arborescence attribute in Mesh database (Lowe and Barnett, 1994). In this database, disease nodes were labeled in a directed acyclic graph (DAG). Diseases that relate to the same genes are likely to have similar phenotypes; as such, disease functional similarity is calculated by the disease–gene connections (Luo et al., 2017; Jiang L. et al., 2018). In addition, we adopted disease functional similarity information from previous literature (Jiang L. et al., 2018) and utilized SD₂ ∈ R^nd×nd to denote disease functional similarity. The element SD₂(d_i, d_j) represents the value of similarity between disease d_i and disease d_j.

MicroRNAs with similar functions have a high probability of being related to diseases that are similar, and vice versa (Lu et al., 2008) (Goh, Cusick et al., Wang, Zaman et al., Lu, Zhang et al. 2008) (Goh, Cusick et al., Sanghamitra, Bandyopadhyay et al., Lu, Zhang et al. 2008) (Goh, Cusick et al., Lu, Zhang et al. 2008, Bandyopadhyay, Mitra et al. 2010)[38-40]. Therefore, we downloaded the miRNA function similarity information from http://www.cuilab.cn/files/images/cuilab/misim.zip (Wang et al., 2010). The miRNA sequence similarity information was acquired from the miRBase database (Kozomara and Griffiths-jones, 2013). For convenient and efficient calculation, we constructed the matrix SM₁ ∈ R^nm×nm and SM₂ ∈ R^nm×nm to store the miRNA functional similarity and sequence similarity data, respectively.

Gaussian interaction profile (GIP) kernel similarity was used to represent miRNA and disease similarity (Van et al., 2011; Chen et al., 2017b). First, we denoted vector K(d_i) to represent the interaction profile of disease d_i in accordance with whether d_i had a verified association with each miRNA. Similarly, we denoted vector K(m_i) to represent the interaction profile m_i in accordance with whether m_i had a verified association with each disease. The equation to calculate GIP kernel similarity for diseases is as follows:

where SD₃(d_i, d_j) indicates the GIP kernel similarity between disease d_i and disease d_j, ρ_d is applied to control kernel bandwidth. ρ_d is obtained by normalizing the original bandwidth ρd′ to the average number of verified associations with miRNAs per disease, as follows:

Similarly, the equations to calculate GIP kernel similarity for miRNAs are as follows:

where SM₃(m_i, m_j) indicates the GIP kernel similarity between miRNA m_i and miRNA m_j; ρ_m is also employed to control kernel bandwidth.

A flow chart of ILPMDA is shown in Figure 2. After we obtained disease semantic similarity, disease functional similarity, and disease GIP kernel similarity, the similarity kernel fusion (Jiang L. et al., 2018; Jiang et al., 2019) was implemented to integrate them into ultimate disease similarity. Similarly, miRNA functional similarity, miRNA sequence similarity, and miRNA GIP kernel similarity were integrated into ultimate miRNA similarity by implementing the SKF method. The specific integration process of disease similarity matrices is described in the following discussion.

First, three different disease similarities were treated as original disease similarity kernels, which were defined as SD_n, n = 1, 2, 3 in the above sections. The similarity of each disease was normalized using the following equation:

where P_n indicates the normalized kernel that satisfies the condition of ∑_dk∈D P_n(d_i, d_j) = 1, and D={di}i=1n⁢d indicates the set of diseases.

Second, the neighbor-constraint kernel for each original disease kernel was constructed by the following equation:

where C_n (d_i, d_j) denotes a neighbor-constraint kernel that satisfies the condition of ∑_dk∈D C_n(d_i, d_j) = 1 and N_i denotes the collection of all neighbors of disease d_i, including itself.

Third, the normalized kernels and neighbor-constraint kernels were integrated as follows:

where Pnl+1 represents the value of the nth kernel after l+1 iterations, Pt0 represents the initial value of P_t, and the weight parameter β ∈ (0, 1) is used to balance the rate. After Pnl+1, n = 1, 2, 3 is obtained, the overall kernel SD^∗ can be calculated by the following formula:

Fourth, a weighted matrix W is applied to further eliminate noise in the overall kernel SD^∗. The process for constructing W is as follows:

Last, the final disease similarity kernel SD ∈ R^nd×nd can be computed as follows:

Similarly, we could obtain the final miRNA similarity kernel as SM ∈ R^nm×nm.

In order to make the experimental data more accurate and improve the prediction accuracy, we applied the method of weighted k-nearest known neighbor (Ezzat et al., 2017) to the adjacency matrix A ∈ R^nm×nd. A(d_i) = (A_i1, A_i2,, A_ind) and A(m_j) = (A_1j, A_2j,, A_nmj) indicate the interaction profiles for d_i and m_j, respectively. The procedure for using the WKNKN algorithm included several steps.

First, the similarity between each disease and its k-nearest verified diseases was employed to construct the interaction profile. For example, the k interaction profiles between miRNA d_r and its KNNs are represented by AD(d_r), which is obtained by the following formula:

where N_d = ∑_1≤i≤K SD(d_i, d_j) is the normalization term. The weight coefficient ρ_i = tⁱ⁻¹*SD (d_i, d_r) is employed to control the similarity between d_i and d_r, where 0≤ t ≤ 1 is the corresponding balance parameter. Similarly, the k interaction profiles between disease m_r and its KNNs are represented by AM (m_r), which can be obtained as follows:

Second, the AD and AM are combined to obtain the matrix A_md, which indicates the new miRNA–disease interaction profile. The specific process is depicted by the following formula:

Last, the matrix A_md is employed to update the original matrix A, and the corresponding formula is displayed as follows:

When the label propagation (Zhu and Ghahramani, 2002) method was implemented for the disease similarity network SD, we applied A(i,:) to represent the initial label of disease node d_i, where A(i,:) denotes the ith row of the miRNA–disease association matrix A. In addition, this label information is propagated among neighboring nodes in similarity network SD. Thereafter, the label information of each disease node can be updated depending on the label information accepted from the neighboring nodes. However, according to the assumption that local neighboring nodes with high similarity scores are more reliable than remote neighboring nodes with low similarity scores, we employed the KNN algorithm to sort KNNs for each disease node. Hence, Q_i was utilized to denote the nearest neighboring nodes set of disease node d_i. Then, its local affinity could be calculated as follows:

where sd denotes the local affinity matrix of the disease. Similarly, we could obtain the miRNA local affinity matrixsm. Thereafter, we constructed novel disease and miRNA weighted similarity networks that are more suitable for implementing the label propagation algorithm.

After obtaining the disease similarity matrix sd ∈ R^nd×nd, miRNA similarity matrix sm ∈ R^nm×nm, and the miRNA–disease association matrix A ∈ R^nd×nm, we applied the bidirectional label propagation algorithm. The implantation of directional label propagation can be divided into three major steps.

In the first step, updating the label information of a specific disease d_i is affected by two parts of the labels. This involves absorbing labels from neighboring nodes and retaining previous labels. Dit is employed to represent the label of disease d_i after t rounds of updating; then, Dit can be calculated using the following equation:

where γ ∈ [0, 1] is employed to balance the rate between absorbing labels from neighboring nodes and retaining previous labels, and Di0=A⁢(i,:) represents the initial label information of disease d_i.

For all disease nodes, we could obtain their label vectors D1t, D2t, …, Dndt after t rounds until convergence. We constructed the formula Dt=(D1t; D2t; …; Dndt, and Equation (16) could be rewritten as follows:

The iteration of the above equation can be considered convergent. As such, when the difference between the last label matrix D^t and the former label matrix D^t−1 is lower than the predetermined threshold, the process of iterative updating will stop. Hence, we assumed that the iterative process stopped after kd rounds of iterative updating. One miRNA–disease association score matrix F_d could be obtained by the below formula:

In the second step, for a random miRNA m_i, the process of updating its label information is also affected by absorbing labels from neighbors with γ probability and remaining previous labels with 1−γ probability. Mit is used to represent the label of miRNA m_i after t rounds of updating and Mi0=A⁢(:,i) denotes the initial label information of miRNA m_i. Then, Mit can be calculated using the following formula:

In all the miRNA nodes, we could obtain their label vectors M1t, M2t, …, Mndt after t rounds until convergence. We constructed the formula Mt=(M1t; M2t; …; Mnmt, and Equation (19) could be rewritten as follows:

The iteration of the above equation can be considered convergent. Therefore, when the difference between last label matrix M^t and the former label matrix M^t−1 is lower than the predetermined threshold, the process of iterative updating will stop. After km rounds of iterations, we could obtain another miRNA–disease association score matrix:

In the third step, for the purpose of generating accuracy prediction results, the score matrices F_m and F_d can be integrated by an average ensemble method to obtain a final miRNA–disease association score matrix, as follows:

where hyper-parameter α is set to 0.5, which is employed to balance the score matrices F_d and F_m. The adjacency matrix F contains the predictive association score between miRNAs and diseases.

We utilized global LOOCV and 5-CV to evaluate the prediction performance of ILPMDA. As cross-validation is currently general practice in predicting miRNA/circular RNA (circRNA)/long noncoding RNA (lncRNA)–disease associations, this approach was selected for evaluation (Shen et al., 2017; Li et al., 2019; Wang et al., 2020). In the global LOOCV method, each verified association in the HMDD v2.0 database served as a test sample; the rest of the verified associations served as training samples, and the unknown associations were regarded as candidate samples. Similarly, in the 5-CV method, all verified associations in the HMDD v2.0 database were randomly divided into five parts in a random way; one group was treated as test samples, while the others were treated as training samples. The candidate is composed of unknown miRNA–disease associations. We applied repeated segmentations on verified positive samples multiple times to reduce potential deviations. We then implemented the model to calculate the score list of all associations and rank the candidate sample scores with the sample score in both LOOCV and 5-CV. If the ranking of the test sample was higher than the given threshold, the model was regarded as a successful prediction model. In addition, we could draw the ROC curve by employing the true-positive rate (TPR, sensitivity) against the false-positive rate (FPR, 1 − specificity). Sensitivity refers to the percentage of test samples with values higher than the threshold, and specificity refers to the percentage of negative associations with values lower than the threshold. The equations used to calculate FPR and TPR are demonstrated as follows:

where FP indicates the quantity of samples that are negative samples but are considered as positive samples, TP denotes the quantity of samples that are positive samples and are considered as positive samples, and FN and TN are identified as the opposite of FP and TP, respectively. Hence, we could utilize the AUC values between 0 and 1 as evaluative criteria; a greater AUC indicated that the model had better prediction performance. Simulation results showed that ILPMDA achieved AUCs of 0.9751 and 0.9501 for the global LOOCV and 5-CV methods, respectively, demonstrating that ILPMDA achieved excellent prediction performance.

We also optimized two significant parameters γ and α in the ILPMDA model, which were utilized to balance the rate of absorbing labels from neighboring nodes and retaining previous labels and the rate between the score matrix F_d and score matrix F_m. In order to select the best value of parameter γ, we applied the AUCs of 5-CV to evaluate γ ∈ {0, 0.1, 0.2,,1}. According to the evaluation results demonstrated in Figure 3A, ILPMDA could obtain the highest AUC score when γ = 0.2. Thus, more previous label information for the node should be preserved when updating its label information. In addition, we also applied the AUCs of 5-CV to evaluate α ∈ {0, 0.1, 0.2,,1} for selecting the best value of parameter α. According to the evaluation results shown in Figure 3B, ILPMDA could obtain the highest AUC score while α = 0.3. This finding shows that the weight assigned to score matrix F_d should be greater than that of score matrix F_m. The reason for this may be related to the number of collected miRNAs being higher than the number of collected diseases. In conclusion, the parameters γ and α in ILPMDA were set as 0.2 and 0.3, respectively.

We first compared ILPMDA with several recent computational models [MSCHLMDA (Wu et al., 2020), GRL_2,1-NMF (Gao et al., 2020), ICFMDA (Jiang Y. et al., 2018), and SACMDA (Shao et al., 2018)] to demonstrate its superior performance via the global LOOCV and 5-CV methods. Here, multi-similarity-based combinative hypergraph learning for predicting miRNA–disease association (MSCHLMDA) applied the KNN and k-means algorithms to establish different hypergraphs, which were combined to predict potential miRNA–disease associations. GRL_2,1−NMFutilized the Laplacian regularized L_2,1-nonnegative matrix factorization method to observe unknown miRNA–disease associations; ICFMDA incorporated similarity matrices to improve the collaborative filtering method for predicting more newer miRNA–disease pairs; and SACMDA utilized disease and miRNA information to construct a heterogeneous graph and used short acyclic connections in the heterogeneous graph to infer miRNA–disease associations. As illustrated in Figures 4A,B, ILPMDA achieved AUCs of 0.9751 and 0.9501 via global LOOCV and 5-CV, respectively. Both ranked the highest when compared with MSCHLMDA, GRL_2,1-NMDA, ICFMDA, and SACMDA, which achieved AUCs of 0.9287, 0.9280, 0.9072, and 0.8777 in global LOOCV and 0.9263, 0.9276, 0.9046, and 0.8773 in 5-CV, respectively.

In addition, we compared ILPMDA with several LP-based models [MCLPMDA (Yu et al., 2019), HLPMDA (Chen et al., 2018b), and MDLPMDA (Qu et al., 2019)] to evaluate its prediction ability in the frameworks of global LOOCV and 5-CV. Here, MCLPMDA applied the matrix completion method to deal with similarities. Then, label propagation was applied to novel miRNA and disease similarity and association matrices to observe disease-related miRNAs. HLPMDA applied the heterogeneous label propagation algorithm on a multi-network of miRNAs, lncRNAs, and diseases to predict unobserved miRNA–disease interactions. MDLPMDA utilized the matrix decomposition method to obtain a novel association matrix with less noise and then applied the label propagation method to common miRNA and disease similarity matrices and the novel association matrix to infer unverified miRNA–disease associations. As illustrated in Table 1, MCLPMDA, HLPMDA, MDLPMDA, and ILPMDA achieved AUCs of 0.9320, 0.9218, 0.9211, and 0.9501, respectively, via the 5-CV method. By means of the global LOOCV method, MCLPMDA, HLPMDA, MDLPMDA, and ILPMDA achieved AUCs of 0.9410, 0.9232, 0.9222, and 0.9751, respectively. According to the above analysis, the prediction performance of ILPMDA is greater than that of previous computational models.

Moreover, we compared SKF with similarity network fusion (SNF) (Chen et al., 2019b) and average fusion to verify the superiority of SKF for integrating biology data. SNF fuses multiple pieces of complementary data to obtain integrated information, and the AF integrates different data by averaging multiple similarity matrices. In order to make a fair comparison, we used SNF and AF to replace SKF to process similarity data, while the remainder of the model remained unchanged. We also used the AUC values of 5-CV to evaluate the data integration capability of SKF, SNF, and AF. As shown in Figure 5, SKF, SNF, and AF obtained AUCs of 0.9501, 0.9364, and 0.9302, respectively. These findings demonstrate that SKF performs better than SNF and AF in terms of integrating different similarity information.

Furthermore, we also investigated the effect of the WKNKN algorithm for known miRNA–disease associations on model performance. We implemented two methods of global LOOCV and 5-CV, and then plotted the ROC curves, as shown in Figure 6. In ILPMDA, WKNKN considers the sparsity of the original association matrix, thereby improving the prediction performance of the model. By contrast, ILPMDA without WKNKN disregards the sparsity of the association data; thus, the predictive performance is also reduced. Based on the results, the AUC values of ILPMDA based on global LOOCV and 5-CV were 0.9751 and 0.9501, respectively. On the other hand, the AUC values of ILPMDA without WKNKN based on global LOOCV and 5-CV were 0.9494 and 0.9354, respectively. It is apparent that ILPMDA with WKNKN has higher AUC values compared to that without WKNKN.

To demonstrate the prediction ability of ILPMDA, we implemented two common human diseases (colon neoplasms and prostate neoplasms) to perform a kind of case study. For a random disease, known associations of whole diseases in the HMDD v2.0 database were considered as training samples, while unknown associations were treated as candidate samples. We ranked the predicted association score of the candidate samples after performing ILPMDA; then, the top 50 candidate associations with the specific disease were selected and confirmed by the miR2Disease and dbDEMC v2.0 databases (Jiang et al., 2009; Yang et al., 2016). After we compared the information of the HMDD v2.0 database with that of the miR2Ddisease and dbDEMC databases, we found that 232 of the 5,430 verified associations in the HMDD v2.0 database also appeared in the miR2Disease database; meanwhile, 546 of the 5,430 verified associations in the HMDD v2.0 database also appeared in the dbDEMC v2.0 database. Furthermore, because only candidate samples for the specific disease were ranked and verified, the prediction list had no overlapping miRNAs in the training samples.

Colon neoplasms is acknowledged as the third gastrointestinal disease in the medical field (Torre et al., 2015; Bao et al., 2017). In addition, several potential miRNA–colon neoplasm connections have been observed in previous experiments (Hiroko et al., 2014; Rotelli et al., 2015), including miR-17, miR-21, and miR-31. These studies have shown that miRNAs can be utilized as key biomarkers for colon neoplasms. Hence, observing miRNA–colon neoplasm interactions can contribute to the diagnosis and treatment of colon neoplasms. After we ranked the prediction results of our model based on prediction score, 48 of the top 50 miRNAs were confirmed to be related to colon neoplasms according to the miR2Disease and dbDEMC v2.0 databases (Table 2).

Prostate neoplasms is regarded as the disease with the highest incidence rate among men; furthermore, these have been observed to associate with a part of some miRNAs in clinical experiments (Goto et al., 2015). For example, the expression of miR-183 in prostate cells and tissues is significantly higher than that in corresponding normal prostate cells and tissues. In conclusion, prostate neoplasms can be treated by inhibiting miR-183 expression (Ueno et al., 2013). After we ranked the prediction results of ILPMDA according to the prediction score, 46 of the top 50 miRNAs were confirmed to be associated with prostate neoplasms through the miR2Disease and dbDEMC v2.0 databases (Table 3).

Next, we carried out another case study on breast neoplasms to illustrate the applicability of ILPMDA to new diseases. Breast neoplasms are often seen as a common disease in females, which have great negative effects on women’s health. Several miRNAs associated with breast neoplasms have been found by biological experiments in previous years (Fu et al., 2011). For example, the downregulation of miRNA-140 promoted cancer stem cell formation in basal-like early stage breast cancer (Li et al., 2014). In this case study, we deleted all miRNA–breast neoplasms association information from the HMDD v2.0 database, so breast neoplasms would be considered as a new disease without related miRNAs. We also used dbDEMC v2.0 and miR2Disease databases to validate predicted miRNAs related to breast neoplasms that were acquired by operating ILPMDA model. As shown in Table 4, 48 out of top 50 ranked miRNAs were verified by dbDEMC v2.0 database and miR2Disease database. Consequently, ILPMDA could be implemented to observe unverified miRNA–new disease associations.

According to the above analysis, the case studies of colon neoplasms and prostate neoplasms further demonstrate the utility of our model in predicting unknown miRNA–disease associations.