In this section, we will describe our proposed method BCGP that we use to generate pre-trained embeddings for lncRNAs, circRNAs, or miRNAs and fine-tune with general Graph Neural Networks to recover the unseen associations between them. Additionally, we integrate a novel training method, Neural Common Neighbour (NCN) (Wang X. et al., 2023) to further enhance the performance of the fine-tuning GNNs. In Section 2.2, we will present the notations used in this article and briefly describe our research problem formulation. Second, in Section 2.3, we demonstrate the pre-training stage of BCGP, where we use the k -mers method to split all RNA sequences into fragments of equal length, then we use Bidirectional Encoder Representations from Transformers (BERT) (Devlin et al., 2018), which is a transformer-based contextualized language representation model to generate pre-trained embeddings. Third, in Section 2.4, we describe a general fine-tuning method to incorporate a GNN to leverage the pre-trained embeddings obtained from Section 2.3. To boost the performance of the fine-tuning GNN, we also present the integration of a novel training method named Neural Common Neighbour. We use Figure 1 to illustrate the overview of our BCGP method. The detailed mathematical formulation of our model is provided in Equations 1–13.

To initialize our task, we first define all entities and their associated information. In particular, we denote sets of lncRNAs, circRNAs and miRNAs as V lnc , V circ and V mi respectively. Considering the associations between lncRNAs, circRNAs, and miRNAs we construct an undirected graph G = { V , E , A } , where V is the set of all RNAs ( V = V lnc ∪ V circ ∪ V mi ) , E ⊆ V × V represents the set of associations between RNAs and adjacency matrix A ∈ R N × N is a symmetric matrix, which is defined as follows: A u v = 1 if  u , v ∈ E , 0 otherwise , (1) where N = | V | is the number of all types of RNAs. The degree of node u is d ( u , A ) ≔ ∑ v = 1 N A u v . The neighbours of node u are the nodes connected to u , which is defined as N e ( u , A ) ≔ { v ∣ v ∈ V , A u v > 0 } . For brevity, we use N e ( u ) to denote N e ( u , A ) since A is fixed. The common neighbour refers to the nodes connected to i and j : N e ( i ) ∩ N e ( j ) . Let S be the collection of RNA sequences, including all lncRNAs, circRNAs and miRNAs. Specifically, S = { s 1 , s 2 , … , s N } , where each s i represents a unique RNA sequence.

In this next section, we will provide a detailed description of our proposed method, including the techniques in the pre-training and fine-tuning stages.

During the pre-training stage, our goal is to effectively project the information of RNAs from their nucleotide sequences to latent embeddings. Specifically, we aim to simplify computational demands while enriching the semantic content of the embeddings, which captures the intrinsic patterns and relations within RNA sequences. Such embeddings not only preserve the biological significance and genetic information of RNA sequences but also serve as the initial features of nodes for the subsequent fine-tuning stage so that those sequence information can be encoded into the graph representation.

Instead of treating each base (A, C, G, U) as an individual token, given an RNA sequence s i , we leverage the k -mers tokenization to segment s i into overlapping and equal-length segment, with k indicating the length of each segment. Let C = { C 1 , C 2 , … , C N } be the collection of RNA k -mer sequences, and C i can be denoted as follows: C i = s i , j : j + k ∣ j = 1,2 , … , M , (2) where M = L i − k + 1 represents the total number of k -mers of s i and L i represents the length of s i . For example, the RNA sequence ‘UAACAC’ can be tokenized to a sequence of four 3-mers: UAA, AAC, ACA, CAC. This method not only enables the implementation of sequence embedding algorithms, but also deepens the understanding of richer contextual information for each nucleic acid sequence.

After tokenizing each RNA sequence into the overlapping segment using k -mer tokenization, we utilize BERT to capture both contextual and structural information from the whole corpus of RNA sequences with an attention mechanism. As a widely used transformer-based language representation model, BERT enables the generation of contextualized representations that consider the global context of the entire sequence, allowing the identification of the intricate patterns and relationships within the RNA sequences.

