The prediction of RNA–ligand binding sites is a binary classification task at the residue level. The model predicts the binding of each nucleotide (residue) in each RNA sequence to an organic small molecule or metal ion. Each residue is designated as binding or non-binding. A nucleotide is labelled as a binding site if the distance between any of its heavy atoms and any ligand heavy atoms is less than or equal to 4 Å, following the definition from previous research [12, 13, 16].

A total of six publicly available RNA–ligand complex structure datasets are utilized in this study: HARIBOSS, TR60, RB9, TL12, JL10, and TE18. Panei et al. [17] extracted RNA-small molecule complexes from the PDB to construct the HARIBOSS dataset (Accessed February 2025), composed of 862 structures. Subsequently, these complexes were clustered based on their sequence and structural similarity with RNA. Su et al. [11] also collected 712 RNA-small molecule crystal structures and used TM-scoreRNA to calculate their structural similarity, resulting in 78 representatives. The dataset is further segmented into two distinct sets: the TR60 set and the TE18 set. These sets are utilized for the training and testing of RNAsite. Gao et al. [12] generated RB9, TL12, and JL10 datasets to assess their binding site prediction program, ZeSTa. The RB9 originated from RB19 by removing ten structures from the dataset that overlapped with TR60. From the RNA-small molecule structures deposited after January 2021, JL10 is characterized by junction loop structures, which exhibit high structural complexity. In contrast, TL12 structures are devoid of the loop, resulting in low structural complexity. Among the six datasets, HARIBOSS and TR60 were merged to generate training, test, and validation sets. Non-biological small molecules commonly used as crystallization additives or cryoprotectants, including water, sulfate, phosphate, glycerol, ethylene glycol, and polyethylene glycol (HOH, SO₄, PO₄, GOL, EDO, PEG), were excluded from the ligand set to prevent artifactual contacts. RNA chains longer than 161 nucleotides were subsequently removed, resulting in a final dataset of 432 unique RNA chains. The dataset was segmented using a pre-determined random seed to create three subsets with a ratio of 8:1:1 for training, validation, and internal testing, respectively.

The HARIBOSS set was utilized to further evaluate the generalizability of the proposed model. Zhu et al. [13] provide the four structure split sets based on structural similarity with a TM-Score threshold of 0.5. Sixty models combining RNA LMs and loss functions were trained and tested using the same protocol as the combined sets of HARIBOSS and TR60.

The remaining four, RB9, JL10, TL12, and TE18, were reserved as test sets. To ensure independent evaluation, any overlapping sequences between the training and test sets were removed. The problem formulation and dataset design enable quantitative assessment of both prediction accuracy and generalizability under practically meaningful scenarios.

To predict whether each nucleotide in an RNA sequence is a binding site or not, a residue-level binary classification model was designed. Fig. 1 illustrates the model architecture of CoBRA. The program’s input is a query RNA sequence, which is embedded using pre-trained RNA LMs. During the training process, these embeddings are maintained as a frozen state and do not undergo updates. The prediction model employs an MLP architecture comprising five fully connected layers. Each layer is followed sequentially by layer normalization, a ReLU activation function, and dropout with a probability of .1, promoting stable training and improved generalization. The final output layer produces a two-dimensional logit for each residue, corresponding to the binding and non-binding classes. These are converted into probabilities via a SoftMax function. All inputs are standardized to a maximum sequence length of 161 residues. Sequences of shorter length are padded with zeros in the input matrix X and with −1 in the target vector y. A masking strategy is implemented to ensure that the padding regions do not contribute to the loss computation or affect model training.

The model training was executed using the PyTorch [18] framework. The optimization process was executed employing the AdamW [19] optimizer, with an initial learning rate set to 5 × 10⁻⁴. The cosine annealing method was employed for the purpose of learning rate scheduling. All experiments were run with a batch size of 4, a dropout probability of .1, for 100 epochs, and with the random seed fixed to 42 to ensure reproducibility.

