We implemented LOOCV and five-fold cross validation to evaluate the performance of GRMDA. Both of LOOCV and five-fold cross validation are implemented using the recorded miRNA-disease associations in HMDD V2.0 database. During LOOCV, each one of the known miRNA-disease associations will be left out in turn to be considered as test sample. After calculating association scores of all the miRNA-disease pairs by GRMDA, we compared the score of the test sample with all the candidate pairs including all the miRNA-disease pairs which have no known associations to observe whether the rank of the test sample was above the threshold given in advance. Moreover, we plotted the true positive rate (TPR, sensitivity) vs. the false positive rate (FPR, 1-specificity) at different thresholds to obtain the Receiver operating characteristic (ROC) curves, which were shown in Figure 1. Sensitivity means the percentage of the positive samples correctly identified among all the positives; specificity refers to the percentage of negative samples correctly identified among all the negatives. Area under the ROC curve (AUC) is calculated as an index of the prediction ability of GRMDA, the value of which is between 0 and 1. Higher the AUC is, better the prediction performance will be. If AUC is smaller than 0.5, it means that the model performs not better than random prediction. As a result, GRMDA obtained the AUC of 0.8272 in the LOOCV as shown in Figure 1. The AUCs of WBSMDA and RKNNMDA in LOOCV are 0.8030, 0.7159 respectively. Therefore, according to the LOOCV results of these methods, we can intuitively observe the improvement of predicting the miRNA-disease associations with GRMDA.

During five-fold cross validation, we firstly randomly divided the known miRNA-disease associations into five parts with the same size. Then, each one of the five parts was treated as test samples in turn and the other four parts were treated as training samples. All of those miRNA-disease pairs that have no confirmed associations are candidate samples. After applying GRMDA, every score of test samples would be taken out to be compared with all the scores of candidate samples. Then we can get the rankings of test samples. In order to make the validation more accurate, we have repeated this procedure 100 times. Compared with RKNNMDA whose average AUC was 0.6723 ± 0.0027, average AUC of GRMDA in five-fold cross validation was 0.8080 ± 0.0024. The result confirmed that the GRMDA superior to RKNNMDA is able to predict miRNA-disease associations.

Based on two well-known miRNA-disease association databases, namely dbDEMC (Yang et al., 2010) and miR2Disease (Jiang et al., 2009), we studied Lymphoma to examine the practicability of GRMDA. In the end, we counted the number of the verified miRNAs in the top 10, top 20, and top 50 ones to evaluate the effectiveness of GRMDA.

Lymphoma is a group of blood cell tumors that develop from lymphocytes (a type of white blood cell) (Anagnostopoulos et al., 2000). According to the type of oncocyte, Lymphoma is divided into Hodgkin lymphoma (HL) and Non-Hodgkin Lymphoma (NHL) (Good and Gascoyne, 2008). About 90 percent of people who suffer Lymphoma have NHL (Alizadeh et al., 2000). Recent experimental studies showed the effect of re-expression of miRNA-150 on the formation of EBV-positive Burkitt lymphoma (Chen et al., 2013). A distinct set of five miRNAs (miR-150, miR-550, miR-124a, miR-518b, and miR-539) was shown to be differentially expressed in gastritis as opposed to MALT lymphoma (Thorns et al., 2012). After applying GRMDA method on Lymphoma, we got the result that 8 out of top 10, 17 out of top 20, and 44 out of top 50 potential miRNAs in the prediction result list for Lymphoma have been experimentally verified according to dbDEMC and miR2Disease (see Table 1). Compered with RKNNMDA and WBSMDA whose confirmed results are respectively 29 and 42 within top 50 predicted miRNAs for Lymphoma, GRMDA presents a more powerful predictive ability.

In order to help scientists to make use of our method and predictive results more efficiently, we have provided the prediction list of the whole potential miRNAs associated with all the human diseases and their association scores predicted by GRMDA (see Supplementary Table 1).

To estimate the applicability of GRMDA on the new diseases which do not have any known associations with miRNAs, we set all of the associations which involve the test disease as unknown ones. After implementing GRMDA, we obtained the ranking of the miRNA-disease association prediction scores. We use Breast Neoplasm as an example, the predicted result of which is shown in Table 2. From the result, we can see that 10, 20, and 50 related miRNAs out of the top 10, 20, and 50 have been confirmed by at least one of the three databases HMDD, dbDEMC and miR2Disease. The result that all the top 50 associations had been confirmed means that our method has a wonderful performance in this aspect. For example, hsa-mir-302b is ranked at top 1, which exhibits high frequency genomic alternations in human Breast Neoplasm (Zhang et al., 2006).