Given a sequence of k -mers C i derived from an RNA sequence s i , we first initialize each k -mer into a high-dimensional vector through an embedding process, resulting in the embedding matrix X i ∈ R M × D , where D denotes the dimension of the embedding vectors. To obtain contextual and informative embedding Z i , BCGP performs the multi-head self-attention mechanism on X i , which is defined as: MultiHead X i = ⊕ head i 1 , head i 2 , … , head i H W O , (3) head i h = softmax Q i h K i h T d k V i h , (4) where Q i h = X i W h Q , K i h = X i W h K , V i h = X i W h V , (5) where ⊕ denotes the concatenation operation, W O , W h Q , W h K and W h V are the learnable parameters for linear projection for the h -th head. Here, H represents the number of attention heads, and d k is the scaling factor used to maintain numerical stability and facilitate stable gradients during training, often set to D / H . After applying the multi-head self-attention mechanism multiple times, we can derive the contextual embedding Z i from the output of the last layer, which is denoted as: Z i = MultiHead Last X i . (6)

Building upon the aforementioned self-attention mechanism, we adapt the Masked Language Modeling (MLM) to train the BERT model on RNA sequences. For each RNA sequence, we randomly select regions constituting 15 % of the sequence and mask contiguous k -mers within these regions. Using the surrounding context, the model is then trained to predict the masked k -mers. The training objective to minimize the cross-entropy loss which is defined as follows: L MLM = − ∑ i = 1 T y i ′ ⁡ log y i , (7) where y i ′ represents the one-hot encoded ground-truth vector for the masked k -mers and y i denotes the predicted probability distribution over the k -mer vocabulary for each of the T masked positions. Specifically, the predicted probability distribution y i for a masked position is defined as follows: y i = Softmax z i W c + b c , (8) where z i is the i -th encoded representation from Z i , W c and b c are the parameters of the linear classifier respectively.

After obtaining the initial RNA embeddings for lncRNAs, circRNAs and miRNAs from BERT during the pre-training stage, we construct a heterogeneous graph. This graph integrates known associations between lncRNAs and miRNAs, as well as the high-quality associations between circRNAs and miRNAs, with these relationships represented as edges. The pre-trained RNA sequence embeddings are utilized as node features in this graph. To obtain effective node embeddings, we leverage GNNs to capture the intricate and complex relationships between entities.

The current de facto design of GNNs follows the message passing framework (Kipf and Welling, 2016), which is based on the core idea of recursive neighborhood aggregation. Specifically, for an L -layer GNN, the representation learning function of the l -th layer is represented as: b i l = A G G R E G A T E h i l − 1 ∪ h j l − 1 : j ∈ N e i , (9) h i l = C O M B I N E h i l − 1 , b i l , (10) where b i ( l ) is a message vector computed from the representations of the neighbors N e ( i ) from the previous layer i.e., l − 1 -th layer, N e ( i ) is a set of nodes adjacent to v i , h i ( l ) is the representation of node v i at the l -th layer with h i ( 0 ) = z i , and A G G R E G A T E ( ⋅ ) and C O M B I N E ( ⋅ ) are the component functions of GNN layers. It is worth noting that our proposed BCGP approach is a general framework that can be incorporated with various GNNs.

To further boost the performance of association prediction, we leverage the Neural Common Neighbour (Wang X. et al., 2023) to capture more refined node features such as multi-hop structure and attribute information. Specifically, after obtaining the representation h i of node v i , instead of directly using the node representation for link prediction, BCGP focuses on the pairwise relationships between nodes, specifically leveraging the common neighbours of each pair nodes under consideration. For a target link between node i and node j , BCGP sums up the representation of their common neighbours obtained from the GNNs. This emphasizes the structural context and the shared neighbourhood, which are pivotal for predicting the existence of a link. Formally, the pairwise representation e i j for a potential link between node i and node j can be represented as: e i j = ∑ u ∈ N e i ∩ N e j GNN θ u , A , Z , (11) where θ represents the parameters of the GNN, N e ( i ) ∩ N e ( j ) denotes the set of common neighbours between node i and node j . The aggregated pairwise representation e i j is then used to compute the probability y ̂ i j of a link between node i and node j . This is achieved by passing e i j through a final prediction layer, such as a fully connected layer with a sigmoid activation function, which is denoted as follows: y ̂ i j = σ W n e i j + b n , (12) where σ denotes the sigmoid function, and W n and b n are learnable parameters of the prediction layer.

To train the model on the link prediction (LP) task, we use the binary cross-entropy loss, which is denoted as: L LP = − 1 | M | ∑ i , j ∈ M y i j ⁡ log y ̂ i j + 1 − y i j log 1 − y ̂ i j , (13) where M represents the set of all node pairs, which contains both positive (existing links) and negative (non-existing links) samples, | M | denotes the total number of node pairs in the set M , y i j is the ground truth label for the link between nodes i and j , where y i j = 1 if a link exists and y i j = 0 otherwise.