Ten pre-trained RNA LMs were employed to pick the best model for residue-level binding site prediction. The models were pre-trained on diverse RNA types for various learning objectives, resulting in differences in embedding dimensionality and representational properties. Table 1 lists the RNA LMs, their pre-training targets, and the embedding size.

The models can be categorized into two groups based on pre-training sets. The first category is the models trained on non-coding RNAs (ncRNAs), composed of ERNIE-RNA [20], RiNALMo [21], RNABERT [22], RNA-FM [23], and RNA-MSM [24]. This group captures generalizable structural patterns across various ncRNA families. ERNIE-RNA [20], comprised of 12 Transformer blocks, is pre-trained via masked language modeling (MLM) on 20 million ncRNAs from RNAcentral with structural information. RiNALMo [21] is also pre-trained on 36 million ncRNAs, employing MLM with 33 Transformer blocks. RNABERT [22] adopts the pre-training BERT algorithm to 762 K ncRNAs. It also encodes the characteristics of the RNA family and structure. RNA-FM [23] is built upon 12 bidirectional Transformer encoder blocks that produce an L × 640 embedding matrix for input length L, pre-trained on 23.7 million ncRNAs. RNA-MSM [24] follows the MSA Transformer architecture with ten blocks to learn two-dimensional co-evolutionary signals from homologous multiple sequence alignments of 3932 Rfam families.

The members of the second class are pre-trained on pre-mRNA or mRNA untranslated (UTR) regions. SpliceBERT, SpliceBERT-510, SpliceBERT-H510 [25], UTRLM-MRL, and UTRLM-TE_EL [26] are members of the class. The category specializes in modeling sequence features in post-transcriptional regulatory regions. SpliceBERT [25] is a BERT-based model with six Transformer encoder layers pre-trained on 2 million RNA sequences from 72 vertebrate species. Additionally, two variants of the model were also employed: SpliceBERT-510, an intermediate checkpoint, and SpliceBERT-H510, pre-trained exclusively on human data. UTR-LM [26] is a six-layer Transformer encoder pre-trained via semi-supervised masked nucleotide reconstruction, 5′ UTR secondary-structure prediction, and minimum-free-energy regression on 5′ UTRs from multiple species. Two task-specific variants, UTR-LM–TE_EL (translation efficiency and mRNA expression-level prediction) and UTR-LM–MRL (mean ribosome loading prediction) were also used to generate embeddings.

All language models were used in their pre-trained form with frozen parameters; no fine-tuning or weight updates were performed during training. Each RNA sequence was transformed into a residue-level embedding sequence using the corresponding model, then padded to a fixed length of 161 residues before being fed into the MLP classifier. Embedding dimensionality varied across models, ranging from 120 to 1280. The predictive performance of different embedding types was compared to assess their impact on RNA–ligand binding site prediction. Additionally, a series of models was constructed and tested to connect the two LMs by concatenating their embeddings. Furthermore, one-hot encoding of bases was employed in order to provide a baseline.

Compared to the RNA sequence length, the proportion of binding nucleotides (positive class) is low, leading to a severe class imbalance problem. To resolve the issue and quantitatively evaluate how different training objectives affect performance, we employed six loss functions.

Each loss function represents a different optimization strategy. Binary cross-entropy (BCE) serves as the conventional baseline for binary classification, encouraging predicted probabilities to converge to the ground truth. Class-balanced focal loss [27] and Tversky loss [28] are designed to address the class imbalance. The class-balanced focal loss emphasizes hard-to-classify samples by assigning them higher weights. In contrast, the Tversky loss applies asymmetric weighting to false positives and false negatives. By giving different penalties, the loss function tries to improve model sensitivity for an imbalanced dataset. Dice loss [29] and Lovász hinge loss [30] focus on residue-level spatial alignment and structural consistency. Dice loss maximizes the overlap between predicted and true positive regions, whereas Lovász hinge loss optimizes the Intersection over Union directly.