Finally, we used HMDD V1.0 to test GRMDA and observe whether our method has a good robustness by observing whether our method can keep a good performance in other dataset. According to the statistical results, we can see that 7, 16, and 42 respectively out of top 10, 20, and 50 miRNAs predicted to be related to the Esophageal Neoplasms have been confirmed by three databases mentioned above (see Table 3). For example, hsa-mir-196a which ranks the second in the top 50 has been confirmed that its binding-site SNP (rs6573) can regulate RAP1A expression, which contributes to the risk and metastasis of esophageal squamous cell carcinoma (Wang et al., 2012).

In conclusion, the prediction performance of GRMDA is satisfactory. Because of that, we can foresee that with the development of experimental tools and the improvement of experimental measures, more and more miRNA-disease associations predicted by our method will be confirmed in the medical laboratory.

We downloaded data about human miRNA-disease associations between 383 diseases and 495 miRNAs from HMDD V2.0 database, which includes up to 5430 associations. The miRNA-disease association matrix A was built, where the entity A(i,j) will be 1 if miRNA m(i) and disease d(j) compose an association, otherwise 0. Variables nm and nd denote the number of miRNAs and diseases respectively.

The assumption that miRNAs with functional similarity tend to be associated with diseases that have phenotypical similarity is our basis to calculate the functional similarity score between miRNAs (Wang et al., 2010). Combining known miRNA-disease associations and disease similarity, the functional similarity between two miRNAs could be obtained through measuring the similarity between two disease sets that are related to the two miRNAs. The data came from http://www.cuilab.cn/files/images/cuilab/misim.zip, according to which we constructed the matrix FS to represent miRNA functional similarity, in which the entity FS(m(i), m(j)) denotes the functional similarity score between miRNAs m(i) and m(j) with value from 0 to 1.

We use Directed Acyclic Graph (DAG) whose descriptor is DAG(D) = (D,T(D),E(D)) to represent each disease, in which T(D) is the node set composed of node D itself and its ancestor nodes, E(D) is the edge set consisting of the direct edges from parent nodes to child nodes (Wang et al., 2010). The formula to calculate the semantic value of disease D is shown as below:

Where Δ is the semantic contribution factor, whose value is between 0 and 1. For example, for a certain disease D, it contributes to the semantic value of itself with a value of 1. The farther the distance is from disease D to the disease d in T(D), the less the semantic contribution of d to D will be. Moreover, contributions from diseases in the same layer to the semantic value of disease D would be equal. The way to calculate semantic similarity between disease d (i) and d (j) comes from a reliable assumption that the larger part of the sharing of DAGs of two diseases, the larger the semantic similarity between them will be. The formula shown below is the semantic similarity between disease d(i) and d(j):

In this section, we calculated the disease semantic similarity following the method given in the reference (Xuan et al., 2013). The method in Disease semantic similarity model 1 has a good performance, however, it has a weakness. For example, if different diseases d₁and d₂are in the same layer of DAG(D), then as a result of disease semantic similarity model 1 defined in the section above, d₁and d₂ have the same contribution to the semantic value of disease D. But in certain circumstances, d₁ may appear in less disease DAGs than d₂. If that happens, it is easy to realize that d₁ is more specific than d₂ and should have a higher contribution to the semantic value of D. Therefore, we defined the contribution of disease t in DAG(D) to the semantic value of disease D as follows:

The semantic value of disease D in model 2 is calculated in the similar way as in equation (2). The way to calculate semantic similarity between disease d(i) and d(j) also has same formation with method 1. The formulas are shown below:

SS2 is the disease semantic similarity matrix calculated based on model 2 and its entity SS2(d(i),d(j)) in row i column j is the disease semantic similarity between disease d(i) and d(j) based on disease semantic similarity model 2.

It is observed that functional similar miRNAs always tend to be associated with similar diseases. Based on this observation, we could utilize the topologic information extracted from the known miRNA-disease association network to compute the Gaussian interaction profile kernel similarity for diseases. First, we defined a binary vector IP(d(i)), the same value as the ith column of our miRNA-disease association matrix A, to represent the interaction profiles of disease d(i). The formula to calculate Gaussian interaction profile kernel similarity between disease d(i) and d(j) was shown below:

The effect of parameter γ_d is to control the kernel bandwidth. γ′d is usually set as 1. γ_d is calculated by normalizing γ′d by the average number of known miRNA-disease associations for all diseases. KD is the Gaussian interaction profile kernel similarity matrix for diseases.

