RNA structure drives interaction with proteins

1Since the central dogma was proposed in 1950, the main role attributed to RNA has been to act as the intermediate between DNA and protein synthesis. Yet, more than 70% of the genome is transcribed and just a small part codes for proteins,, which indicates that the majority of RNAs could have different biological roles. During the past decade many efforts were made to develop methods to study RNA isoforms: sequencing has been essential for detection of RNA species and recent developments provided a great deal of data on polymorphisms, expression and half-lives of all types of RNAs, generating a valuable resource to understand their cellular functions and regulation. Although a number of techniques identified biological characteristics such as cellular location and secondary structure,, the characterization of the interaction network remains one of the most urgent challenges,. To this aim, computational methods are being developed to identify physicochemical features of the transcripts, their conservation between species and, most importantly, binding partners that are also active in the cellular environment.

2RNA is involved in many cellular processes such as control of gene expression, catalysis of various substrates, scaffolding of complex assemblies, and molecular chaperoning. Its ability to act as a hub of cellular networks is at the centre of an active research field and has already led to the discovery of diverse ribonucleoprotein (RNP) assemblies,. A number of membrane-less organelles contain specific mixtures of RNAs and RBPs (RNA-binding proteins) that, due to their intrinsic lability, are difficult to characterize. In most cases, liquid-like RNP assemblies, or condensates, such as P-bodies and stress granules, exchange components with the surrounding content and adapt to the environmental condition in a dynamic way. Within these phase-separated assemblies RNA plays a central role: whereas a polypeptide of 100 amino acids can interact with one or two proteins, a chain of 100 nucleotides is able to bind to 5–20 proteins, thus providing an ideal platform or scaffold for interactions,. Not surprisingly, changes in the interactions within RNP granules leading to liquid-to-solid phase transition are associated with the development of several human diseases, including neurological disorders and different types of cancer. In RNP condensates such as stress granules, regulation of protein and RNA contacts is primarily controlled by HSP70 and co-chaperones that act as versatile elements promoting assembly and disassembly of complexes.

3In this large spectrum of activities, RNA structure controls the precise binding of proteins by creating spatial patterns and alternative conformations for the interactions to occur. Known complexes in which the structure regulates protein binding include transfer RNAs (tRNAs) whose three-dimensional conformation facilitates the codon/anticodon interaction and the ribosomal RNA (rRNA) scaffold that sustains the ribosome. Importantly, the structure of a messenger RNA (mRNA) defines its lifecycle, recruitment of ribosomes and response against environmental changes. There are several cases of nucleotide chains of non-coding RNAs acting as scaffolds for protein complexes: structured domains in NEAT1 attract paraspeckle components and repeat regions in XIST sequester proteins to orchestrate X-chromosome inactivation. By contrast, poorly structured snoRNAs have been shown to facilitate the assembly of other transcripts.

4Of both coding and non-coding transcripts, RBPs are known as the major regulators and are classified as single-stranded RNA (ssRNA) and double-stranded RNA (dsRNA), depending on their binding preferences. Here we investigated the relationship between RNA structure and ability to interact with RBPs. At the transcriptome level, we find that the amount of RNA secondary structure correlates with the number of protein interactions. We propose several possible implications of this relationship: a link to RNA types and biological roles; a connection to regulatory networks; and the ability to modulate phase separation. Based on our observations, we also demonstrated that this RNA property can be exploited in vitro to tune the contact network of a protein aggregate.

5With the aim of studying how RNA structure influences protein binding, we measured the amount of double-stranded regions of the human transcriptome (Fig. 1a). We first grouped the RNAs, as detected by enhanced crosslinking and immunoprecipitation (eCLIP) approach, in classes based on the structural content measured by ‘parallel analysis of RNA structure’ (PARS) (Supplementary Fig. 1a and Fig. 1b). PARS is an experimental technique that distinguishes double- and single-stranded regions of RNA using the catalytic activity of two enzymes, RNase V1 (able to cut double-stranded nucleotides) and S1 (able to cut single-stranded nucleotides) and for which positive scores indicate double-stranded regions (see Eq. (1) in Methods). We then used catRAPID predictions of protein–RNA interactions (available from the RNAct database that contains both proteome-wide and transcriptome-wide calculations) and compared the interaction scores of different groups (HS, high structural content, vs. LS, low structural content) (Fig. 1b). The catRAPID algorithm estimates the binding potential through van der Waals, hydrogen bonding and secondary structure propensities of both protein and RNA sequences (total of 10 properties), allowing identification of binding partners with high confidence. Indeed, as reported in a recent analysis of about half a million of experimentally validated interactions, the algorithm is able to separate interacting vs. non-interacting pairs with an area under the curve (AUC) receiver operating characteristic (ROC) curve of 0.78 (with false discovery rate (FDR) significantly below 0.25 when the Z-score values are >2). Comparison of RNA groups with different structural content shows a consistent trend in which higher structural content in RNA molecules results in higher protein interaction scores (Fig. 1b). As for the PARS data, we note that the amount of double-stranded regions correlates weakly (<0.10; Pearson’s) with RNA length and GC content, indicating that these two factors positively contribute to the secondary structure by increasing the size of the conformational space as well as the overall stability.

