GlyinsRNA: a webserver for predicting glycosylation sites on small RNAs

ABSTRACT

Versatile RNA modifications play important roles in post-transcriptional regulations of gene expression, among which glycosylation modifications on small RNAs emerge as a novel clade whose characteristics need further interrogations. Here, we demonstrated that the sequence pattern around RNA glycosylation sites was not random and could be exploited for glycosylation site prediction. A machine learning predictor, GlyinsRNA, which integrated multiple RNA sequence representation encodings, was established. GlyinsRNA achieved AUROC (area under the receiver operating characteristic curve) of 0.7933 and 0.7979 in five-fold cross-validation and independent tests, respectively. GlyinsRNA was implemented as an online webserver, where both the predicted glycosylation sites and the overrepresented RNA-binding protein (RBP)-related motifs were annotated to facilitate the users. GlyinsRNA webserver is freely available at http://www.rnanut.net/glyinsrna.

KEYWORDS:

• RNA modification
• RNA glycosylation
• modification site prediction
• webserver
• machine learning
• RNA sequence analysis

Introduction

Various RNA modifications are extensively involved in the post-transcriptional gene expression regulation processes like RNA localization, RNA degradation, alternative splicing and translational efficiency regulation [1]. Although more than one hundred types of RNA modification had been discovered [1], just recently, the glycosylation modification has been shown to further extend the scope of RNA modification world [2]. Unlike most of the known types of RNA modifications, in glycosylation, instead of typical small chemical group (e.g. methyl- or acetyl-group), a much larger carbohydrate residue (e.g. sialic acid and fucose) was conjugated onto the RNA substrates. In the leading reports from Flynn et al., RNA glycosylation was mainly enriched on small RNAs and plays roles on RNA localization and RNA recognition, suggesting its functional importance in cell signalling [2]. However, the mechanism-of-action for RNA glycosylation remains to be fully elucidated, of which whether and how the sequence pattern contributes to this process has not been investigated.

In this study, we demonstrated that the sequence context around RNA glycosylation sites contains predictive information for RNA glycosylation, where both the positional sequence pattern and the RNA-binding protein (RBP)-related motifs are distributed in a non-random fashion. Based on these observations, GlyinsRNA, a sequence-based machine learning predictor of glycosylation sites on small RNA, was established.

Results and discussion

Non-random sequence pattern around RNA glycosylation sites

The sequence context around one site can be represented by a sequence window of 51 nucleotides (51 nt) with the modification sites at the centre. Figure 1(a) shows the positional sequence pattern of RNA glycosylation sites, where several positions, especially those proximal to the central glycosylation site (e.g. positions 0, +1, +2, +4 and +5), show non-random nucleotide propensities. On the other hand, however, no positional nucleotide propensity is observed for most of the other positions. One possible explanation is that there are position-independent sequence patterns (like short motifs) in addition to the observed positional sequence patterns. To test this possibility, we scanned the sequence windows with the pre-computed binding motifs of known RBPs. Indeed, the motifs of seven RBPs (RPS11, LIN28B, CPSF-100, IGF2BP1, FAM120A, DDX54 and NOP58) were significantly overrepresented in the sequence windows of the glycosylation sites in comparison with those of the non-modified sites (Figure 1(b), bar plots), and most of the results could be still observed with shorter window sizes until those of 37 nt (Supplementary Fig. S1), indicating plausible functional links between glycosylated RNA and these RBP regulators. Weak or moderate position preferences were observed when checking the positional distribution of few overrepresented motifs along the sequence windows (Figure 1(b), violin plots). For example, RPS11 and LIN28B have a tendency of binding to the upstream regions of glycosylation sites. Nonetheless, no strong position-biased distributions could be observed for any of these motifs. These results suggest the largely position-independent, short motif-like sequence patterns around glycosylation sites.

Figure 1. The positional and position-independent sequence pattern around the RNA glycosylation sites. (a) The positional sequence pattern around glycosylation sites. (b) The overrepresented RBP-related motifs found in the flanking sequence windows of glycosylation sites

Sequence-based prediction of RNA glycosylation site by machine learning models

Based on the above observations, we employed multiple sequence encoding schema to depict this non-random sequence pattern and trained random forest models based on these sequence feature encodings. We firstly optimized the window size k for each encoding. The performances of prediction models changing with the window size were summarized in Supplementary Fig. S2. The optimal window sizes of models by six encodings (PSNSP, PSDSP, One-hot, KNF, KSNPF, and PseKNC) were 51 nt, 51 nt, 51 nt, 17 nt, 19 nt and 37 nt, respectively.

Then, we constructed the models using the optimal window size for each encoding. As shown in Figure 2(a), the best-performed model in cross-validation is that based on the PSDSP encoding (AUROC = 0.7704). We further tested if combining multiple models could achieve better performance. As a result, a predictor combining three encodings (PSDSP, One-hot and KNF) was found to be the most competitive one (AUROC = 0.7933). Interestingly, One-hot [3,4] and KNF [5,6] encodings were designed to represent positional and position-independent sequence patterns, while PSDSP [7] encoding considered both aspects. Thus, the selection of the best encoding combination again suggests that both the positional and the position-independent sequence patterns provide predictive information for glycosylation sites.

Figure 2. Predictor performance evaluation results. (a) The ROC curve showing the performance of single and combined random forest models in five-fold cross-validation. (b) The ROC curve showing the performance of predictor on independent test dataset

To further validate the predictor performance, we tested the combined predictor on the independent test dataset. Samples in this dataset were not used for model training nor for model parameter and model combination selection and therefore could serve as a more rigorous and unbiased validation of the performance. On the independent test dataset (Figure 2(b)), the combined predictor (AUROC = 0.7979) also outperformed the best single-encoding-based model (AUROC = 0.7711). We also tried to evaluate the prediction model performance by using the precision–recall curve and found similar tendency (Supplementary Fig. S3).

