piRNA-disease association prediction based on multi-channel graph variational autoencoder

Benchmark datasets

The benchmark dataset is constructed from piRDisease v1.0 (Muhammad et al., 2019), which is the same version used by iPiDi-PUL (Wei, Xu & Liu, 2021), iPiDA-sHN (Wei, Ding & Liu, 2020), and DFL-PiDA (Ji et al., 2021), which offers the information of piRNA-disease associations. There are 7,939 associations in this database, labeled by manual work. In this study, we used the same data preprocessing steps as iPiDA-sHN (Wei, Ding & Liu, 2020). These steps include filtering redundant and non-human piRNAs, removing piRNAs that have no sequence information in the database. After applying these steps, we obtained a benchmark dataset that is identical to the datasets used by iPiDi-PUL (Wei, Xu & Liu, 2021), iPiDA-sHN (Wei, Ding & Liu, 2020), and DFL-PiDA (Ji et al., 2021), containing 5,002 experimentally validated associations with 4,350 piRNAs and 21 diseases. The benchmark dataset is a publicly accessible resource at http://bliulab.net/iPiDi-PUL/dataset/.

Computation method

This study proposes a MC-GVAE model to predict the potential associations between piRNAs and diseases. The model is a pipeline with four sequential steps: data preprocessing, similarity network construction, representation learning based on Multi-channel graph VAE, and prediction. In the step of data preprocessing, the piRNA sequence expression and the adjacency matrix of piRNA-disease associations are extracted from the piRDisease V1.0 database. The adjacency matrix is as follows:

(1) $\mathrm{A}=\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \dots & a_{1,n}\\ a_{2,1}& a_{2,2}& \dots & a_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m,1}& a_{m,2}& \dots & a_{m,n}\end{array}\right]$where ${\mathrm{a}}_{\mathrm{i},\mathrm{j}}=1$ if the i-th RNA is associated with the j-th disease, otherwise ${\mathrm{a}}_{\mathrm{i},\mathrm{j}}=0$. Then, the k-mer algorithm is used to calculate the sequence features of each piRNA, and the Pearson product-moment correlation coefficient is used to obtain the sequence similarity matrix of piRNAs. Meanwhile, based on the adjacency matrix of piRNA-disease associations, the Gaussian kernel function is used to calculate the GIP similarity matrix of piRNAs and diseases. In addition, the directed acyclic graph (DAG) of each disease in the associated matrix is obtained from the disease descriptor vocabulary MeSH, and the semantic value of each disease is calculated, resulting in the semantic similarity matrix of diseases. In the step of similarity network construction, four different similarity subnetworks are constructed based on the four different similarity features obtained in the previous step. In the step of representation learning, four identical VGAE models are used to separately learn the feature representations of piRNAs and diseases from four different similarity subnetworks, and then concatenate the four different feature representations to better mine the associations between piRNAs and diseases. In the step of prediction, a three-layer fully connected neural network is used to classify and predict all sample pairs, obtaining the final prediction results.The flowchart of MC-GVAE’s computation process is shown in Fig. 1. The following sections will give a detailed introduction to each step’s specific implementation.

Similarity network construction

To construct the similarity network, we computed the similarity features between piRNAs and diseases from four different perspectives, namely piRNA sequence similarity, disease semantic similarity, piRNA GIP kernel similarity, and disease GIP kernel similarity. These similarity features can reflect the different correlation patterns between piRNAs and diseases, which can help our model to better mine the associations between piRNAs and diseases.

Representation learning based on multi-channel graph VAE

A graph VAE is a kind of variational autoencoder applied to graphs composed of GCN and VAE and is an unsupervised learning framework. Therefore, the graph VAE not only learns graph structure features but also generates data distribution features of piRNA and disease. In the section that follows, GCN and graph VAE will be briefly introduced.

Predicting piRNA-disease association based on neural network

Firstly, the representation of piRNA-disease association pairs is obtained by concatenating four features $Z_{\mathrm{P}\mathrm{S}}$, $Z_{\mathrm{P}\mathrm{G}}$, $Z_{\mathrm{D}\mathrm{S}}$, and $Z_{\mathrm{d}\mathrm{G}}$, where $Z_{\mathrm{P}\mathrm{S}}$ denotes the distribution representation learned from the piRNA sequence similarity network, $Z_{\mathrm{D}\mathrm{S}}$ from the disease semantic similarity network, $Z_{\mathrm{P}\mathrm{G}}$ from the piRNA GIP kernel similarity network, and $Z_{\mathrm{d}\mathrm{G}}$ from the disease GIP kernel similarity network:

(15) $$F=\left(Z_{PS},Z_{DS},Z_{PG},Z_{dG}\right)$$

Then, the final representation $F$ was fed into a three layers full connection neural network classifier, and get the prediction score ${\hat{y}}_i$ as follows:

(16) $${\hat{y}}_i=sigmoid\left(MLP\left(F\right)\right)$$where $sigmoid()$ is an activation function to transform ${\hat{\mathit{y}}}_{\mathit{i}}$ to be a number between 0 and 1. In this study, the piRNA-disease association pairs with ${\hat{\mathit{y}}}_{\mathit{i}}>0.5$ as 1 (positive sample). The otherwise labeled as 0 (negative sample). However, unlike graph VAE, only the binary cross-entropy was utilized as loss function for prediction.

Experimental setup and hyperparameter configuration

In the experiments, a five-fold cross-validation is conducted in the training stage in order to make this model more reliable. More specifically, samples in the training data are divided equally into five subsets, four of which are utilized for training and one of which for testing. Note that each subset has one chance to be the test data. The evaluation metrics that are mostly employed in association prediction tasks, including accuracy, recall, precision, F-1 score, and Matthews correlation coefficient (MCC).

To visually illustrate the performance, this study employs the receiver operating characteristic (ROC) curve and the precision-recall (PR) curve. The ROC curve reflects the relationship between the true positive rate and the false positive rate of the model at different thresholds, while the PR curve reflects the relationship between the precision and recall of the model at different thresholds. In addition, the Area Under Curve (AUC) and the area under PR curve (AUPR) were also employed as evaluation metrics, which are frequently used to measure the performance of classification or prediction tasks. AUC and AUPR are obtained by calculating the area under the ROC curve and the PR curve, respectively, and they can comprehensively reflect the performance of the model at different thresholds.

The experiment is based on the Python deep learning framework PyTorh v1.11.0, and the GCN model is performed by adopting Deep Graph Library (DGL) v0.8.1, a highly-performant package for graph neural networks (Wang et al., 2019b). The computational method of this article was performed on the Ubuntu 20.04 with two 1080Ti GPUs. Moreover, several important hyperparameters settings of this study were presented. According to the previous work (Kipf & Welling, 2016b), GCN was also adopted as an encoder, and the experiment demonstrates that if the number of hidden layers is over 2, the prediction performance of the model will be declined. The dimension of the vector in hidden layer1 and hidden layer2 were set as 48 and 16, and the learning rate is determined as 0.001. The training epoch is detected as 50, and the Adam optimizer and loss function with cross entropy and KL divergence were employed to update the parameters of MC-GVAE in each epoch. The following Table 1. provides a detailed overview of the model architecture and hyperparameter settings for GVAE and the classifier.

Predicting piRNA-disease association based on neural network

A benchmark dataset (Wei, Xu & Liu, 2021) based on the piRDisease database v1.0 was made use of. The benchmark dataset included 10,004 piRNA-disease association pairs, in which the positive sample pairs contained 5,002 associations with 4,350 piRNAs and 21 diseases. In addition, 5,002 negative sample pairs were randomly selected from unknown piRNA-disease association sample pairs. Finally, the detail of the results were shown in Table 2, and the results of AUC and AUPR are shown in Figs. 2 and 3, respectively. Table 2 shows the average performance of model MC-GVAE on the benchmark dataset based on five-fold cross-validation.

Table 2 shows that MC-GVAE achieved high levels of performance on all metrics, with average values of 0.9310, 0.9247, 0.9020, 0.9390, 0.8601, 0.8977, and 0.8070, respectively. This indicates that MC-GVAE can effectively learn the complex features of piRNA and disease, and use the graph structure information to improve the prediction accuracy. From the description of Fig. 2, the ROC curve is close to the upper left corner, the P-R curve is close to the upper right corner, and the curves of each fold are very close, which shows that MC-GVAE has excellent accuracy and robustness in distinguishing positive and negative samples and precisely retrieving positive samples.

The effect of GCN encoding module