Lastly, we introduce a composite loss that combines triplet center loss (TCL) [31] and class-balanced focal loss to simultaneously improve class separation in the embedding space and address the class imbalance. This design is inspired by Wang et al. [32], demonstrated strong performance in protein-small molecule binding site prediction using a combination of TCL and class-balanced focal loss to enhance both feature discrimination and class imbalance handling. The total loss is defined as shown in Equation 1.

L _focal and L_tcl are class-balanced focal loss and TCL loss, respectively. λ is used as a weight to balance between the losses. L_focal, a class-balanced focal loss, assigns higher weights to address class imbalance, as shown in Equation (2).

p_i^t denotes the predicted probability of the true class t, and (1-p_i^t)^g serves as a modulating factor that down-weights well-classified samples and focuses the training on hard examples. The focusing parameter g controls the relative importance of easy and hard samples. The term E_n, represents the effective number of samples for each class and balances the impact of class frequency.

TCL learns a center vector for each class and encourages embeddings of the same class to cluster together while enforcing a margin-based separation between different classes (Equation 3).

Here, f is the embedding vector of the input, c_y and c_1 − y denote the center vectors of the true and opposite classes, respectively, and m is the margin. The hyperparameters were set as g = 3, b = 0.999, l = 0.2, and m = 4, following the optimal configuration reported by Wang et al. [32].

All metrics are calculated based on the values of true positive (TP), false positive (FP), true negative (TN), and false negative (FN) in the confusion matrix. Precision is the ratio of TPs to all predicted positives (TP + FP), while recall is the percentage of TPs to all actual positives (TP + FN). F1-score is a harmonized mean of precision and recall that evaluates the balance between the two metrics. The Matthews correlation coefficient (MCC), a class imbalance robust metric, includes all elements of the confusion matrix, calculated as follows:

The value of MCC can range from +1 (perfect prediction) to −1 (perfect mismatch).

The area under the receiver operating characteristic curve (AUROC) quantifies the model’s discrimination ability by measuring the area under the ROC curve (TP rate versus FP rate), and area under the precision–recall curve (AUPRC) measures the area under the precision–recall curve, summarizing the trade-off between precision and recall across all thresholds.

To analyze and characterize the three-dimensional RNA structures quantitatively, the Laplacian-based curvature descriptor (LN) [12] is employed. This metric quantifies the curvature at each nucleotide by measuring the deviation of its C3′ atom coordinates from the relative positions of surrounding residues. It thereby captures local structural distortion and reflects how naturally a nucleotide is embedded within the overall structure.

For each nucleotide, a coordinate vector was defined by the C3′ atom position. Then a Gaussian kernel-based weighting function was constructed from the all residue pairwise Euclidean distance matrix following Equation 5.

Here d is a single-scale parameter, set as the median of all pairwise distances between nucleotides in the structure. p_i and p_j denote the Euclidean coordinates of the C3′ atoms of nucleotides 𝑖 and 𝑗. This choice of scale balances local and global shape sensitivity, avoiding excessive localization or over-smoothing.

Using the weight function, the Laplacian norm value LN_i for nucleotide i is defined as Equation 6.

This value reflects how far nucleotide i lies from its Gaussian-weighted center, capturing the relative positional characteristics within the RNA structure. Higher LN values indicate convex or protruding regions, while lower values suggest concave or densely packed areas. The computed LN values were utilized for analysis and structural visualization.

To evaluate the computational cost, all sequences were embedded using the frozen RiNALMo model on an NVIDIA RTX A6000 GPU. The resulting 1280-dimensional embeddings were used to compare the cost of logistic regression (LR) and random forest (RF). The hyperparameter settings for LR and RF are provided in Supplementary Table S1.

A total of 60 models were evaluated and compared, incorporating 10 RNA LM embeddings and 6 loss functions. Table 2 lists the top 10 MCC models and the best one-hot encoding model, benchmarked on the test set. Supplementary Table S2 showed the performance of individual models.