In this study, we focus on leveraging a comprehensive dataset to model lncRNA-miRNA and circRNA-miRNA associations. We have constructed a total of four datasets: lncRNA-miRNA association 1 (LMA1), circRNA-miRNA associations 1 (CMA1), lncRNA-miRNA association 2 (LMA2), and circRNA-miRNA associations 2 (CMA2). For the lncRNA-miRNA associations in LMA1, based on previous research (Wang Z. et al., 2023), we utilized the LncACTdb 3.0 (Wang P. et al., 2022) database. From this database, we extracted 1,057 experimentally verified lncRNA-miRNA associations, containing 284 lncRNAs and 520 miRNAs. Regarding the circRNA-miRNA associations in CMA1, building upon prior studies (Guo et al., 2024), we used CircBank (Liu et al., 2019) as the primary data source, obtaining 20,771 high-quality associations involving 3,802 circRNAs and 1,273 miRNAs. The lncRNA-miRNA associations in LMA2 were sourced from lncRNASNP V3.0 (Yang et al., 2023), comprising a total of 8,502 lncRNA-miRNA associations, including 467 lncRNAs and 254 miRNAs. As for the circRNA-miRNA associations in CMA2, they were derived from the dataset 1 of the KGANCDA (Lan et al., 2022), with a total of 702 circRNA-miRNA associations, encompassing 471 circRNAs and 439 miRNAs. The sequences of lncRNAs were sourced from LNCipedia (Volders et al., 2019) and NONCODE (Zhao et al., 2021), the sequences of circRNAs were obtained from CircBase (Glažar et al., 2014), and the sequences of miRNAs were acquired from miRBase (Griffiths-Jones et al., 2007).

To rigorously evaluate the performance of pre-training within our method, we compare several pre-training methods commonly used for initializing node embedding in GNN. We include two baselines, namely, the Random Embedding and the Adjacency Matrix Embedding methods which do not consider any RNA sequence information. The Random Embedding method initializes node embeddings with random values generated from a Gaussian distribution, while the Adjacency Matrix Embedding method leverages the adjacency matrix, representing the relationships between nodes in a lncRNA-miRNA-circRNA association graph, to generate node embeddings. Furthermore, we also compare three sequence-specific pre-training methods, including Text2vec (Mikolov et al., 2013), Doc2vec (Le and Mikolov, 2014) and HyenaDNA (Nguyen et al., 2024). Firstly, Text2vec employs a “bag of words” model for the k -mer representations of RNA sequence, converting each lncRNA, circRNA or miRNA sequence into a numerical representation based on k -mer occurrence frequency. Meanwhile, Doc2vec extends the Word2vec algorithm that generates document-level embeddings of the RNA sequences. Additionally, HyenaDNA leverages pre-trained HyenaDNA to generate the embeddings of RNA sequences at the nucleotide level, capturing the long-range dependencies within the RNA sequences.

The experimental results on the LMA1 and CMA1 datasets presented in Table 1, demonstrate the effectiveness of BCGP-BERT. For lncRNA-miRNA prediction, BCGP-BERT achieves the highest scores across all evaluation metrics, with an F1 score of 0.439, AUC of 0.904, AP of 0.435, and NDCG of 0.816, outperforming all other pre-training methods. Similarly, for circRNA-miRNA prediction, BCGP-BERT excels in all metrics, with an F1 score of 0.572, AUC of 0.948, AP of 0.712, and NDCG of 0.957. These results suggest that BCGP-BERT enhances association prediction by learning contextualized representations that capture the global context of entire sequences. Furthermore, the findings also highlight the effectiveness of using the BERT model trained on task-related data through Masked Language Modeling, enabling BCGP-BERT to discover the intricate pattern of RNA sequences and demonstrate its robustness in predicting RNA interactions. Additional pre-training results on the LMA2 and CMA2 datasets are provided in Supplementary Appendix S3.