GlyinsRNA: webserver for small RNA glycosylation site prediction

Based on the above performance evaluation results, we finally selected the combined predictor to build our RNA glycosylation site prediction tool, i.e. GlyinsRNA. GlyinsRNA has been implemented as an online webserver, which requires only RNA sequences as its input. Users could also select one stringent threshold according to their expected specificity level. Using a low-specificity threshold means to report more low-scored prediction results and therefore allow a higher tolerance to false-positive results. Here, we provided six alternatives for specificity level to users, including 99.9%, 99.5%, 99.0%, 95.0%, 90.0% and 85.0%. The primary output of the webserver is the position and score of each predicted glycosylation site. Besides, the sequence context around each predicted glycosylation site is also displayed, where both the glycosylation sites and overrepresented RBP-related motifs (if any) are highlighted by the distinct colours (Supplementary Fig. S4). We hope this easy-to-use webserver would facilitate further experimental investigations of small RNA glycosylation in the future.

Materials and methods

Dataset

The genomic coordinate annotation of glycosylation site was obtained from GEO (accession: GSE136967). The sites exclusively observed on the glycoRNA library (i.e. RNA library with glycosylation modification) were selected as the positive samples, while those only observed in the input library (i.e. background RNA library) but not in the glycoRNA library were assigned as the unmodified negative samples. Compared to the random negative sites, such selection could exclude unrelated sites and therefore is more rigorous. We primarily obtained 2720 positive and 5078 negative sites, which were randomly separated into training dataset and independent test dataset with the proportion of 7:3 (which is same for both the positive samples and the negative samples) to construct and assess the predictor. After mapping the sample sites to small RNA sequences, we employed CD-HIT-EST tool [8] to remove the redundant sequence samples with 90% identity threshold, by which 857 positive samples and 1068 negative samples were retained. Among all retained sample sequences, 114 positive and 186 negative samples were in the independent test set, while the rest were in the training set (with no overlap between the training and the testing samples).

Sequence pattern analysis

The sequence context around one site can be represented by a sequence window of k nucleotides with the modification/unmodified sites at the centre. The window size k was preliminarily optimized as 51 nt according to the cross-validation performance of the best-performed encodings. Meanwhile, we also explored the enrichment of RBP motifs with other window sizes and observed overrepresentation of most of the seven RBP motifs in shorter windows until those of 37 nt (Supplementary Fig. S1). If a sequence window extends beyond the border of an RNA, gap characteristics were filled to ensure all sequence windows were aligned to the centre. To investigate the positional sequence pattern of glycosylation sites, the aligned positive and negative sequence windows were analysed by TwoSampleLogo server [9]. Besides, we screened the significantly overrepresented RBP motifs in the positive sequence windows (in contrast to the negative sequence windows) by the AME tool of MEME Suite [10]. The Bonferroni adjusted P-value <0.05 was adopted as the significance threshold for the motif overrepresentation. The known RBP motifs were derived from the high-throughput CLIP-Seq experiment datasets, with the detailed description available in the previous work [11].

Predictor construction and evaluation

Here, the One-hot [3,4], PSNSP [12], PSDSP [7], KNF [5,6], KSPNF [4,13] and PseKNC [14] encodings were employed to depict the sequence pattern around glycosylation sites. These encodings are briefly introduced as following:

1) One-hot encoding translates four nucleotides, A, C, G, T, into binary numeric vectors, 0001,0010,0100,1000.

2) Position-specific nucleotide sequence profile (PSNSP) encoding displays the difference in frequencies of four nucleotides between positive samples and negative samples at each position of the sequence window.

3) Position-specific dinucleotide sequence profile (PSDSP) encoding is similar to PSNSP encoding, while it compares the frequency difference in dinucleotide between positive samples and negative samples.

4) K-nucleotide frequencies (KNF) encoding depicts the frequencies of continuous k-nucleotide on the sequence. The selection of K was preliminarily optimized as 2, 3 and 4.

5) K-spaced nucleotide pair frequencies (KSPNF) encoding also describes the frequencies of dinucleotide on the sequence. The dinucleotide is split by K-random nucleotides, and selection of K was preliminarily optimized as 0, 1, 2, 3 and 4.

6) Pseudo k-tuple nucleotide composition (PseKNC) encoding is the combination of the frequency and the physico-chemical property of dinucleotide or trinucleotide.

For each encoding, a random forest model was trained, and the prediction results from each were combined by the weighted summing schema where the weights of all models sum to 1. To determine the optimal window sizes of sequences, we trained models with various windows from 13 nt (larger than the RBP motif length which is 12 nt) to 51 nt with the step of 2 and selected optimal window sizes for each encoding. The best combination of random forest models was selected based on five-fold cross-validation performance. Besides, we also tested the models on the independent test set to avoid overestimation of performance. The performance was measured by typical area under the receiver operating characteristic curve (AUROC) metric. Finally, the predictor was implemented as an online webserver which was constructed under the ‘Linux + Apache + PHP’ framework.

Author contributions

CC established and tested the models. CC constructed the online webserver. CC and XW performed the sequence analysis. YZ and CC prepared and revised the manuscript. YZ supervised the study.

Acknowledgments

We thank Prof. Qinghua Cui, Yiran Zhou and Wanqing Zhao for their suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The RNA glycosylation site data from Flynn et al. is publicly available in public GEO series (accession: GSE136967) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136967here

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work has been supported by the National Natural Science Foundation of China [31801099 and 32070658 to Y.Z.].