In order to explore the impact of different encoders, we make a comparative experiment, where two encoders are built: one is based on a GCN model, and the other is implemented with a linear layer neural network. The decoder and the training objective are kept unchanged. AUC and AUPR scores of the two models VAE and MC-GVAE on the benchmark dataset are shown in Figs. 4 and 5, respectively.

It can be seen from the Figs. 4 and 5 that the AUC and AUPR scores of MC-GVAE are higher than those of VAE, respectively, under the same benchmark dataset based on 5-fold cross validation, which also means that the prediction performance of MC-GVAE is significantly better than that of MC-GVAE. It indicates that the GCN encoder can effectively obtain the structure information of four views graph, and have an important impact on MC-GVAE model.

Performance comparison of graph representation learning

To investigate the performance of different graph neural networks in modeling the similarity knowledge with four similarity networks, three models, GCN, GAT and Graph Transformer model—Graphormer (Ying et al., 2021), were employed. A two-layer GCN model was built to obtain low-dimensional feature representations of each view, and the numbers of the hidden units in the first and the second layer were 48 and 16, respectively. A single-layer GAT with an eight-head self-attention mechanism was also applied to learn the feature representations of four different views. Furthermore, the Graphormer model was also applied to learn feature representations, utilizing the GraphormerLayer provided in the DGL library. The hidden layer was designed with 48 units, num_heads set to 3, and a Dropout rate of 0.1 for the feature layer.

The number of parameters in all models is the same in the hidden layer, and the dimension of the input and the output vector were set to 48 and 16, respectively. The learning rate of the model is 0.001 during the training process.

Using a five-fold cross-validation method, four GNN models were trained on the benchmark dataset and their prediction performance was compared using common evaluation indicators, including AUC, AUPR, precision, recall, F1-scores, and MCC, for the purpose of comparing the performance of representation learning. Finally, the average of these indicators was listed as experimental results in Table 3.

Experimental results indicate that the graph variational autoencoder (GVAE) model surpasses other models across all evaluation metrics, particularly in terms of area under the curve (AUC), area under the precision-recall curve (AUPR), precision, recall, F1 score, and MCC. These findings affirm the efficacy and superiority of the graph VAE in the task of predicting piRNA-disease associations. Although the Graphormer, as a Graph Transformer model, demonstrates potential abilities in processing graph-structured data, it remains inferior to the GVAE models, which is built for capturing complex relations over graphs. Moreover, the generative learning approach of the graph VAE offers enhanced interpretability in understanding the associations between piRNAs and diseases. Consequently, this research underscores the significance of graph VAE in such tasks and provides new directions for future studies.

Ablation study

The feature representation of two different views was concatenated as the final classification feature of piRNA-disease association pairs in this model. In order to explore the contribution of the features of two views to the overall performance of MC-GVAE, an ablation study was performed under the benchmark dataset on five-fold cross validation. Some variants of MC-GVAE model were as follows:

• MC-GVAE without graph VAE: the graph VAE is removed, and four similarity features are directly fed into a three-layer full-connected neural network to predict the piRNA-disease association.

• MC-GVAE without NN: it only uses the inner product to reconstruct two prediction matrices based on the features of two views and gets the final prediction matrix by averaging two re-construction matrices.

• MC-GVAE without association: only two similarity calculation data, piRNA sequence and disease semantics, are employed to build the input data to learn node features of piRNA and disease.

• MC-GVAE without attribution: only the GIP similarity data for piRNA and disease are employed as the input data to learn node features of piRNA and disease.

Figure 6 shows the results of the models mentioned above. It can be seen that the original model still achieves the best performance, while MC-GVAE without graph VAE has the lowest AUC and AUPR performance. It can also be seen that, the performance of MC-GVAE without association view and MC-GVAE without attribution view based on different views by using the same computing framework is similar, and both achieve remarkable results. In addition, the best performance based on two different view features was achieved by this model MC-GVAE.

To quantify the impact of each component, statistical significance tests were conducted for the four variants and MC-GVAE. As demonstrated by the T-test results in Table 4, the MC-GVAE without graph VAE component and the MC-GVAE without NN component resulted in a significant decrease in performance, with statistical values of 61.905 and 68.509, respectively, and p-values of 0.000, indicating the importance of these two components to the model’s performance. In contrast, the variants without association and without attribution did not show statistically significant performance decreases, with statistical values of 0.849 and −0.382, respectively, and p-values of 0.396 and 0.703.

