In this study, we proposed a new computational method, HGATMDA, which combines graph embedding and graph attention network (GAT) methods to identify miRNA-disease associations. As shows in Figure 1, HGATMDA consists of four parts. Firstly, HGATMDA builds a heterogeneous graph including disease-disease sub-graph, miRNA-miRNA sub-graph and miRNA-disease sub-graph. Secondly, HGATMDA adopts a novel weighted DeepWalk to obtain dense representations of diseases and miRNAs from the disease-disease and miRNA-miRNA sub-graphs. Thirdly, it applies GAT to learn node features from the miRNA-disease sub-graph. Lastly, HGATMDA introduces a fully-connected (FC) neural network as an effective classifier for inferring the potential disease-related miRNAs.

HMDD v2.0 (Wang et al., 2010) and HMDD v3.2 (Huang et al., 2019a) databases, the experimental verified human miRNA-disease associations, are used in this paper. We directly downloaded data from http://www.cuilab.cn/hmdd, and constructed a miRNA-disease graph with known associations between miRNAs and diseases. We denote the association matrix as ANm×Nd, and N_d is the number of diseases, N_m is the number of miRNAs. The value of A_ij ∈ {0, 1} indicates whether there is a known connections existed between miRNA m_i and disease d_j. Here, N_d = 383 and N_m = 495, and 5,430 known interactions between miRNAs and diseases in the HMDD v2.0. We found that the names in HMDD v3.2, including diseases and miRNAs, did not exactly match those in HMDD V2.0. For consistency, we mapped the items in HMDD v3.2 to HMDD v2.0 as the previous work (Wang et al., 2010).