MiRNA Gaussian interaction profile kernel similarity matrix is calculated in a similar way with disease Gaussian interaction profile kernel similarity:

IP(m(i)) has the same value with the ith row of our miRNA-disease association matrix A to denote the interaction profiles of miRNA m(i).

In this work, the integrated disease similarity S_d were constructed based on disease semantic similarity SS and Gaussian interaction profile kernel similarity KD. Specifically, if disease d(i) and d(j) have semantic similarity, for simplicity, we assume that two types of semantic similarities between them are equally important. Then the final integrated similarity is computed directly using the average of SS1(d(i),d(j)) and SS2(d(i),d(j)), since both of two types of disease semantic similarity are calculated based on Directed Acyclic Graph (DAG) of diseases. Otherwise, if disease d(i) and d(j) do not have any semantic similarity, the integrated disease similarity equals to the Gaussian interaction profile kernel similarity which is taken as a supplement to the semantic similarity. The formula was shown as follows:

Furthermore, by supplementing the miRNA functional similarity with miRNA Gaussian interaction profile kernel similarity, we obtained the integrated miRNA similarity as follows:

We use a graph regression (Hu et al., 2015) among G_r, G_d, and G_a which represent graphs about miRNA similarity network, disease similarity network and miRNA-disease association network respectively, to predict unknown associations between miRNAs and diseases (see Figure 2). Because the graph regression is synchronously performed in miRNA similarity space, disease similarity space and miRNA-disease association space, we can get the following formula:

The first three items in formula (13) denote three low-rank decompositions to map G_r, G_d and G_a in three spaces respectively. The first item helps to decompose A into two parts, each part represents the information of G_a in miRNA or disease aspect. The second item is used to convert G_r to feature matrix about miRNA and the third item is used in the same way for generating feature matrix about disease. The fourth item means a regression between miRNA-disease association space and miRNA similarity space, from which we can get the regression matrix which can connect G_a and G_r. The fifth item represents a regression between miRNA-disease association space and disease similarity space, from which we can get the regression matrix which can connect G_a and G_d. To be specific, we map the miRNAs and diseases in G_a into an nm × r miRNA associating matrix A_r and a nd × r disease associating matrix A_d respectively. We map miRNAs in G_r into an nm × p miRNA latent feature matrix F_r. We map diseases in G_dinto a nd × q disease latent feature matrix F_d. In the end, the p × r matrix B_r and the q × r matrix B_d are the corresponding regression coefficient matrices.

It is an intricate problem to minimize the objective function as a whole. However, since our goal is to make regression between the latent association spaces and similarity spaces, it is reasonable to divide the solution into two steps: obtaining the latent matrices and regressing between the latent matrices. To make the question easier, we assume that the five items in the formula (13) are independent. Then we can easily solve the above optimization problem by minimizing the items individually. We applied SVD for low-rank decompositions to generate A_r, A_d, F_r, and F_d respectively in the following way:

SVD is an important and widely used method to decompose matrixes. SVD can condense the size of data and extract the possible association attributes between what the column and row represent. Σ is a diagonal matrix in which each value on the diagonal represents the importance of its mapped attribute. We can omit some attributes whose value is too small in Σ. However, if too small number of attributes are selected, the subtle information may be lost (Sharma, 2016). In this work, according to the particular data structure, we have retained about 45% (Franceschini et al., 2016) components for disease similarity space and miRNA similarity space. As a result, the parameters were set as q = 170 and p = 220 respectively. And for association space, the number of selected attributes was moderately set as r = 180. We operated a canonical correlation analysis on the principal components of miRNA latent feature matrix F_r and miRNA associating matrix A_r, as well as disease latent feature matrix F_d and disease associating matrix A_d respectively to check for the mutual correlation between them. The results were shown in Supplementary Figures 1, 2. After that we utilized PLS regression on the latter two items in formula (13) to generate B_r and B_d individually. For the purpose to preserve the predictive ability with small noise in PLS (Kreeger, 2013), the percentage of components to keep was set as 90% both for regression between association space and miRNA similarity space and regression between association space and disease similarity space.

According to the previous formulas, we know that F_r represents the features of miRNA, F_d represents the features of diseases, B_r represents the relation between A and F_r, B_d represents the relation between A and F_d. Then, BrBdT builds a bridge between the features of miRNAs, the features of diseases and the associations between them. In the end, the confidence scores of miRNA-disease pairs to be potential associations are calculated in the following formula:

where, C is the confidence score matrix and C(i,j) represents the associations core of miRNA m(i) and disease d(j). The higher the score is, the more likely the association exists.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.