These results underscore the significance of the graph VAE and n three-layer full-connected neural network components in our model. The original MC-GVAE model achieved optimal performance in integrating feature representations from two different views, demonstrating that our proposed multi-view computational framework can significantly enhance the performance of piRNA-disease association prediction tasks.

Comparison with other methods

In this part, the performance of MC-GVAE and comparative methods including iPiDi-PUL (Wei, Xu & Liu, 2021), iPiDA-sHN (Wei, Ding & Liu, 2020), and DFL-PiDA (Ji et al., 2021), based on the same benchmark dataset under five-fold CV.

• iPiDi-PUL is a widely used and cited method for piRNA-disease association prediction, which also adopts four types of features like our method. Specifically, it adopts an ensemble learning method, which employs multiple random forests with different depths to predict piRNA-disease association. It extracts key features by principal component analysis (PCA) (Abdi & Williams, 2010), and gets various prediction scores via multiple random forests classifiers. Finally, averaging all classifier scores is used by the model. Different from our study, it does not use a graph-based learning model to capture the complex relationships between piRNAs and diseases, but relies on PCA to extract features. Choosing iPiDi-PUL as a comparative method may help to study the performance of graph-based learning models against models of reduction dimensions in representation learning.

• iPiDA-sHN is a state-of-the-art method for piRNA-disease association prediction. It is a two-step PUL model. It trains the initial classifier with positive samples and randomly selected negative samples, and then trains another classifier using positive samples and negative samples acquired in the previous step. Support Vector Machine (SVM) is leveraged to build the classifier. Different from our method, it does not use a graph variational autoencoder, but a CNN, to learn the latent representations of piRNAs and diseases. Choosing iPiDA-sHN as a comparative method may help to study the performance of graph variational autoencoder against other graph-based learning models.

• DFL-PiDA is constructed from a convolutional denoising autoencoder neural network combined with an extreme learning machine (ELM) to predict disease-related piRNAs. It utilizes a single hidden layer feedforward neural network to build the ELM. Different from our method, it use a convolutional denoising autoencoder to learn feature representations. Choosing DFL-PiDA as a comparative method may help to study the performance of graph variational autoencoder against other autoencoder models.

Our MC-GVAE method was performed on the Ubuntu 20.04 with two 1080Ti GPUs environment. The results are shown in Table 5. We can see that, the MC-GVAE model outperforms other methods in terms of both average AUC and AUPR metrics, demonstrating that it can effectively predicts piRNA-disease associations and has higher prediction accuracy and stability. Meanwhile, there are some possible reasons of low performance of the three models: (1) iPiDi-PUL employs PCA to extract key features, which may overlook some relevant information; (2) iPiDA-sHN adopts SVM to construct the classifier, which may be influenced by the sample imbalance; (3) DFL-PiDA merely fuses four similarity matrices, which fails to learn the rich feature representations of piRNA and disease. In addition, iPiDA-sHN and DFL-PiDA outperform iPiDi-PUL, which show that neural networks are effective than traditional models.

Case study

In this part, the prediction performance of MC-GVAE was further evaluated, three diseases, renal cell carcinoma, Alzheimer’s disease, and cardiovascular diseases with the largest number of associations in the association matrix were selected as case studies. Specifically, the association information of three diseases-related was firstly removed from the benchmark dataset, and the data remained was employed as the positive samples for training. Then, three diseases-related associations were used as positive samples in each prediction process, and the negative samples from unknown associations of three diseases with other piRNAs. Finally, the prediction scores between three diseases and other known piRNAs were obtained, respectively, and the top 20 piRNAs-related with the highest prediction score were selected for analysis.

Tables 6–8 shows the top 20 renal cell carcinoma-related, the top 20 predicted Alzheimer’s diseases-related, and the top 20 predicted cardio-vascular diseases-related piRNAs which achieves the high scores of disease-associated prediction, respectively. As shown in Table 6, only one association is still unconfirmed in renal cell carcinoma. In Table 7, eight piRNA candidates of the top 20 in Alzheimer’s disease are verified by piRDisease v1.0 database, and the other five candidates’ piRNAs are confirmed by published literature. However, seven piRNA candidates with the possible association are unconfirmed. In Table 8, all piRNA candidates in top 20 are predicted for cardio-vascular diseases. In other words, those piRNA-disease associations predicted only according to the known knowledge by our model are verified by biological experiments. It indicates that iPiDA-VAGE is effective in piRNA-disease association prediction.