To evaluate the performance of BCGP integrated with different fine-tuning methods, we conduct association prediction experiments using six different GNNs. Specifically, we compare GAT (Veličković et al., 2017), GATv2 (Brody et al., 2021), FiLM (Brockschmidt, 2020), GraphSAGE (Hamilton et al., 2017), SGC (Wu et al., 2019), and GCN (Kipf and Welling, 2016) using the BCGP-BERT pre-training framework. All methods share the same embedding dimension. As detailed in Table 2, for the lncRNA-miRNA association prediction, GCN outperforms all other fine-tuning methods, achieving an F1 score of 0.439, AUC of 0.904, AP of 0.435, and NDCG of 0.816. For the circRNA-miRNA association prediction, GCN also demonstrates a robust performance with an F1 score of 0.572, AUC of 0.948, AP of 0.712, and NDCG of 0.957, while SGC closely follows with competitive performance. These results demonstrate that GCN is the most effective fine-tuning method on BCGP-BERT for predicting lncRNA-miRNA and circRNA-miRNA interactions. Further fine-tuning results on the LMA2 and CMA2 datasets are provided in Supplementary Appendix S3.

Next, we compare our BCGP to the state-of-the-art SPGNN (Wang Z. et al., 2023), which utilizes the k -mer technique, Doc2Vec model and fine-tuning with GNN for RNA association prediction. Additionally, we include GCNFormer (Yao et al., 2024), which leverages graph convolutional networks and transformers for predicting lncRNA-disease associations, as a baseline. The experimental results in Table 3 illustrate that our BCGP consistently outperforms SPGNN across all metrics for both lncRNA-miRNA and circRNA-miRNA association predictions. For lncRNA-miRNA associations, our BCGP achieves an F1 score of 0.439, AUC of 0.903, AP of 0.435, and NDCG of 0.812, outperforming SPGNN’s scores of 0.430, 0.894, 0.419, and 0.828, respectively. Similarly, for circRNA-miRNA, BCGP obtains an F1 score of 0.572, AUC of 0.948, AP of 0.712, and NDCG of 0.957, significantly outperforming SPGNN. Additionally, we evaluate the effect of the NCN technique on the performance of BCGP. Compared to the BCGP without NCN, the integration of NCN shows improvements in both lncRNA-miRNA and circRNA-miRNA association predictions. These results demonstrate that BERT has a stronger capability to capture contextual information than the classic Doc2vec method, which allows our method to effectively capture complex and intricate relationships between RNA sequences, leading to more accurate predictions. Furthermore, integrating the NCN technique enables our method to learn more refined node representations, which ultimately enhances the prediction accuracy in complex RNA interactions.

In this section, we investigate the impacts of two important parameters of the BCGP, including the k -value in k -mers and the RNA embedding vector size in the pre-training stage.

Prediction of new circRNA-miRNA and lncRNA-miRNA interactions can reveal new biomarkers, identify therapeutic targets, and enhance understanding of the regulatory mechanisms of biological networks. To validate the practicality of our method, we select two miRNAs, namely hsa-miR-143 and hsa-miR-6808-5p, to verify the prediction results of miRNA-lncRNA and miRNA-circRNA associations generated by BCGP. External literature is employed to validate the predictions for miRNA-lncRNA associations, while associations recorded in CircBank (Liu et al., 2019) are used for verifying miRNA-circRNA associations. For the target miRNAs, 8 out of the top-10 predicted miRNA-lncRNA associations are confirmed in Pubmed, and 9 out of the top 10 miRNA-circRNA associations are validated in Circbank.

Using the BCGP model, we predict lncRNAs linked to hsa-miR-143. As illustrated in Table 4, the top-10 predicted lncRNAs are MALAT1, MEG3, NEAT1, UCA1, DANCR, HOTAIR, TUG1, GAS5, KCNQ1OT1, and MIAT, with 8 of these associations being validated in external literature. For instance, MALAT1 the top-ranked lncRNA, was shown in a study by Chen et al. (2017) to regulate ZEB1 expression by sponging miR-143-3p and promoting the progression of Hepatocellular Carcinoma. Additionally, Dong et al. (2020) has demonstrated that MEG3 overexpression inhibited LPS-induced injury in PDLCs by deactivating the AKT/IKK pathway by sponging miR-143-3p.

Similarly, we predict circRNAs that are potentially associated with hsa-miR-6808-5p validated by CircBank. The results presented in Table 5 reveal that the top-10 ranked circRNAs as hsa_circ_0082878, hsa_circ_0020316, hsa_circ_0049111, hsa_circ_0057955, hsa_circ_0037997, hsa_circ_0000726, hsa_circ_0049109, hsa_circ_0049112, hsa_circ_0016773, and hsa_circ_0085900. Upon searching CircBank for hsa-miR-6808-5p, we find 9 out of the top 10 predicted results in the CircBank dataset. The results of the case study indicate that BCGP possesses commendable practicality.