All the top 10 models demonstrated superior performance to the baseline model (MCC: 0.099, AUROC: 0.568), confirming that the performance improvements of CoBRA primarily originate from the contextual representations learned by RNA language models. The highest MCC is obtained by combining ERNIE-RNA with TCL focal as the language model and loss function, respectively. The model has also been determined to be second in terms of AUPRC and precision among the 60 models considered. It was observed that ERNIE-RNA, RNA-FM, and RiNALMo emerged as top contenders, appearing in the top ten MCC models on three separate occasions. The three LMs were pre-trained from ncRNAs. Conversely, the LMs pre-trained on mRNA-UTR regions are not included among the top 10, except for SpliceBERT-510. The comparatively diminished performance of mRNA-pretrained models might be due to their constrained exposure to a variety of RNA structures, which could result in diminished generalizability in binding site prediction when compared to ncRNA-pretrained counterparts. A survey of the top 10 MCC models reveals that 9 of them utilize entropy-like loss functions, including focal, TCL focal, and BCE. Except for the three aforementioned functions, the Lovasz hinge with RiNALMo is ranked sixth in terms of MCC.

The distribution of MCC, AUPRC, and AUROC for all 60 models is illustrated in Fig. 2 and Supplementary Table S3 across the LMs and loss functions. As observed by the top 10 MCC models, RiNALMo, ERNIE-RNA, and RNA-FM, which were pre-trained on ncRNAs, exhibited the highest average MCC, AUROC, AUPRC, and precision values (Fig. 2A). The five LMs, trained on mRNA-UTR regions, demonstrated consistent performance, with MCC ranging from 0.292 to 0.346. Conversely, the performances of LMs trained on ncRNAs exhibit significant variability. As observed in the top 10 models, ERNIE-RNA, RiNALMo, and RNA-FM, which were trained on more than 20 million sequences, exhibited superior performance to mRNA-UTR models, while RNABERT and RNA-MSM showed lower MCC values than UTR models. In terms of MCC, the best model using RNABERT (0.331) performed worse than the worst model using RiNALMo (0.372).

With regard to loss functions, models that have been trained with BCE and focal loss have demonstrated a consistent and superior performance across key evaluation metrics in comparison to other loss functions. Conversely, models employing dice loss and Lovasz hinge loss exhibited unstable convergence and suboptimal performance in most cases (Fig. 2B).

The LMs with the top three MCC values, ERNIE-RNA, RNA-FM, and RiNALMo, were subsequently selected for the concatenation experiment (Supplementary Table S4). The combination of ERNIE-RNA and RNA-FM yielded the highest MCC among the concatenated models. However, none of the models performed MCC higher than the best model of single LM embedding. This implies that the embeddings might contain analogous information, rendering their combination ineffective.

Based on the single-embedding models from the main experiments, a selection of the top five MCC configurations was made, and subsequent benchmarking was conducted on the external validation sets.

To assess the generalizability of the models, the sixty models were retrained and evaluated using the RNABind structure-based split dataset. This ensures strict non-redundancy between the training and the test sets at both levels. The results of the top five models are given in Supplementary Table S5. Compared to the methods that combine sequence information as LMs and structural information, CoBRA demonstrated lower AUROC values ranging from 0.605 to 0.657, while RNABind, RLBind, Ret, and MVRBind reported AUROC values ranging from 0.671 to 0.776. On the other hand, CoBRA and RNABind demonstrated similar performance in terms of AUPRC.

We conducted a benchmarking study of the top five MCC models on various test sets composed of RNA-compound crystal structures. RB9, JL10, TL12, and TE18. Table 3 presents the performance of the models and other RNA-compound binding site prediction methods. The CoBRA models are designated CoBRA-M1, M2, and so forth, following their MCC ranks.

While the CoBRA-M1, ERNIE-RNA combined with TCL focal model demonstrated the highest MCC, AUPRC, and precision among the five models evaluated on the test set, its performance was not optimal on the four benchmark sets. Conversely, CoBRA-M3, utilizing RiNALMo as the RNA LM and BCE as the loss function, exhibited the highest performance in MCC and AUROC on average of the four test sets. For precision and recall, M1 and M2 performed best, respectively.

