dataesr/scientific-paragraphs-categorization huggingface dataset:
RNA recognition motifs (RRM) are the most prevalent class of RNA binding domains in eukaryotes. Their RNA binding preferences have been investigated for almost two decades, and even though some RRM domains are now very well described, their RNA recognition code has remained elusive. An increasing number of experimental structures of RRM-RNA complexes has become available in recent years. Here, we perform an in-depth computational analysis to derive an RNA recognition code for canonical RRMs. We present and validate a computational scoring method to estimate the binding between an RRM and a single stranded RNA, based on structural data from a carefully curated multiple sequence alignment, which can predict RRM binding RNA sequence motifs based on the RRM protein sequence. Given the importance and prevalence of RRMs in humans and other species, this tool could help design RNA binding motifs with uses in medical or synthetic biology applications, leading towards the de novo design of RRMs with specific RNA recognition.
The interactions between proteins and RNAs are key to many different biological processes and crucial for the proper functioning of the cells. The RNA recognition motif is the most prevalent protein carrying such functions in eucaryotes, so understanding how this protein motif interacts with different RNA sequences is of immediate relevance. However, a general recognition code between this motif and the RNA is not known yet. We have performed a computational analysis to understand how the recognition process works, identify the main binding mode between the protein and the RNA and which are the most relevant amino acids involved in the recognition. This analysis allowed us to build a predictor that estimates the binding between an RNA recognition motif and any RNA sequence. We have named this method RRMScorer and it only needs the sequence of both the protein and the RNA to predict the results.
One of the main limitations of this method is that they rely on docking models and a detailed calculation of all the atomic interactions, and therefore have a strong dependence on the precise structural data that the potential was based on  [29]. The recent release of RoseTTAFoldNA  [30]  has also implied a huge advance in the field, providing high accuracy models with atomic resolution for protein/nucleic-acid complexes, which is extremely useful for proteins where it is clear which RNA the protein binds, but that is not always the case. Therefore, a fast and interpretable method that works at the sequence level and that is applicable to genome scale studies to identify possible interactions, or that can be used in computational screening in protein design, is not available. Here we introduce RRMScorer, a scoring method that overcomes these limitations by carefully combining information on protein and RNA, structure and sequence data, whilst minimising bias in both areas.
Q15717, P09012, P26599) recognize very different RNA sequences  [19]  whilst sharing the canonical binding mode. This meta-analysis of available RRM-RNA information brings us closer to understanding how this versatile motif works. RRMScorer provides a novel means to decipher a general recognition code for canonical RRMs by using a completely different approach from other RBP predictors, where complex deep-learning networks are trained using high-throughput RNA binding data to determine the consensus for different RRM domains  [21, 22]. Instead, our method relies on a carefully curated alignment for all the structural information available, that is then translated into a single score that estimates how likely is any residue in a specific position of the RRM to interact with any nucleotide of the RNA.
Source data To gather the RRM data available we started from the Pfam database  [31] , where 19 families containing RRMs were identified (S1 Table ) and further validated through visual inspection of representative structures for each family. All the available entries on the Protein Data Bank belonging to those families were retrieved and split into the different RRM domains. The structures can be found in proteins containing a single RRM domain to larger proteins containing up to four RRM domains. 1259 RRM domain structures were so identified in the complete RRM dataset (S1 Dataset), from which the domains in complex with RNA were extracted. 271 entries remained after removing the complex structures where less than three nucleotides are recognized (Fig  2 ). The complete RRM dataset and the RRM-RNA complex dataset (S2 Dataset) are available from the supplementary information folder in the Bitbucket repository.
All entries are identified by their UniProt code, RRM number, PDB Id. and chain, starting and ending positions of the RRM by PDB and UniProt, and UniProt starting and ending positions matching the sequence included in the file. The latter numbering is required to include some extra residues at both the C-and N-terminus regions that might still be relevant for RNA binding. From these RRM-RNA complexes, all protein-RNA interactions were computed using an in-house script and stored in a text file also available from the Bitbucket repository (S3 Dataset). An amino acid residue and nucleotide were considered to interact if any atom from the residue and the nucleotide were less than 5 Å away from each other. This is a broadly used interaction definition to keep strong interactions such as hydrogen bonds or electrostatic interactions, while still accounting for hydrophobic interactions that can occur at distances of 3.
Data cleaning and alignment procedure For the analysis, it was necessary to reduce the bias in the original complete RRM dataset. The protein sequences were extracted from the 1259 RRM domain structures (S1 Dataset) and to eliminate nearly identical RRM sequences in this set a sequence identity threshold of 99% was applied, after which 356 sequences remained in the reduced RRM dataset (S4 Dataset). This set still showed a strong bias towards certain RRM families, especially the RRM_1 Pfam family (PF00076), with 314 entries out of the 356 belonging to this family. To overcome this bias a set of 19 representative RRM domains was defined using CD-HIT  [33]  with a 30% sequence identity cut-off, and validating that we were selecting entries from different Pfam families, even though larger families such as RRM_1 are still repeated.
This representative RRM set (S2 Table ) contains very diverse RRM sequences for which the structure is available, and served as the core for generating the master multiple sequence alignment (MSA) (Fig  2). The sequence and structure information from the reduced RRM dataset was fed into PRO-MALS3D  [34]  to generate the alignment (Fig  2 ). Besides using structure and sequence, this tool employs secondary structure predictions to generate the MSA. For our study, we used the structures of the previously selected pool of 19 representative RRMs, and the sequences for the remaining 336 domains, totalling 356 RRM sequences. After manual checks to avoid gaps in the MSA for the 6 main secondary structures (β1-α1-β2-β3-α2-β4) due to some unusually long β-strands and/or α-helices, 347 RRMs remained in the clean alignment.
The clean alignment was further enhanced to improve the often poor alignments for the loops and terminal regions. Based on the principle that for amino acids in the loop regions the most important characteristic is how they are connected to the fixed secondary structure elements, these regions were 'squeezed' so all the gaps are placed in the middle of the loop regions, or at the extremes for the N and C terminal regions, in case of shorter loops or terminal regions, respectively (Fig  3 ). This alignment, hereinafter referred to as the master alignment, is then used to generate an HMM useful to quickly align other RRM sequences, such as the RRM-RNA structures dataset (Fig  2 ). The master alignment (S5 Dataset) and the.
RRM-RNA complexes similarity matrix To identify and cluster the entries where RNA molecules bind in a similar orientation to the RRMs, we developed a method to compute a pairwise similarity score based on the amino acid positions of the RRM that recognize the nucleotides of the RNAs (S1 Fig). This was calculated for the 271 RRM-RNA structures available (S7 Dataset). Some entries contain long RNAs (up to 3000 nucleotides), for which the RNA sequences were truncated to keep only the part that binds to the RRM. All RNA sequences were pairwise compared, by sliding them with respect to each other and checking whether their nucleotides bind similar amino acid sequence positions in the RRM master MSA. The number of matching positions between the two RNAs (positions of the RRM that both nucleotides are interacting with) are divided by the number of unique positions between both nucleotides, i.
Similarity score ¼ P i¼n i N Matching positions N Unique positions N Aligned nucleotides For each compared RNA pair, only the alignment with the highest similarity score is retained. A similarity matrix was then constructed from these best scores that includes all the RRM-RNA complexes, and which is later used to identify the different binding modes. The matrix is available as a CSV file from the supplementary information folder in the Bitbucket repository (S7 Dataset). To select a homogeneous cluster, we grouped all the entries that have a minimum score of 0. 25 with at least 25% of the complexes in the cluster. We explored the effect of changing the cutoffs on the cluster generation (S2 Fig) , and selected a combination of values that balances the variability and the similarity within the clusters, so allowing further meaningful analysis of the complexes.
The first cluster defined, hereinafter referred to as cluster 0, retained 187 entries and it was the largest found (Fig 2) using those cutoffs. Other smaller clusters were generated, but in this work, we will focus on the canonical binding mode represented by cluster 0. After visual inspection of some of the entries within cluster 0 we verified that the RNAs were bound similarly, and that we essentially captured the canonical binding mode  [1]. The entries raising the lowest similarity score were PDB Id. 6g90 (chain B-RRM1) and PDB Id. 3nnh (chain B-RRM1), with a similarity score of 0. 083. Based on their best alignment from the score calculation, a UC RNA fragment from 6g90 aligns with a GU RNA fragment from 3nnh. To verify they still shared a similar binding mode both complexes were superimposed (S3 Fig) and the RMSD was calculated using the heavy atoms in the sugar for the four aligned nucleotides, obtaining a value of 1.
RNA sequence alignment To conduct the RNA alignment, we used the same method as for the RNA binding modes identification, based on the amino acid sequence positions of the RRM that the nucleotides of the RNA are in contact with. When comparing two RNAs, the sliding window position that generates the highest score corresponds to the best possible alignment for those sequences. To align the 187 RNAs included in cluster 0, we selected the medoid, the entry with the highest similarity scores with respect to all other entries (PDB Id. 3hhn, chain D). For the 186 remaining entries we found the best alignment against the medoid and generated the RNA MSA. Gaps were added at the 5' or 3' ends of the RNA sequences when required to ensure that nucleotides in the same position were properly aligned. This is necessary due to the different length of the RNA fragments considered, which ranges from 3 to 11 nucleotides.
Following the previous example, this is the number of times any nucleotide except adenine is observed in RNA position 1 divided by the number of times an adenine is observed in that position. RRMScorer takes the information from the interactions observed in all available RRM-RNA complexes in cluster 0 (training set). However, selecting all available interactions for the scoring would bias this approach to the amino acids and bound RNAs for the most studied RRMs, which are overrepresented in our dataset. To generalize the approach as much as possible, we first selected a subset of nucleotide and residue positions that interact with each other in at least 20% of the RRM-RNA complexes having different UniProt identifiers in cluster 0, so limiting the analysis to key binding positions. These retrieve all the RRM positions already well known for its importance with respect to RNA binding, such as solvent-exposed RNP1 and RNP2 residues.
Validation with RNAcompete data RNAcompete is an in vitro technique to quickly analyse the RNA binding preference of RNA binding proteins (RBP)  [19]. The experimental RNA preferences for all the available RRMs are downloaded from the CISBP-RNA Database  [37] , resulting on information for 171 different proteins. The RNA preferences are presented as a matrix of frequencies observed for each nucleotide along 7 different positions (8 positions in a few particular cases). We reconstructed all the possible sequence combinations for those 7-mers and computed the associated bits the same way it is done for classic sequence logos  [38] , capturing the information content of each nucleotide in a particular position. As RRMScorer can score up to 5 nucleotide long RNA fragments, we keep the 5-mer with the highest average bits value from the RNAcompete 7-mers.
Considering that in the RNAcompete assay they used an RNA pool comprising ~240,000 short fragments (30-41 nucleotides) that guarantees that any 9-mer is at least repeated 16 times, we can assume that any 5-mer sequence is also present multiple times in the pool. Therefore, there are 1024 (4 5 ) theoretical different 5-mers, from which we can extract the associated bits value for each of the 171 RRM-containing proteins. The 171 proteins in the dataset are divided into 3 categories based on their RBP composition. This must be done to properly align the sequences to the master alignment and analyse the resulting scores for each individual domain. There are 47 proteins that consist of a single RRM, 118 have multiple RRMs and the 6 remaining ones have at least one RRM but in conjunction with other RBPs such as KH domains or zinc fingers.
The single RRMs are simply aligned to the master alignment (HMM) and the scores are calculated for all the possible RNA 5-mers (1024 sequences). For proteins with more than one RRM, we cut the sequence of each individual RRM and computed the score for the RNA 5-mers, we keep the highest score to correlate with the bits data. All the available data for each RRM domain, bits average value and RRMScorer predictions are available from the supplementary information folder in the Bitbucket repository (S10 Dataset).
Score confidence The confidence of the scores is calculated based on the final binding scores between the RRM and the RNA in the complex we are scoring, and how that compares with the scores from the training and randomized sets. The scores for both sets were fit into a gaussian kernel density estimator (KDE) using Scikit-learn  [39]. We then can compute the likelihood for the final binding score of an RRM-RNA complex to fall on either the training or the randomized regions. The ratio between those probabilities is calculated, normalized from 0 to 1 and provided as confidence score. Low values means that the RRM-RNA complex does not show very favourable contacts according to our dataset, which places it in the predominant randomized region. Correspondingly, high values (near 1) mean that the contacts (or absence thereof) were observed very often in our training set, and we can be confident about the prediction.
RRM alignment analysis and representation The RRM master alignment (S5 Dataset) was validated by checking the number of gaps in the main secondary structure elements, and comparing the sequence logos with prior information about the β3-strand and β1-strand, also known as ribonucleoprotein domain 1 (RNP1) and RNP2  [1] , respectively. The alignment for the RRM domains from the RRM-RNA complexes (S2 Dataset) was similarly validated, and the interactions previously extracted for all the complexes (S3 Dataset) were mapped into the alignment to analyse the RNA contacts along the protein structure. The frequency of the contacts and gaps is calculated as a percentage for each of the positions in the alignment of the 271 RRM-RNA complexes (Fig  4 ). As expected, most of the contacts occur in the β-strands, the 1st, 3rd and 5th loop, and the C-terminal region, which further validates the RRM-RNA alignment.
Internal validation We performed an internal validation by computing the scores for the entries in cluster 0, referred to as the training set, and comparing them with the scores from a set of randomly generated RNA sequences, referred to as the randomized set. The scores were calculated as explained in the score validation section on Materials and methods, taking out the entry from the dataset before computing the matrices with the scores. Following this procedure, we ensure that we are not biasing the scoring. The distribution of the scores for the training and randomized entries (187 entries in each set) are compared and plotted  (Fig 6). A separation between the sets is observed confirming that the scores can discriminate with certain confidence between the RNA sequences that would bind to an RRM from the ones that would likely not bind.
Validation with RNAcompete data Using the RNAcompete data we were able to evaluate the RRMScorer method with a large set of experimental RNA binding preferences. The data were processed to translate the RNA binding preferences into bits values that then we can compare with our predictions (S10 Dataset). RNA fragments with a higher average bits value should be better binders for that particular RRM, and consequently should correlate with higher scores from our method. The correlations between the scores and bits values are shown for the three different categories (Fig  7 ). While we can make a one-to-one connection between the nucleotide preferences and the single RRMs, in case of multiple RRMs, and especially when other RBPs are present, the observed preferences might be a combination of different specificities of each individual domain.
The authors determined the binding affinities between the RNA mutants and the MSI1 RRM1 using fluorescence polarization assays. The mouse variant was used on the experiments, which is identical to the RRM1 human variant. We then used RRMScorer but removed the MSI1 RRM1 entry from the training set to make the test-case as agnostic as possible. For this validation we used a window size of 3 nucleotides to calculate the scores, as the authors claim that this is the length that MSI1 specifically recognizes. We correlated the experimental results with the 3-mer fragment raising the highest score for each of the 36 RNA mutants (Fig 8 ). A table with the scores and affinity data is available from the supplementary information file (S3 Table ). The correlation proves that RRMScorer successfully distinguishes between high affinity RNAs containing the UAG core motif (bottom right corner) that are in the low nM range (40-200 nM), from the ones without the conserved three-nucleotide motif, whose affinities drop to the μM range (1000-2500 nM).
The guanine in the UAG motif is the most relevant nucleotide with respect to the binding affinity, when mutated to any other nucleotide the resulting interaction is on average 37-fold weaker than with the WT RNA (S3 Table ). Using as a reference the PDB entry 2rs2 for MSI1 RRM1, two residues are involved in the specific recognition of this guanine, K21 and F65, interacting with their respective sidechains (Fig  9A )  show a preference for any nucleotide (Fig  9C ). This was expected as p-stacking interactions are not nucleotide-specific. By computing the scores for all sliding windows in the RNA to find the highest possible score, we also predict which part of the RNA is more likely to bind the RRM. This is particularly useful to spot likely binding regions in long RNAs and the prediction confidence (Fig  10 ).
Validation with SRSF1 The complex of the human prototypical SR protein SRSF1 RRM1 with RNA was not available when we generated the master alignment. Therefore, it was not used to develop RRMScorer and it can be used as an external validation test case. The structure was recently released as part of a work by A. Clery et al. [24]  where the authors identified that this RRM has a strong preference to bind cytosines (PDB Id. 6hpj). They also engineered a variant to gain the ability to bind uridines, mutating a single residue in the β4 strand from glutamate to asparagine (E87N). We want to assess whether our method can separate the RNAs that bind from the ones that do not bind, and whether it can capture the difference in binding of the engineered variant. The closest RRM sequence in cluster 0 has 51% sequence identity (SRSF3 RRM1), but a visual inspection of the complex shows a similar binding mode, which is crucial for the reliability of RRMScorer.
91 to -0. 82 (Table  1 ), which is a significant change considering that only one residue is mutated from the 20 RRM positions used for the scoring. This score now falls in the likely-binder region (  the glutamate sidechain interacts with the amino group of the cytosine, while not being able to recognize a uracil in this position, as stated by the authors  [24]. The single point mutation E87N that enables this RRM to bind a uracil is in perfect agreement with our scores (Fig 11B ), in fact, an asparagine is the only residue showing a positive score for uracil (from the residues observed in our training dataset).
With our statistical approach we can score any RRM-RNA complex while easily tracking the individual scores of each residue-nucleotide pair. This information can be useful to rationally design new RRMs. Still, limited data availability is the main reason why our method is restricted to the canonical binding mode. We generated clusters for other binding modes but the number of complexes available was very limited, which makes it difficult to define well-grounded binding preferences. Finding validation sets was also challenging, because big changes on either the protein or the RNA sequence can completely change the RNA binding mode, so invalidating use of our method. From the protein structure side, the data availability is no longer a limitation after the AlphaFold Protein Structure Database release  [44].
Even though it does not solve the RNA recognition problem  [45] , current challenges purely based on protein structure, such as assessing the preferred RNA binding mode of an RRM, might be solved soon, although the current inability of such methods to cover dynamics and multiple conformations remains a bottleneck to be solved. The more generic approach we present here might therefore be more applicable at this point in time; by 'averaging' the limited available information we reduce complexity but enable interpretability. The amino acid representations extracted by unsupervised language models  [46]  might in this context also be able to provide generalisations of amino acid behaviour that are applicable for improved prediction of RRM-RNA binding. S7 Dataset. RNA binding similarity matrix for the 271 RRM-RNA complexes.
The corresponding positions in the RRM alignment for K21 and F65 are β1-1 and β3-3, respectively (Fig 1B and S6 Dataset), and the guanine from the UAG motif corresponds to the RNA position 4 in the RNA alignment (Fig 5 and S8 Dataset). In agreement with the experimental observations, a lysine in position β1-1 only shows a positive score for guanine (score of 0. 48), and negative scores when binding any other nucleotide (-0. 84, -1. 64 and -0. 58 for adenine, cytosine and uracil, respectively, Fig 9B). On the other hand, a phenylalanine in position β3-3 does not.