6We repeated the analysis with a unrelated approach, RPISeq, which predicts protein–RNA interactions using sequence patterns in nucleotide and amino acid sequences. RPISeq is comprised of two methods based on support vector machines (RPISeq-SVM) and random forest (RPISeq-RF). Due to specific computational requirements, we applied RPISeq to an ensemble of RBPs (50 proteins with sequence similarity <0.85; http://cd-hit.org/) against the HS and LS set from the tails of the structural content distribution (100 transcripts) to estimate the binding probabilities (Supplementary Data 1). In both cases, the HS set (RF 0.80, SVM 0.71) is predicted to bind with significantly higher probabilities than the LS set (RF 0.70, SVM 0.54; p value <10–5; Kolmogorov–Smirnov (KS) test; Supplementary Fig. 1b–c), in agreement with catRAPID analysis (Fig. 1b). Thus, our analysis suggests that the RNA structure content has effect on the interaction with proteins.

7To match our predictions with experimental data, we investigated all RBP–RNA interactions revealed by enhanced CrossLinking and ImmunoPrecipitation, eCLIP (118 RBPs; see Methods). eCLIP provides protein contacts on target RNAs at individual nucleotide resolution through ligation of barcoded single-stranded DNA adapters. In agreement with catRAPID predictions (Fig. 1b), eCLIP binding scores correlate with PARS secondary structure, which indicates that the RNA propensity to interact with proteins is proportional to the amount of structure measured transcriptome wide (Fig. 1c). We note that the CLIP-seq approaches in general favour detection of single-stranded (SS) RNA at the expense of double-stranded (DS) RNA and the eCLIP dataset is not enriched in double-stranded RNA-binding proteins (9 out of 118 are assigned according to UniProt as dsRNA binding, 12 out of 118 as ssRNA binding, using available GO annotations), which indicates that our results are not biased by the protein types used in our analysis.

8To further corroborate that the trend is genuine and not only intrinsic to PARS measurements, we analysed the protein-interacting potential of the entire human transcriptome against the RNA secondary structure measured with the dimethyl sulfate modification (DMS) technique (differently from PARS, high values indicate single-stranded regions; Fig. 1d). This method of assessing RNA structure employs deep sequencing to detect unpaired adenosine and cytidine nucleotides. Once more, the analysis shows that the RNA secondary structure of the human transcripts is tightly correlated with protein-binding abilities.

9We also used the POSTAR database (containing >1000 CLIP-seq datasets; http://lulab.life.tsinghua.edu.cn/postar/) to retrieve the RNA-binding preferences of human proteins (103 experiments, 85 different RBPs) measured with PAR-CLIP, high-throughput sequencing-CLIP (HITS-CLIP) and individual nucleotide resolution CLIP (iCLIP). Due to intrinsic differences in the CLIP approaches (and other factors, such as the cell lines employed), each experiment reports different protein–RNA interactions. Yet 77% of the RBPs have preference for highly structured RNAs for at least one of the experimental methods (DMS or PARS; Fig. 1e).

10Given possible technical biases of high-throughput experiments, we decided to verify the reproducibility of the trend by investigating the correlation between RNA structure and protein interactions in low-throughput analyses. We first studied the interactome of eight large (>1000 nt) RNAs whose protein partners have been identified by microarray, a crosslinking-free approach,, (see Methods). In parallel, we estimated the structural content of each transcript using the CROSS algorithm that was previously trained on SHAPE data to predict the double-stranded propensity at nucleotide level resolution. Our results presented in Fig. 1f indicate that highly structured transcripts have more protein contacts than poorly structured transcripts, which is fully compatible with the findings presented in our previous analysis (Fig. 1b–e).

11We corroborated our observations through the study of RNP complexes deposited in the Protein Data Bank (PDB) database (X-ray resolution <2 Å; Supplementary Data 2; see Methods), which is comprised of 196 distinct RNA–protein pairs (>20 species) analysed with different techniques (mainly X-ray and nuclear magnetic resonance (NMR)) by different laboratories. Measuring the amount of RNA intra-contact (i.e. amount of RNA structure) and inter-contact (i.e. amino acid) per nucleotide chain, we found a striking correlation of 0.78 between the two variables, which provides compelling evidence of their tight relation (Fig. 1g; see Eqs. (2) and (3) in Methods).

12Thus, independently of the experiment (PARS, DMS, microarray, X-ray, NMR, eCLIP, PAR-CLIP, HITS-CLIP and iCLIP), the algorithms employed (catRAPID and RPISeq or CROSS to mimic SHAPE data) or organism (PDB database), we found a correlation between number of protein interactions and RNA structural content.

13We next investigated if the tight link between secondary structure and number of protein interactions is a property of specific RNA types (Fig. 2a). To this aim, we compared the secondary structure and protein interactions of transcripts ranked by sequence similarity using the CD-HIT algorithm (http://cd-hit.org/). With a threshold of 85% similarity, we found 22 clusters (total of 55 transcripts) with at least one RBP contact revealed by eCLIP. We then calculated the correlation between DMS signal and eCLIP protein interactions for each cluster and obtained a negative correlation in 64% of cases. This finding indicates that between two similar transcripts the one with higher structural content is more likely to have a larger number of protein interactions.

14The two transcripts sharing the highest similarity (99.31%) are the γ-globins HBG1 and HBG2 (haemoglobin subunits γ1 and γ2) that are expressed in fetal liver, spleen and bone marrow (NCBI Gene ID: 3048). The γ-globin variant with higher structure (HBG1) has a significantly larger number of protein interactors (HBG1, average DMS signal of 0.04, 29 interactors; HBG2, average DMS signal of 0.07, 14 interactors; p value = 0.003; KS test; Fig. 2b). While the nucleotide composition of the two transcripts remains nearly the same (HBG1:280c, 463c, 514t, 552a, 575g; HBG2: 280t, 463g, 514g, Δ552a, 574a), the differences between HBG1 and HBG2 are concentrated in regions where the secondary structure is altered (Supplementary Fig. 2). These results indicate that protein interactivity is tightly associated with conformational changes in elements of secondary structure. Interestingly, the increased double-stranded content in HBG1, especially in the 3′-UTR, is accompanied by an accumulation of translation regulatory elements (Fig. 2b) and a concomitant decrease in expression (NCBI Gene ID: 3048).

15We then wondered whether specific RNA structures are involved in protein regulation. We divided the human transcriptome in different classes and analysed their secondary structure as detected by two independent experimental techniques, PARS and DMS. Both techniques show that protein-coding RNAs have the largest structural content (Fig. 2c, Supplementary Table 1). Although part of the mRNA structure is concentrated in the UTRs, when these are excluded, the distribution of the structural content does not change substantially (Pearson’s correlation between transcripts with and without their UTRs = 0.94; Supplementary Fig. 3). The RNAs known to interact with proteins, such as small nuclear RNA (snRNAs) and small nucleolar RNAs (snoRNAs), show the highest amount of structure, whereas RNAs targeting complementary regions in nucleic acids such as antisense, miRNAs and a number of long intergenic non-coding RNAs (lincRNAs), feature the smallest amount of structure (Supplementary Table 1).

16In agreement with our findings, Seemann et al. previously observed a tight relationship between protein binding and conservation of structural elements in mRNAs, which occur to a lesser extent in long non-coding RNAs. Although lincRNAs show a lower amount of double-stranded regions (lowest in PARS, third-lowest in DMS), we note that some of them, such as for instance NEAT1 and XIST, are able to scaffold protein assembly through structured domains. As there is an ongoing debate on the structural differences between coding and non-coding transcripts, and our analysis of DMS and PARS data reveals contradicting results for specific RNA types, we suggest further investigations in future studies (Fig. 2c; Supplementary Table 1).

17To investigate functional differences between highly and poorly structured RNAs, we analysed GO terms associated to the least and most structured RNAs (100 LS vs. 100 HS transcripts) using the cleverGO approach. While the LS set (14 non-coding RNAs and 86 mRNAs) is not associated with specific semantic similarity clusters (total of 36 terms with p value <0.05; Bonferroni test), the HS set (100 mRNAs; total of 395 terms with p value <0.05 and 103 terms with p value <0.01; Bonferroni test; Fig. 2d) includes 20 distinct clusters. The five main categories associated with the clusters and covering at least a quarter of the entries are: (i) complex protein regulation (49/103), (ii) nucleoside metabolic process (39/103), (iii) cellular response (29/103), (iv) gene expression (29/103) and (v) protein targeting (28/103). We also repeated the GO term analysis using as a background the 25% higher expressed transcripts and obtained similar results (K562 strain GENCODE, Methods, Supplementary Fig. 4).

18The cluster analysis reveals the intriguing finding that transcripts with strong structural content interact more with polypeptides and code for proteins involved in regulatory functions and in the formation of complex contact networks. Given the relationship between RNA structure and number of protein interactions (Fig. 1), one preliminary interpretation of our results is that a high degree of control is required for genes that coordinate the activity of a large number of cellular networks. Thus, our analysis suggests a ‘recursive’ property: highly contacted transcripts code for highly contacting proteins (Fig. 2e),.

19To understand the molecular basis of the structure-driven interactivity of RNA molecules, we analysed which physicochemical properties of the proteins better discriminate the HS and LS sets. We studied all 10 variables used in the catRAPID algorithm (Fig. 2f), and removed them one by one to estimate the impact on the prediction of RNA–protein interactions. We found that the capacity to distinguish between the least and most structured RNAs (100 HS and LS transcripts; Supplementary Data 3) sets is more affected when the polarity (p value = 0.28; KS test) and α-helical propensity (p value = 0.06; KS test) are removed (Fig. 2f). The property that more significantly affects the HS binding propensity is polarity, which is enriched in structurally disordered proteins and anti-correlates with hydrophobicity that is key in macro-molecular recognition (Supplementary Table 2). As for the α-helical propensity, we note that the helices are the most frequent structural elements involved in the formation of contacts with double-stranded regions and occur in dsRBD and zinc fingers (Supplementary Table 3). Our observation suggests a possible co-evolution between proteins and RNAs: while the RNA adopts complex shapes to expose binding regions, proteins change their structural content. In agreement with the key lock theory, we propose that natural selection favours highly structured RBPs as interactors of dsRNAs.

20We validated the importance of protein polarity and helical structure by comparing three datasets of well-studied RBPs (human and yeast)– and two sets of proteins retrieved from UniProt (all organisms) as either exclusively ssRNA binders (453 proteins) or dsRNA binders (390 proteins; Supplementary Data 4). Analysis of biophysical properties with the cleverMachine approach revealed that ssRNA binders and dsRNA binders differ for two properties: disorder and α-helix content (Fig. 2g). The comparison of the two sets, one against the other, indicate that RBPs interacting with highly structured RNAs are structured and hydrophobic, while disordered and polar RBPs associate with less structured RNAs (Supplementary Fig. 5). Thus, our analysis further expands what was previously reported for protein–protein interaction networks, in which structural disordered regions have been shown to play a central role, and suggests new rules for nucleotide base pairing with amino acids.

21The analysis of the human transcriptome and across organisms indicate that highly structured RNAs are prone to interact with polypeptides and, in turn, code for proteins involved in biological processes associated with large and complex contact networks. To better investigate the structure-driven protein interactivity of RNA molecules, we focused on a class of transcripts coding for proteins interacting with several partners. The natural choice for this analysis is the molecular chaperones, as they promote folding into the native state and organize the assembly of phase-separate RNP assemblies, thus fulfilling the ‘recursive’ property presented in Fig. 2d. eCLIP data show that most of the RNAs coding for human chaperones are involved in interactions with multiple proteins (Supplementary Fig. 6). We found a significant correlation between protein–RNA and protein–protein interactions annotated in BioGRID (Fig. 3a). This result confirms that transcripts bound by many RBPs also code for highly contacted proteins.

22To understand if the correlation between protein–protein and protein–RNA interactions is a general property or simply a feature of the chaperone family, we analysed interactions of the transcriptome ranked by PARS scores and 24 mRNAs coding for chaperones for which PARS data are available (Genecards; https://www.genecards.org; ‘HSPs’ set; Methods, Fig. 3b). We found a positive correlation between the amount of RNA structure and the number of BioGRID interactors of the encoded proteins (Supplementary Fig. 7a–b). Thus, our calculations agree with the GO analysis (Fig. 2d) and suggest a relationship between mRNA and their coding partners: highly structured RNAs code for highly interacting proteins.

23The data presented so far suggest that RNAs related by type (e.g. miRNA, snRNA) or function (e.g. coding for chaperones) share similar structural characteristics (Fig. 2). Thus, it should be possible to estimate differences in the interaction network of two unrelated transcripts by analysing their structural content, and vice versa. To test this hypothesis, we selected the highly structured HSP70 transcript (HS RNA, log of PARS score of −1.3 corresponding to 26% of double-stranded content, Fig. 3c) coding for a chaperone essential to regulate protein complex assemblies such as clathrin coats and stress granules,. As a control we chose the RNA coding for BRaf that is less structured (LS RNA, score of −2.8 indicating 6% of double-stranded content according to PARS, Fig. 3c–e) and encoding for an oncogene involved in transmission of chemical signals from outside the cell to the nucleus (the structural comparison is confirmed by CROSS predictions and DMS experiments, as shown in Supplementary Fig. 8).

24We found that HSP70 has a larger number of partners (30 RBPs identified by eCLIP) than BRaf (9 eCLIP RBPs, 6 in common to HSP70, Supplementary Fig. 9), which is perfectly in agreement with the structure-driven protein interactivity property. In keeping with the trend of Fig. 1b, catRAPID indicates that proteins have a larger propensity to bind to HSP70 than BRaf (Fig. 3f). Moreover, the highly structured HSP70 codes for a protein with a higher number of interactors (244 BioGRID physical interactors), while the poorly structured BRaf has a protein product binding to a smaller set of molecules (88 BioGRID physical interactors). Our observations suggest that an RNA with a large number of interactions is prone to act as a network regulator: we speculate that, because of the higher interactivity, HSP70 transcript could perform as a chaperone depending on the context.

25Thus, we hypothesize that a structured RNA, because of its higher protein-interacting potential, is able to affect the protein interaction network more than a poorly structured RNA. In a proof-of-concept experiment, we used a chemical compound, biotinylated isoxazole (b-isox) to induce formation of a liquid-to-solid phase transition of a protein assembly, that we incubated with either HS (HSP70) or LS (BRaf) transcripts (Fig. 4a and Supplementary Fig. 10). We observed that HS altered the composition of the protein aggregate more than LS RNA (Fig. 4b and Supplementary Data 5). Indeed, when HS RNA was added, a significant change of concentration was observed for 29 proteins (Fig. 4c; 21 ‘released’ set, black dots, and 8 ‘kept’ set, red dots in Fig. 4b), while only nine proteins were identified in the LS RNA case. Thus, the composition in presence of LS RNA remained similar to that of the background control (‘static’ set, grey dots in Fig. 4b).

26We reasoned that the competition of RNA with the b-isox precipitate contact network, could be the result of either direct or indirect protein–RNA interactions (Fig. 5a). Yet, catRAPID predictions support the hypothesis of a direct effect: an increase in the experimental stringency (Supplementary Fig. 11; Methods) is also associated with an increase in the theoretical predictive power (Fig. 5b). In accordance with our previous analysis of RNA-binding preferences, proteins released upon HSP70 incubation result significantly deprived of polarity (Fig. 5c). Thus, our experiment suggests that the structure-driven protein interactivity of RNA molecules is active at every level, promoting individual interactions and altering the composition of condensates (Fig. 2e).

27Owing to recent advances in high-throughput sequencing, it is now possible to collect information on the majority of RNAs. Large-scale experiments unveiled many functions of transcripts,: chromatin modification, protein assembly and phase separation, among others. At the molecular level, there are still many questions to be addressed in order to understand the full picture.

28Here we focused on the relationship between RNA secondary structure and ability to interact with proteins. It is widely accepted that the structure of a molecule determines all aspects of its life, from stability to function. Yet, to the best of our knowledge, we are the first to report that the structural content of an RNA is intimately connected with the number of protein binders. We demonstrated the solidity of this observation by analysing PARS, DMS, crystals, protein microarrays and CLIP datasets, and also carrying experimental work on RNP aggregation. Our results are not completely unexpected, since lack of RNA structure is linked to more flexible and variable conformations and, thus, a shorter residence of proteins. Moreover, it should be considered that presence of a native fold favours the formation of stable and well-defined binding site that promotes functional roles and, in turn, evolutionary selection. Thus, our findings are reminiscent of the ‘target accessibility’ model that links protein–RNA recognition to the secondary structure of physical contacts,.

29The trends presented in this work suggest the existence of a layer of regulation that directly associates an RNA with its protein product,. In agreement with our GO ontology analysis, it has been previously reported that the transcripts with high conserved secondary structure are enriched in regulatory processes in plants. Indeed, a tight relationship exists between the number of protein contacts of transcripts and the participation of encoded proteins in the network, which reveals an important level of transcriptional regulation for highly connected genes: proteins involved in many interactions are encoded by RNAs bound by several proteins. Our observation indicates an important functional connection, since highly contacted proteins participate in many cellular processes and require tight control at the post-transcriptional level,,.

30By studying transcripts sharing high sequence similarity, such as HBG1 and HBG2, we found that the structure-driven protein interactivity is finely regulated by the exact position of nucleotides. Through the study of a phase-separated RNP assembly we proved that a highly structured RNA such as HSP70 is able to transform the interaction network by competing with pre-existent protein interactions. The main effect observed was the release of proteins from the aggregate; proteins computationally and experimentally tested to be direct interactors of HSP70.

31While the role of HSP70 protein as a protein chaperone is well documented and there are reports on its binding to hydrophobic peptide domains to prevent aggregation and facilitate protein folding, very little is known about the property of the mRNA. Our data suggest an intriguing activity of the HSP70 transcript as a chaperone solubilising the protein assembly, which creates a connection between RNA and protein activities. Our findings are in agreement with previous reports indicating that RNAs are directly involved in RBPs assembly. Indeed, ribosomes have been shown to be powerful co-factors aiding the folding of polypeptide chains as they emerge from their channel. Thus, we speculate that other RNAs can act as chaperones assisting the assembly of proteins. As each transcript is continuously handed off from one protein to another, we expect a mutual chaperoning effect of proteins on RNA and RNA on proteins, which is likely the result of the co-evolution between the two molecules.

32In the future, kinetic analyses tracking the RNA–protein association will be needed to further elucidate to which extent protein partners actively contribute to RNA structure formation. Our findings are reminiscent of the lock-and-key model in the field of enzymology: the structure of both, enzyme and substrate, are key determinants of their association. Yet, structure contributions are not trivial in the case of RNP associations because the combination of different nucleotides bears an obvious specificity-determining potential. While unfolded regions promote protein–protein assembly and disordered proteins exploit short motifs to ensure high connectivity, the reduced nucleotide alphabet and its complementarity suggest that nature favours structure to connect RNAs with proteins.

33The observations presented here, from transcriptome-wide to single molecule analyses, indicate that RNA controls gene regulation at multiple levels. The correlation between RNA structure and number of protein interactions could lead to the discovery of RNA functions that are presently unknown. As demonstrated in the case of protein aggregation, there are RNA-based mechanisms that control phase separation and could be important for the formation of membrane-less organelles. Overall, the complexity and diversity of protein–RNA networks reported here open the avenue for the investigation of regulatory processes.

34To profile the secondary structure of human transcripts, we used PARS data,. PARS distinguishes double- and single-stranded regions using the catalytic activity of two enzymes, RNase V (able to cut double-stranded nucleotides) and S (able to cut single-stranded nucleotides),. Nucleotides with a PARS score higher than 0 indicate double-stranded conformation, while values lower than 0 are considered single stranded,. Transcripts with all nucleotides undetermined were discarded from our analysis. To measure the PARS structural content, for each transcript, we normalized by length and computed the fraction of double-stranded regions over the entire sequence. Given the stepwise function ϑ(x) = 1 for x > 0 and ϑ(x) = 0 otherwise, we computed the fraction of structured domains as:where V(i) and S(i) are the number of double- and single-stranded reads. The top (larger fraction of double-stranded nucleotides) and bottom (larger fraction of single-stranded nucleotides) transcripts are listed in Supplementary Data 1.

35To measure the secondary structure content of the human transcripts without the 5′- and 3′- UTR, we retrieved the corresponding locations of the 5′- and 3′-UTR from Ensembl database (section Database of RNA types) and repeated the same procedure described above eliminating the PARS values corresponding to the UTRs.

36We obtained secondary structure data from in vitro DMS modification published by Rouskin et al. The number of reads of each transcript was normalized to the highest value (as in Rouskin et al. Fig. 2a) and then averaged. Transcripts with available DMS data for <10% of the entire sequence length were excluded from the analysis.

37The full list of protein–RNA interactions calculated with catRAPID is available in the RNAct database. These interaction predictions were calculated over several months on a shared set of 80 HP BL460c nodes with two Intel Xeon E5-2680 2.70 GHz CPUs and 120 GB of usable DDR3-1600 memory each, using eight cores per cluster job, corresponding to 120 years of calculations in a high-throughput cluster. catRAPID omics was used to compute the interaction propensities of 100 LS and HS against the human RNA-binding sub-proteome. catRAPID omics uses large pre-compiled protein libraries and ranks interactions based on the propensity score, with strong predictive power as well as presence of motifs and RBDs. We used catRAPID omiXcore to calculate the interaction propensities of HSP70 mRNA with the ad hoc set of proteins present in the b-isox precipitate.

38RPISeq is a sequence-based predictor of protein–RNA interactions based on two distinct classifiers, SVM and RF. RPISeq is available only through webserver submission and for a limited number of inputs (1 RBP against 100 transcripts per batch run). To circumvent this practical problem, we used CD-HIT (similarity clustering threshold of 85%) to select a representative group of 50 distinct RBPs out of our pool of 579 RBPs. Similarly, we built a set of 100 transcripts (ranked by PARS scores and divided into two groups, as reported in Results, ‘Highly structured RNAs bind a large amount of proteins’: highest values or HS, and lowest values or LS). Using both RPISeq approaches (RPISeq-SVM and RPISeq), we predicted that HS has higher interaction propensities (RPISeq-SVM: HS median: 0.80, LS median: 0.70; RPISeq-RF: HS median 0.71, LS median: 0.54; Supplementary Table 2). We estimated the interaction propensities on the same dataset using catRAPID and observed a similar trend (HS median: 0.02; LS median: −0.17; Supplementary Fig. 1a).

39We used human protein–RNA interactions that were identified through eCLIP experiments in two cell lines, K562 and HepG2, dataset downloaded 5 February 2018. The dataset contains interactions for 118 unique RBPS: 92 proteins in the K562 cell line and 76 in the HepG2 cell line. We used input-normalized eCLIP data (SM input normalization). We processed the eCLIP data using bioconductor package for R (https://www.bioconductor.org/) and mapped the genomic locations to UCSC transcript Identifiers.

40To measure the fraction of protein binders for each transcript, we applied stringent cut-offs [−log10(p value) >5 and −log2(fold_enrichment) >3] as in a previous work. The interactions retrieved after applying the above thresholds represent the top 5% of all interactions in the eCLIP database. The number of protein binders for the transcripts at the top and bottom PARS are listed in Supplementary Table 1.

41iCLIP, PAR-CLIP, and HITS-CLIP interaction data were obtained from the POSTAR database (http://lulab.life.tsinghua.edu.cn/postar/) using POSTAR’s thresholds (p value <0.01 for Piranha, score <0.01 for CIMS, score >0.5 for PARalyzer) and mapped to human GENCODE release 27 transcripts (which directly correspond to Ensembl release 91). POSTAR provides the results of two alternative CLIP-seq data processing methods for each CLIP-seq method (two of Piranha, CIMS, and PARalyzer). We prioritized Piranha since it had been applied to all three experimental methods in POSTAR. In total, we retrieved data for 103 experiments on 85 distinct RBPs (23 studied with HITS-CLIP, 58 with PAR-CLIP and 22 with iCLIP).

42To measure the RNA structure preference of the 85 RBPs, we used DMS data (see Methods: RNA secondary structure measured by DMS) and compared the secondary structural content of the interactome of each RBP with the transcriptome as background. An RBP has (i) high structured RNA preference when the median secondary structure content of its transcriptome is higher than that of the background with p value ≤0.01, (ii) ‘low structured’ RNA preference when the median secondary structure content of its transcriptome is lower than that of the background with p value ≤0.01 and (iii) ‘no preference’ for the rest of the cases.

43We predicted the secondary structure of transcripts using CROSS (Computational Recognition of Secondary Structure. CROSS was developed to perform high-throughput RNA profiling. The algorithm predicts the structural profile (single- and double-stranded state) at single-nucleotide resolution using sequence information only and without sequence length restrictions (scores >0 indicate double-stranded regions). In our analysis, we used the CROSS icSHAPE method that is trained on icSHAPE data.

44To find an alternative link between RNA secondary structure and the number of RBP contacts, we analysed protein-array data,. The GEO entry that we used in our analysis includes the protein interactomes of eight RNAs. For each RNA, we computed its number of interactors after applying a threshold on the signal (affinity over background >5, average of two experiments). Then, we computed the percentage structural content of the list of RNAs using CROSS. We then defined the percentage structural content as the percentage of nucleotides having a positive score (double-stranded tendency).

45From the PDB database we retrieved all structures that contain complexes composed of protein and RNA molecules, excluding those that contain DNA or RNA/DNA hybrid and with X-ray resolution ≤2.0 Å (exact Search Parameter: ‘Chain Type: there is a Protein and a RNA chain but not any DNA or Hybrid and Resolution is between 0.0 and 2.0’). These complexes account for 196 PDB identifiers (Supplementary Table 1). For all protein–RNA pairs reported in those complexes, we consider as internal contacts the total number of times that two phosphate atoms of two distinct nucleotides are in a distance <7 Å and as external contacts the total number of times that a phosphate atom and a carbon A atom are in a distance <7 Å. We normalized those numbers by the length of the nucleotide and protein sequence:

46Using physicochemical properties, the multiclever Machine, approache discriminates ssRNA and dsRNA binders from RBPs. The datasets employed in this analysis comprise 453 ssRNA binders and 390 dsRNA binders, as well as three different RBPs sets. For each feature, enrichment or depletion of a set is indicated with a specific colour: green indicates that ssRNA or dsRNA binders are enriched with respect to the general RBP sets; red means depletion with respect to the general RBP sets; yellow indicates no significant differences between the sets (p values <10–5; Fisher’s exact test).

47For sets and results see Data availability section.

48We used the CD-HIT algorithm (http://weizhongli-lab.org/cd-hit/) with similarity threshold of 85% to cluster RNAs for which DMS (see Methods: RNA secondary structure measured by DMS) and eCLIP data are available. CD-HIT returned 22 clusters (total 55 transcripts) meeting the criterion that a cluster should be populated with at least two transcripts. Pearson’s correlation was estimated for each individual cluster. HBG1 and HBG2 were the two transcripts sharing the highest similarity (99.31%). HBG1 was found to be significantly more structured (DMS score: 0.04) than HBG2 (DMS score: 0.07) and in agreement with our observations HBG1 has more eCLIP interactors than HBG2 (29 and 14, respectively).

49The gene type information for all RNAs that PARS value was available was extracted from the Ensembl database version GRCh38.p12 (https://www.ensembl.org/index.html).

50We retrieved the heat-shock protein members from Genecards database (www.genecards.org) searching for keyword criteria: ‘HSP heat shock’. This returned as a result 1258 entries. Out of these results, we looked for the term ‘Heat Shock Protein’ in the description and removed entries referring to ‘DNA Heat Shock Proteins’, entries described as ‘Pseudogenes’ and ‘Heat Shock Transcription Factors’. After the above criteria, the dataset was reduced to 27 gene names, and for each of those we retrieved from Ensemble, the corresponding UCSC identifier. As described (see Methods, ENCODE database of protein–RNA interactions (eCLIP)), our compiled database for eCLIP is annotated with UCSC identifiers. We used the UCSC identifiers of the HSP family to retrieve the fraction of Proteins that they bind.

51We used the cleverGO algorithm to analyze GO terms. The algorithm provides dynamically organized visualization of the GO terms and grouping depending on the strength of their internal connections. Separate analyses are generated for biological process, molecular function and cellular component ontologies.

52Parameter of the GO semantic grouping:

53Similar strength threshold: 0.5

54Minimal precision: 0.6

55Minimal level: 0.3

56P value cut-off: 0.01

57For sets and results see Data availability section.

58Protein–protein interaction information was retrieved from BioGRID, version 3.4.163, which contains 312,474 human non-redundant interactions (https://thebiogrid.org/). For each protein of interest we counted the number of unique interactors (physical) as defined by their Gene Name.

59To assess whether ssRNAs exhibit different trends from dsRNAs, we used the Wilcoxon’s test (Mann–Whitney U test). Wilcoxon’s is a non-parametric test used to compare the mean of two distributions without any given assumption about them. To compare properties and measure the difference between RNA sets, we used the KS test. KS test is also a non-parametric test used to compare the distance between two cumulative distribution functions.

60For in vitro transcription, linearized plasmids were purified by phenol/chloroform/isoamylalcohol extraction and ethanol precipitation. Ampliscribe T7 high Yield kit (Epicenter, Madison, WI) was used to transcribe 1 μg of template for 4 h at 37°C. The reaction mix was digested with DNaseI and RNA was purified using MegaClear transcription clean-up kit (Thermo Fisher Scientific). About 30 million HeLa cells were then trypsinized and collected. They were washed two times with ice-cold phosphate-buffered saline, resuspended in 300 µL lysis buffer (20 mM Tris, 150 mM NaCl, 0.5% NP-40, 0.1 mM PMSF, 10% glycerol, phosphatase inhibitor, RNasin (Promega) 1:100) and incubated at 4 °C for 20 min with gentle shaking. The lysate was centrifuged at 21,000 × g, at 4 °C for 15 min and the supernatant was collected. Protein concentration was measured using Bradford reagent and BSA standards and normalized to 2 µg/µL with lysis buffer. B-isox (10 mM stock in dimethyl sulfoxide) was added to the lysate at 100 µM and the mixture was incubated at 4 °C for 1 h with gentle shaking. After 1 h, the lysate with b-isox was splitted into three tubes (3 × 120 µL) and 1 µg HSP70 mRNA, 1 µg BRaf mRNA or the same volume of water and incubated again at 4 °C for 1 h with gentle shaking. The aggregates were then collected by centrifugation at 10,000 × g at 4 °C for 10 min. The supernatant was collected and 2× sample buffer was added. The pellets were washed two times by being resuspended in 50 µL lysis buffer, vortexed, kept on ice for 10 min and centrifuged again at 10,000 × g, at 4 °C for 10 min. The washed pellets were resuspended in 50 µL MS buffer (4% SDS, 100 mM HEPES, pH 7.6, 150 mM NaCl) and frozen before mass spectrometry analysis (or in 50 µL 2× sample buffer for sodium dodecyl sulfate-polyacrylamide gel electrophoresis and Coomassie blue staining). Subsequently, the samples were processed according to the FASP protocol using 30k filtration units (MRCF0R030, Millipore) for data-dependent LC-MS/MS analysis with a Q Exactive Plus mass spectrometer.The MS data were analyzed using the software environment MaxQuant version 1.5.3.30. Proteins were identified by searching MS and MS/MS data against the human complete proteome sequences from UniProtKB, version of November 2015, containing 70075 sequences. Carbamido-methylation of cysteines was set as fixed modification. N-terminal acetylation and oxidation of methionines were set as variable modifications. Up to two missed cleavages were allowed. The initial allowed mass deviation of the precursor ion was up to 4.5 ppm and for the fragment masses it was up to 20 ppm. The ‘match between runs’ option was enabled to match identifications across samples within a time window of 2 min of the aligned retention times. The maximum false peptide and protein discovery rate was set to 0.01. Protein matching to the reverse database or identified only with modified peptides were filtered out. Relative protein quantitation was performed using the LFQ algorithm of the Maxquant with a minimum ratio count of 1. Bioinformatic data analysis was performed using Perseus (version 1.5.2.6). The proteins with minimum three valid LFQ values in at least one group (background/Hsp70 mRNA/B-Raf mRNA) of 4 biological replicates were considered as quantified and used for downstream analysis. Proteins significantly changed in the presence of Hsp70 or BRaf mRNA compared to the background control (water) were identified by two-sample t-test at a p cutoff of 0.05.

61catRAPID omiXcore integrates the interaction propensities of individual protein and RNA fragments into a unique score that was used to compare predictions with experiments. To this aim, we computed the AUC ROC of two protein groups: (i) the positive set was comprised of all proteins that were significantly released in the presence of HS mRNA (HSP70; −log(p value) >1.3 and log(RNA/Background) <0) and not significantly released in the presence of LS mRNA (BRaf; −log(p value) <1.3 or log(RNA/Background) >0); (ii) as negative set we used all the proteins that did not show any significant change in the presence of either HS or LS mRNAs (HSP70 or BRaf, respectively; −log(p value) <1.3) in the b-isox precipitate.

62We used a cut-off to define the experimental stringency: (i) highly stringent cases, all those proteins with minimum fold change in the presence of both the two mRNAs (−log(p value) <1.3 and |log(RNA/Background)| <0.15), (ii) medium stringent cases, those that show modest change (−log(p value) <1.3 and |log(RNA/Background)| <0.30) and (iii) low stringent cases, those that show greater change (−log(p value) <1.3 and |log(RNA/Background)| <1).

63We built the ROC curves and computed the corresponding AUCs for the three different cases (low stringent, medium stringent, high stringent) and we observed higher performances when increasing the experimental stringency.

64Further information on research design is available in the Nature Research Reporting Summary linked to this article.