The Medical Subject Headings (MeSH) dataset contains hierarchical relationships between diseases. Tree numbers in the Mesh headings denotes parent-child associations between nodes in the networks. By using the disease name as nodes and tree numbers as edges, we can build a directed acyclic graph (DAG) with the headings for each disease (https://www.nlm.nih.gov/mesh/meshhome.html). We define DAG_d = (T_d, E_d) as disease d, and T_d and E_d are the nodes and edges in the DAG. Then, we can compute the contribution of disease d_i ∈ T_d to disease d as follows:

where Δ is a weight decay parameter and is set as 0.5 in this paper. We define the semantic value of disease d by the following formula:

Based on the assumption, the more overlap between two DAGs, the more similar the two diseases are. The similarity between disease d_i and d_j can be calculated as follows:

The basic assumption is that miRNAs with similar functions tend to be connected with similar diseases, and vice versa (Lu et al., 2008; Wang et al., 2010). Wang et al. have proposed miRNA MISIM functional similarity, that is, to calculate the similarity score by the related disease DAG between two miRNAs. Thanks to this excellent work, we can directly download the data from website (http://www.cuilab.cn/files/images/cuilab/misim.zip). We then constructed the miRNA functional network, which is denoted as FS with the shape of N_m × N_m, where N_m is the number of miRNAs. The element FS(m_i, m_j) represents the similarity score between miRNA m_i and miRNA m_j.

In order to overcome the sparsity in disease semantic similarity SS and miRNA functional similarity FS, we introduced the Gaussian interaction profile (GIP) kernel similarity. Firstly, we used column A_·i to represent a disease d_i and row A_j· to represent miRNA m_j. Then, GIP similarities between two diseases or two miRNAs are defined as follows:

here, γ_d and γ_d are parameters and can be calculated by following forums:

Based on the previous work (van Laarhoven et al., 2011), we set the bandwidth parameters γd′ and γm′ to 1 in this paper.

We observed that there was no overlap between DAGs for many diseases, resulting in zero for many elements in disease semantic similarity SS and miRNA functional similarity FS. While all entries in the GIP similarities KD and KM are non-zero. Therefore, we integrated GIP similarities with SS and FS as follows:

DeepWalk (Perozzi et al., 2014) is an algorithm that can learn the representations of vertices in graphs, inspired by the well-known unsupervised feature learning framework word2vec (Mikolov et al., 2013). Given a graph G = (V, E), V denotes all vertices in the graph, and E represents the edge or transaction matrix in the graph. In order to generate the corpus in graphs, a vertex v_i ∈ V is uniformly selected as the root, and then a random sampling is used from the neighbors as the next hop. After generating samples, DeepWalk applies SkipGram model to learn the representation of each vertex v_i, denotes as Φ(vi)∈ℝd. Usually, the dimension size d is used as a hyper-parameter.

We applied a variant DeepWalk, called weighted DeepWalk, to learn disease and miRNA representations. The element in SD and SM can be used as the weight of an edge between two vertices of the disease-disease sub-graph or the miRNA-miRNA sub-graph. Here we used the element in SD and SM to denote the probability of the current vertex walking to other vertices in the sub-graph. Given a vertex v_i, we define a walk sequence W_vi = {v_i, ⋯ } to represent the vertices passed by a weighted random walk starting at the vertex v_i. We denote the transition possibility of next hop as follows:

where v_i, v_j ∈ V are the vertices in the sub-graph, and v_j ∈ Ner_vi is in the neighbors of vertex v_i. p_ij ∈ SD or p_ij ∈ SM is the weight of an edge from v_i to v_j.

After this processing, we obtained two samples {W_d1, W_d2, ⋯ } and {W_m1, W_m2, ⋯ } as the corpus, Each sample is a sentence, and each vertex is a word. For learning the representations of diseases and miRNAs, we used DeepWalk to construct two SkipGram models for learning node representations. There are hundreds of diseases and miRNAs in the respective corpus. Therefore, hierarchical softmax wass used for accelerating the training process. We set the dimension size as F in both disease and miRNA representations, and the learned representations of diseases and miRNAs can be denoted as Xd={d1,d2,…,dNd},di∈ℝF, Xm={m1,m2,…,mNm},mi∈ℝF. F is the dimension size.

In this section, we further constructed a GAT model in the miRNA-disease sub-graph with node features learned from the previous section, to refine dense vectors of miRNAs and diseases. In consequence, these node features contain global heterogeneous graph structure information. We denote the miRNA-disease sub-graph as G_m−d = (V, E), and the number of nodes is N_d + N_m, node v_i ∈ V = {v₁, …, v_Nd, v_Nd+1, …, v_{Nd + Nm}}, edges (v_i, v_j) ∈ E = {A_ij = 1 ∈ A}. The initial node features in the miRNA-disease sub-graph are defined as X = [X_d, X_m], and xi∈ℝF is the feature of the i-th node in the graph. GAT uses multi-head attention mechanism to compute the contributions of neighbors of vertex v_i to itself. We denote the input of l-layer of our GAT model as H(l)={h1(l),h2(l),…,hN(l)},hi(l)∈ℝF(l), where N = N_d+N_m is the number of nodes, and F^(l) denotes the input dimension size of each node. Here, H⁽⁰⁾ = X is the input features of the GAT model. We define the output of l-layer as H(l+1)={h1(l+1),h2(l+1),…,hN(l+1)},hi(l+1)∈ℝF(l+1). For each node, we first compute the important score from the neighbor node j to node i

where W^(l) is a parameter matrix W^(l) ∈ ℝ^F(l+1) × F^(l), and a is a one layer feed-forward neural network. Then, we normalize the neighborhood important scores of node i by the softmax function as follows:

where Ni is the neighborhood node set of node i, including node i itself. Finally, we can use these scores to calculate the new features of node i by aggregating information from its neighbors:

where σ denotes a non-linear activation function, such as LeakyReLU. The power of GAT is benefit from the multi-head attention mechanism, we apply K independent attention of node i on its neighborhood, and the output of node features is as follows:

where Λ denotes concatenation or averaging. In our paper, all GAT layers are used concatenation, except averaging for last layer. A GAT carries out the first-order neighborhood information aggregation, and the graph convolution on multi-layers realizes multi-order neighborhood aggregation. As the number of training iterations increases, the node representations can obtain the structure information of the miRNA-disease sub-graph. Together with the node features obtained from miRNA-miRNA and disease-disease sub-graphs, the final node features of mi∈ℝF and dj∈ℝF contain rich structure information of the global heterogeneous graph.

At last, we designed a scoring function that can calculate the correlation score between a pair of miRNA and disease. We first integrated node features with raw features in SD and SM. In order to transform the raw features in SD and SM to the same dimensions as the node features, we introduced projection parameters Wd∈ℝNd×F and Wm∈ℝNm×F. Then, we can define the correlation score as follows:

where FC is a fully-connected neural networks, and the details will be discussed in the section 3.3. g(·) represents an accumulation function, such as concat(·), i.e., a concatenation of node features and raw features, or sum(·), i.e., summation of node features and raw features. Our model is trained by minimizing the cross-entropy loss and L2 regularization.

where Θ denotes parameters of our GAT model, G+ is the set of known associations between miRNA and disease, and G- is the same number set of unknown associations between miRNA and disease from negative sampling, which we will discuss the impact in the later section. X denotes the vector form of the miRNA-disease pairs in G+∪G- and y ∈ Y represents the corresponding label.

We used benchmark datasets of known miRNA-disease associations, including HMDD, dbDEMC v2.0 (Yang et al., 2017), and miR2Disease (Jiang et al., 2008). Mesh dataset is also used for computing disease semantic similarity. To evaluate the prediction performance of our method, we compared the experimental results with other stat-of-the-art methods. Moreover, we also performed case studies on three common human diseases such as breast neoplasms, lung neoplasms and kidney neoplasms for further evaluation. HMDD v2.0 is used for training, and dbDEMC v2.0, HMDD v3.2, and miR2Disease are used for validation. We carried out experimental analysis using five-fold cross validation. Specially, all known associations in HMDD v2.0 are taken as positive samples. We first selected equal number of negative samples from a uniform distribution of unknown associations between miRNAs and diseases in HMDD v2.0 dataset. Then, we randomly shuffled all the samples and divided into five equal parts, four of which are in turn used for training and the rest for validation. In addition, we applied several metrics for evaluation, such as Accuracy, Precision, Recall, F1-score. Moreover, we plotted the receiver operating characteristic curve (ROC) and precision-recall curve (PR). Furthermore, the area under the ROC curve (AUC), and the area under the P-R curve (AUPR) are used for analysis the prediction performance quantitatively.

The code of our proposed method is implemented based on the machine learning library PyTorch v1.6.0 (Paszke et al., 2019). We trained the weighted DeepWalk models using GenSim library (Rehurek and Sojka, 2010). For GAT, we used PyTorch Geometric deep learning library (Fey and Lenssen, 2019). Our experiments are run on the Ubuntu 16.04 operating system, with two Intel Xeon CPUs (2.30 GHz, 16 cores), and two Tesla V100 GPUs. The Adam optimizer is used for training (Kingma and Ba, 2017), with a learning rate of 0.001 and the weight decay is 0.00005.

In our experiments, we performed ablation study to analyze the effect of architectures and hyper-parameters, then selected the appropriate parameters. Details of the choices will be discussed in the section 3.3. We chose the model with window size of 5, walk length of 20 in weighted Deepwalk, 2 layers in GAT, and 3-layer FC as predictor. The results of our best model under five-fold CV based on HMDD v2.0 dataset are shown in Table 1 and Figure 2. HGATMDA achieves the average prediction Accuracy of 87.02%, Precision of 94.07%, Recall of 90.04%, and F1-score of 87.39%. The threshold used in theses metrics is 0.5. The mean values of AUC and AUPR are 94.54 ± 0.34 and 94.05 ± 0.18%, respectively. In summary, the results showed that HGATMDA can significantly promote the performance of predicting miRNA-disease associations.

In order to further demonstrate the prediction performance of our model, we carried out extensive experiments under different architectures and hyper-parameters, analyzed how the design of sub-model and the choice of hyper-parameters have different effect on the final performance of the proposed model. Five-fold CV on HMDD v2.0 dataset was used to evaluate the model sensitivity of architectures and hyper-parameters in the experiments. The following discussion followed the structure in Table 2 and Figures 3, 4. Dot product Predictor denotes the vector dot product of miRNA m_i and disease d_j as the predictor. FC Predictor represents a fully-connected neural networks, as shown in the Equation (15). GAT (untrained) uses randomly weights of neighbors without training, and Raw features denotes integrated similarities of miRNA and disease.

In this section, we compared the prediction performance of our model with other state-of-the-art methods, including PBMDA (You et al., 2017), GRNMF (Xiao et al., 2018), MDHGI (Chen et al., 2018e), BNPMDA (Chen et al., 2018d), MCLPMDA (Yu et al., 2019), NIMCGCN (Li et al., 2020a), and GAEMDA (Li et al., 2020b). We noted that different evaluation metrics and datasets are used in these methods. For fair comparison, AUC values under five-fold CV are selected based on HMDD v2.0 in all these studies for comparison. It is worth noticing that the AUCs reported in these papers are the best values. Therefore, we ran the five-fold CV for 100 times and picked up the best and average AUCs for comparison. The results are shown in Table 3. We can see that our model achieves the best AUC performance among these methods, with the best AUC and average AUC of 94.52 and 93.88%, respectively. In particular, NIMCGCN and GAEMDA are GCN-based methods, and AUCs of our best model are 0.96 and 1.61% higher than these two methods, respectively. This further shows that HGATMDA obtains better performance than other methods.

We performed case studies on three common human neoplasms, including Breast Neoplasms, Lung Neoplasms, Kidney Neoplasms, to further evaluate the prediction performance. In these experiments, HMDD v2.0 dataset is used to train our model, and the validation datasets are HMDD v3.2 (Huang et al., 2019a), dbDEMC v2.0 (Yang et al., 2017), and miR2Disease (Jiang et al., 2008), which are used to verify the candidates. For each specific disease, all known associations in HMDD v2.0 are taken as positive samples, while negative samples are chosen from all unknown miRNA-disease associations except the particular disease related. In the prediction phase, all candidate miRNAs are finally ranked by their correlation scores, which are the last layer of our model with sigmoid activation function, to present how much a miRNA associated with the specific disease. In addition, Top 50 candidates are listed and checked whether they are verified in the validation datasets.

We conducted the first case study for Breast Neoplasms. It is reported that Breast cancer is the leading cancer among women worldwide, and can arise for a wide number of reasons. It usually happens when cells in breast tissue grow and divide out of control. As shown in Table 4, 50 out of the top 50 candidates are confirmed in HMDD v3.2, dbDEMC v2.0, or miR2Disease. The second case study was implemented for Lung Neoplasms. This is a leading cause of cancer deaths among men and women both in United States and the world (Siegel et al., 2020). We applied our model to predict the most relevant miRNAs with Lung Neoplasms. The results are shown in the Table 5, 50 of top 50 selected miRNAs have been verified in the validation datasets. The last case study we selected was Kidney Neoplasms. It is another common disease, and the incidence still continues to increase in the United States (Siegel et al., 2020). We note that there are only 7 known associations related with Kidney Neoplasms in HMDD v2.0 dataset. We listed top 50 related miRNAs in Table 6, and we can see that 50 out of top 50 candidates are confirmed either in HMDD v3.2, dbDEMC v2.0 or miR2Disease. In particular, 6 of top 10 related miRNAs are verified in at least two datasets. Therefore, our model can serve as a powerful and effective tool to infer the potential related miRNAs for specific diseases.