With regard to the MCC from individual benchmark sets, M2 has the highest value for JL10 and RB9, whereas M1 shows the best MCC for TL12. While M3 is only ranked top for TE18 among the five CoBRA models, it is ranked second for JL10, a set with a high structural complexity, and TL12, a set with a low structural complexity. This suggests that M3 has greater generalizability for RNA-compound binding site prediction problems than other models. Consequently, M3 is selected as a representative model for CoBRA, and the results are subjected to further analysis.

A comparative analysis of RNA-compound binding site prediction methods revealed that CoBRA exhibited superior performance, with the exception of TE18 (Table 4). On average, CoBRA exhibited a 15.4%, 9.2%, and 33.8% improvement in MCC, AUROC, and recall, respectively, when compared to ZeSTa, a state-of-the-art program with structural information. In particular, on datasets with relatively simpler structures, such as RB9 and TL12, CoBRA also outperformed other methods, with the exception of RB9 precision and AUROC. In a similar vein, the CoBRA algorithm demonstrated the most optimal performance in all metrics for the highly complex structure set, JL10. Conversely, the TE18 dataset exhibited a decline in performance, with an MCC reduction of 41.90% in comparison to ZeSTa.

One of the successful cases of CoBRA is illustrated in Fig. 3A. A structure of pir-miRNA-300 apical loop fused to ydaO riboswitch scaffold, complexed with c-di-AMP (PDB ID: 6WTR) from the JL10 dataset. The switch regulates the gene expression in bacteria by sensing concentrations of ATP and c-di-AMP. The RNA structure is comprised of two three-way junctions connected by two helices and a large conserved interior loop. The compound works as an off switch when binding to the RNA. The binding site of the molecule is located between the two pseudoknots, stabilizing the RNA structure through intermolecular stacking [33, 34]. A comparison of the two programs revealed that CoBRA exhibited superior performance in comparison with ZeSTa. The MCC, precision, and recall of CoBRA are 0.923, 0.946, and 0.946, respectively, while those of ZeSTa are 0.121, 0.455, and 0.142, respectively. The uniqueness of the binding site characteristics makes the structure-based binding site prediction methods hard to predict, while CoBRA, a sequence-based one, demonstrated superior performance.

The computational cost of CoBRA was then compared with that of the LR and RF models. The generation of the embedding using RiNALMo required a time period of 111 ms per sequence. Utilizing fixed embeddings, the MLP required 109.97 s for training and 0.222 s for testing, while the LR necessitated 44.34 s and 0.008 s, and the RF employed 74.79 s and 0.101 s.

An investigation was conducted to examine the outcomes of TE18, wherein CoBRA exhibited the least optimal performance among the test sets. Dissecting the spatial structure of RNA (DSSR) [35] was employed to analyze the secondary structures of RNA crystal structures. The RNA secondary structures were classified into eight categories: stem, canonical, isolated canonical, internal loop, hairpin loop, bulge, junction, and helix. A subsequent analysis of the prediction performance of CoBRA was conducted, with the structural type serving as the primary variable.

CoBRA demonstrated a consistent level of accuracy, with a minimum of 51% for identifying binding sites across structural types. The highest accuracy was observed in junction sites, with an accuracy of 92%. On the other hand, the performance for identifying non-binding nucleotides exhibited variability according to structure, with higher accuracy in the stem (68%) and junction (67%), but lower accuracy in internal loops (47%) (Supplementary Fig. S1).

We also conducted a residue-level curvature analysis of the TE18 dataset by calculating the LN. As the LN values increase, the geometry of the residue evolves from concave to convex. The LN values of the binding and non-binding sites show distinct distribution (Supplementary Fig. S2), suggesting structural differences between the two classes. The binding site nucleotides have LN values ranging from 1.5 to 20.7, with an average of 10.9. Conversely, the non-binding ones exhibit larger values, ranging from 2.3 to 23.0, with an average of 12.7. It can be inferred that the geometry of RNA compound binding sites is relatively concave, a finding that is also reported by Su et al. [11].

Despite the absence of explicit incorporation of structural characteristics such as concavity as an input feature in CoBRA, a statistically significant relationship was observed between the model predictions and the LN value, as indicated by a point-biserial correlation of −0.171 with a P-value of 2.8e-5. This demonstrates that a model trained solely on RNA sequence embeddings could explain structural characteristics via inherent sequence patterns. We also observed a correlation of −0.247 between LN value and the actual binding label, indicating that actual binding sites tend to be located in structurally concave regions. Correlation analyses with TP/TN/FP/FN from CoBRA prediction demonstrate that TP and TN exhibit relatively high correlation coefficients (TP: r = −0.188, P = 4.0e-6; TN: r = 0.260, P = 1.2e-10). No significant relationship was found for FP (r = −0.027, P = 0.505), and FN exhibited a modest negative correlation (r = −0.126, P = 2.0e-3), indicating a slight reduction in missed positive predictions. A notable example is provided in Fig. 3B (PDB ID: 1FMN), illustrating the substantial correlation between CoBRA prediction and LN values. The flavin binding site of an RNA aptamer, defined by 4 Å from the co-crystalized ligand, exhibited a − 0.541 correlation coefficient (P = 8.0e-4). The right panel displays the LN values of nucleotides, indicated by the color change from blue (concave) to red (convex). CoBRA demonstrates an accuracy in predicting the binding site, with a precision of 0.636. However, the program failed to accurately predict three nucleotides with low LN values (average: 8.36) and erroneously predicted four nucleotides with high LN values (average: 9.89). This finding might show a limitation of only using RNA sequence information.

One of the most severe predicted cases of CoBRA in the TE18 dataset is sisomicin complexed with the bacterial ribosomal decoding site (PDB ID: 4F8U, Fig. 3C). The crystal structure contains homodimer RNA chains, and two sisomicins bind to the interface of the dimer. For prediction, the RNA sequence of a single chain, B chain of the PDB structure, was used as the input. CoBRA demonstrated a successful prediction for one of the binding sites (the red-colored region of the left panel). However, given that the binding site is designated by the ligand with the same chain ID of the input RNA sequence (the blue region), the resulting accuracy for that particular sample was found to be 0.0. A visual inspection of both the A and B chains together (right panel) revealed that the region was misclassified as a false positive. This case demonstrates that the accuracy of prediction may be diminished when binding sites are situated at the interface with other chains of a query sequence. This phenomenon has also been observed in protein-ligand binding site prediction problems using LMs [36]. One potential solution to this issue involves the inference of structural information or the incorporation of multi-chain as an input.

Metal ions such as Mg²⁺ and Na⁺ often bind diffusely across the RNA surface to neutralize the negatively charged phosphate backbone, without forming specific binding pockets. These ions contribute to RNA structural stability and folding, as previously noted by Draper et al. [37]. In contrast, organic molecules typically bind within well-defined pockets, making their spatial localization and prediction of binding sites more tractable. Thus, predicting metal binding sites may be more challenging than predicting organic molecule binding sites due to the nature of metal binding. Also, due to the electrostatic and steric nature of metal ion binding, they may associate with multiple structurally similar regions with comparable affinity, complicating the prediction task, which is reported in MetalionRNA [38].

The binding sites of the test sets were classified according to the type of binding molecule: metal ion or non-metallic compound. The accuracies of metal binding site prediction were 0.647, 0.531, 0.497, and 0.202 for RB9, JL10, TL12, and TE18, respectively. As anticipated, the non-metallic compound category demonstrated higher levels of accuracy, with the values obtained for RB9, JL10, TL12, and TE18 being 0.733, 0.739, 0.690, and 0.468, respectively. Detailed information is summarized in Supplementary Table S6. This discrepancy is likely attributable to the inherent nature of RNA–metal ion interactions. The incorporation of physico-chemical characteristics has the potential to enhance the efficacy of predicting metal binding sites.