Zooming in on protein–RNA interactions: a multi-level workflow to identify interaction partners

Since their discovery and until recently, RNA-binding proteins (RBPs) have been identified by the presence of one or more RNA-binding domains in their sequences [1]. However, concomitantly to a new appreciation for RNA as key biological macromolecule acting at post-transcriptional level [2–4], there has also been a re-evaluation of what constitutes an RBP. Since one of the principal ways by which RNA exerts its function is by the formation of ribonucleoprotein complexes, every protein capable of establishing even weak and extemporary interactions with an RNA molecule may be defined as RBP [5,6]. The interactions of proteins with RNA can be highly dynamic and heavily dependent on the cellular environment [7], which makes the goal of defining the range of affinities and specificities quite challenging [8]. In fact, indiscriminate binding of RNA by RBPs is a quite common phenomenon [9], and the assumption that stronger affinity translates into more relevant biological functions is not necessarily correct [10]. For the scientific community, the revelation of the dynamicity and malleability of the partnership between RBPs and RNA allows for the exploration of new possible interaction mechanisms, networks, genes and protein regulation systems to investigate. It becomes, therefore, increasingly important to complete the catalogue of eukaryotic RBPs, at present containing >315 000 elements (about 6000 orthologues from >150 species) of which 3500 are from H. sapiens and can be divided in conventionally-defined RBPs and several classes of non-canonical RBPs [5,11,12].

Large-scale identification of new potential RBPs can reveal unexpected biological and pathological functions and, when confronted with a novel RBP, the identification of its RNA binding partners is a critical step to define the protein's cellular and molecular roles. To achieve this goal, increasingly sophisticated high-throughput methodologies have been developed, spanning from methods that aim to preserve the native cellular RNA–RBP interactions [13–15] to finely-controlled in vitro techniques that allow to define kinetics and dynamic parameters for each binding pair [16,17]. However, the field that has probably seen the biggest evolution in the shortest time span is the one of computational prediction algorithms [18–22]. The wide range of computational tools available includes several data-driven methods based on learning models, in which the algorithms are trained using experimental outcomes and databases to identify RNA–RBP binding patterns and define the genome-wide profiling of RNA–protein interactions [22].

In this short review, we propose a work pipeline for the identification of the RNA binding partners for novel putative RBPs. Starting from in-cell data harvesting, we would like to guide the reader through the employment of the different tools offered by catRAPID [23,24], our in-house developed RNA–protein interaction prediction algorithm, and propose some indications on how to validate the outcome experimentally. Our wider goal is to support the scientific community in the identification of novel biologically relevant non-canonical RBPs.

The validation and training of computational algorithms for the prediction of protein-RNA interactions is strongly supported by experimental data (Figure 1). Prediction software can be significantly enriched by the output of techniques able to identify a protein's RNA partners within the cellular environment; such procedures are key to defining native interaction pairs and to monitoring responses and variations upon physiological stimuli or under pathological stress.

The main tool to obtain information about the RNA binding partners of a target protein in the cellular environment is immunoprecipitation (IP), a widespread technique to pull down the protein of interest together with its physiological RNA binding partners. RNA immunoprecipitation (RIP) requires incubation of cell lysates with an antibody raised against the target protein [25]. RNA molecules bound to the target protein can then be isolated and analysed to reconstruct physiological native complexes formed within the cell. RIP can be coupled to either microarrays (RIP-Chip) or high-throughput sequencing (RIP-Seq). In either case, RIP can only determine the identity of the RNA molecules associated to the target protein, unless digestion-optimized RIP (DO-RIP) is performed [26]. This variation of RIP introduces an RNase digestion step to preserve only the portion of RNA bound to the protein, allowing for binding-site mapping [27]. If, instead, the interest is focused on identifying multi-subunit ribonucleoproteins, the most appropriate RIP variant may be RIP in tandem (RIPiT), which employs two distinct IP steps performed either with antibodies against different proteins of the complex or with antibodies binding different regions of the same target protein [28]. Information about native protein–RNA complexes formed within the cell can also be obtained by employing affinity tags [29], without having to rely on the antibody's specificity and sensitivity.

To overcome RIP's limitations (enrichment of indirectly bound RNAs, detection of interactions not present in cell but formed after lysis, loss of weaker interactions due to the required stringent washing conditions), cross-linked RNA immunoprecipitation (CLIP) has been developed [30]. CLIP promotes the stabilisation of the bonds between a protein and its interacting RNA, generally by UV radiation [31]. A large amount of CLIP variants is available, and most of them can offer high-resolution results at the single nucleotide level. Pioneers among these are, for example, the high-throughput sequencing CLIP (HITS-CLIP), that enriches the RNA population for sequences corresponding to the RBP binding sites [32]; the individual nucleotide-resolution CLIP (iCLIP) [33,34] and enhanced CLIP (eCLIP) [35], that utilise different oligonucleotide adapter configurations to obtain RNAs of different lengths employed to build the interacting fragments at single-nucleotide resolution; and the photoactivatable ribonucleoside CLIP (PAR-CLIP), that relies on metabolic incorporation of labelled ribonucleoside analogues that yield photo-adducts when cross-linked at selected wavelengths [36]. Radiation-free CLIP variants, such as infrared-CLIP (irCLIP), have also been reported [37]. To overcome the restrictions imposed by the use of irreversible UV cross-linking, formaldehyde-based techniques such as fCLIP, which are more efficient in capturing interactions with dsRNA [38], can also be employed [39]. One limit of CLIP-Seq approaches is that low crosslinking efficiency makes low abundant transcripts difficult to detect, while transcripts present at high amounts are usually over-represented in IP samples. This issue can be partial solved with a proper bioinformatic analysis.

High-throughput sequencing of the RNAs isolated through RIP, CLIP and related protocols yields millions of short ‘reads’, which represent the sequenced portions of cDNA fragments obtained through RNase titration (a step omitted in standard RIP-Seq), followed by reverse transcription and PCR amplification. The RNase treatment allows to obtain fragments long enough to be uniquely mappable but short enough to identify the binding site with the highest possible resolution. Library preparation ends with the production of cDNA fragments flanked by adapters that allow amplification and sequencing, resulting in the generation of short reads that can undergo bioinformatic analysis aimed at identifying RNA targets for the RBP (Figure 2). Reads pre-processing steps, including demultiplexing and adapter trimming, are often required, especially in those cases in which Unique Molecule Identifiers (UMIs) are used [40]. UMIs are random barcodes which identify unique cDNA fragments, allowing to detect and remove PCR duplicates that are commonly produced during CLIP-Seq library preparation [34,35].

To identify the RNA molecules from which they derive, reads are aligned to a reference genome using splice-aware alignment programs commonly used to analyze standard RNA-Seq data, like TopHat2 [41] or Star [42]. If the RBP under investigation binds mature mRNAs, reads can be mapped directly to the transcriptome [43] using a splice-unaware mapper like Bowtie2 [44]. Reads coming from HITS-CLIP experiments have high mutation rates (usually deletions) at the cross-linking site (CIMS, standing for cross-linking induced mutation sites), which are due to residual amino acids hindering the reverse transcriptase [45–47]. Similarly, the use of 4-thiouridine (4-SU) or 6-thioguanosine (6-SG) in the PAR-CLIP protocol leads to a high number of transition events (T to C or G to A, respectively) at the cross-linking sites [48]. Reads mapping can be improved by taking into account the high rate of such mutations. For instance, the splice-unaware BWA aligner [49] has been modified in order to incorporate an error model that favours PAR-CLIP specific transitions [50,51]. Reads produced by iCLIP, and likely eCLIP, experiments do not require such special treatment, since such protocols enrich for cDNAs truncated at the cross-linking site (CITS, standing for cross-linking induced truncation sites), while only a minor proportion of fragments represent CIMS-containing read-through cDNAs [43].

Post-processing of aligned reads is a mandatory step in RIP-Seq and CLIP-Seq data analysis. Reads aligning to multiple genomic positions are usually removed [52]. However, such multi-mapped reads can be used to identify regulatory RNA sites localized within repetitive regions [53]. To filter out PCR duplicates, reads mapping to the same genomic position are collapsed; UMIs, when present, can be used to avoid removing natural read duplicates, common in case of high sequencing depth [54].

In RIP-Seq experiments, target RNAs can be identified either by transcript enrichment analysis of IP versus control samples [55–57], which can be performed using procedures commonly adopted in differential expression analysis, or with ad-hoc peak-calling tools [58,59]. Binding site identification in HITS-CLIP experiments is usually accomplished by means of peak-calling approaches, that identify regions which are enriched in reads with respect to their genomic context (gene, transcript, metagene region) [59], the background represented by control experiments (input, IgG, mock IP) [60], or baseline expression profiles [59]. Single-nucleotide resolution in HITS-CLIP data analysis can be achieved by looking for statistically significant CIMS [59]. Such resolution in binding site identification can also be attained in PAR-CLIP, by searching for transition events that are not likely to be caused by sequencing noise, single nucleotide polymorphisms or contamination [61]. The first nucleotide of reads from iCLIP and eCLIP experiments most often marks the truncation site and it is therefore located one nucleotide after the cross-linking site. Statistically significant CITS can be identified by using tools like iCount [62] and PureCLIP [63].

The outcome of the data analysis obtained by means of CLIP-Seq experiments is of enormous value in approaching the study of a potential RBP with unknown RNA partners. However, there are two major drawbacks in limiting the investigation to CLIP-Seq approaches: the protein regions in direct contact with the target RNAs remain unknown, and there are not sufficient data to speculate on the strength of the interaction for each protein–RNA pair. The use of predictive algorithms such as catRAPID would integrate the results of a CLIP-Seq experiment with this information [24] (Figure 1).

catRAPID is an algorithm able to compute protein–RNA interaction propensities with strong predictive power (area under the Receiver Operating Characteristic curve of 0.78 on >1 000 000 interactions) [64], through the calculation of secondary structures, hydrogen bonding and van der Waals contributions. The algorithm was trained on PDB crystals [24] and was later adapted to predict CLIP-Seq interactions with long non-coding RNAs [65]. For large proteins (>750 amino acids) and RNAs (>1000 nucleotides), the algorithm fragments sequences into overlapping segments and computes the interaction propensity through the analysis of physical–chemical properties and secondary structures of the molecules. According to the chosen implementation (Table 1), the method can either reconstruct the overall interaction propensity score for each protein–RNA pair [66] or rank the fragments according to the predicted interaction strength [67]. The outcome of this analysis allows to map both protein and RNA binding sites and to estimate the overall strength of the interactions [68], overcoming the limitations of CLIP-Seq techniques mentioned above. If used in combination with CLIP-Seq data relative to the protein of interest, catRAPID can be employed to select the best targets based on calculated binding strength, but it could be also useful in predicting putative targets that are not expressed within the cellular system. If no information about the RNA targets is available, catRAPID can represent a promising tool for the investigation of the protein's genome-wide RNA-binding potential against an RNA sequence library. Here we suggest a pipeline that could be followed in such circumstance:

1. If the RNA binding potential of a given protein is unknown, catRAPID signature [71] can be used to predict it, along with the putative RNA-binding regions. This approach is particularly recommended if a potential RBP needs to be selected from a panel of candidate proteins;

2. catRAPID omics [69] can then be used to predict the interactions between the protein and a precompiled or custom RNA library. The result is a ranked list of protein–RNA pairs;

3. If the protein of interest is human, its co-expression with the putative interactors in different tissues can be evaluated using catRAPID express [70];

4. Once the most promising binding partners have been identified, catRAPID strength can be employed to evaluate the strength of each interaction [67];

5. Finally, catRAPID fragments [66] can be run on the highest scoring protein–RNA pairs to predict the binding sites. Interactions with long RNAs can be analyzed using Global Score [65] or omiXcore [68].

By narrowing down the number of potential targets and by suggesting the most likely binding sites, such approach can be employed to guide further experimental and computational analyses. Being a predictive tool, there is always the chance that catRAPID may fail in identifying valid RNA targets. Prediction accuracy depends on the set provided to train the algorithm. As more and more data become available, retraining of the algorithm will be necessary to achieve better performances.

catRAPID is only one of several possible computational methods developed for predicting protein–RNA interactions. We would like to mention here some valid alternatives:

• LncPro [72]: similar to catRAPID in the employment of RNA secondary structure, hydrogen-bonding and van der Waals interaction propensities. It is designed to predict whether a specific long non-coding RNA interacts with one or more protein sequences. Propensities are calculated for both protein and RNA and a probability score ranging from 0 to 100 is generated;

• RPISeq [73]: protein and RNA sequences are encoded into features that are then used to train Support Vector Machine (SVM) and Random Forests classifiers;

• RPI-Pred [74]: it combines RNA and protein sequences with predicted or actual 3D structures. The features are then used to train an SVM classifier;

• iDeepS [75]: a deep learning-based method that exploits convolutional neural networks (CNNs) trained on RNA sequences and predicted secondary structures. At the end of the pipeline, a classification layer is responsible for RBP binding sites prediction. Deep learning models are generated individually for each RBP based on available CLIP-Seq data, allowing the formulation of predictions on a limited set of proteins.

Once the RNA targets of a protein have been determined, further analyses are necessary in order to verify the reliability of the results and to gain insights into the biological function of the RBP. Both tasks can be approached by looking at the function of target RNAs. A common way to do that is to perform an Over Representation Analysis (ORA) [76], which consists of identifying over-represented functional categories in a list of genes. Enrichment is evaluated against a background composed of all the expressed genes. A more sophisticated approach, called Gene Set Enrichment Analysis (GSEA) [76], involves ranking genes based on a certain score and evaluating if some categories are enriched at the top or the bottom of the ranked list. An implementation of this method that is specific for CLIP-Seq data is provided by the Seten tool [77]. This program requires as input a set of CLIP-Seq peaks, each with a score assigned by the peak-caller, but it could also work starting from predicted binding sites, as long as an interaction score is provided.

Another common analysis consists in scanning the identified binding sites in order to detect common patterns highlighting RBP binding preferences (motif analysis). Sequence motifs recurring in large sets of binding sites can be discovered using different tools, like MEME-ChIP [78] and SeAMotE [79]. Both tools start from a set of sub-sequences identified in the positive sequence set and evaluate their enrichment with respect to a control sequence set (unbound RNAs). A more recent tool, named mCross, exploits the single-nucleotide resolution offered by CLIP-Seq techniques to enhance the accuracy of de novo motif discovery [80].

Sequence alone may not be sufficient to fully explain the binding specificity of an RBP: a sequence motif could be accessible only when put in a proper secondary structure context. Tools like GraphProt [81], ssHMM [82] and BEAM [83] are able to detect motifs encoding both sequence and secondary structure information.

To validate the prediction accuracy of the computational analysis proposed so far, it is ideal to evaluate each interacting pair within a controlled environment. A most accurate validation should start from the screening of potential binders, followed by the precise determination of binding sites and kinetic and thermodynamic parameters, and completed with structural insights into the drivers of the interaction (Figure 1). A comprehensive review of the methodology is beyond the scope of this article and for more details we refer to a number of recent reviews of the field [16,17,84–86].

As for other molecular interactions, the binding between a protein and an RNA molecule is kinetically characterised by the rate at which they associate (k_on) and dissociate (k_off). Conventionally, the dissociation constant (K_d), which is the ratio between k_off and k_on at the chemical equilibrium, is used to express the binding affinity: the lower the K_d, the greater the affinity. It is however important to note that binding pairs with the same K_d may have different k_on and k_off and therefore different binding mechanisms. To add another layer of complexity, many RBPs have multiple binding regions that may vastly differ in their affinity towards the same RNA [87].

There are several established methods for determining the kinetic and thermodynamic parameters of binding (Table 2). Techniques such as electrophoretic mobility shift assay (EMSA) and filter binding assay can be useful to estimate binding affinities with basic molecular biology tools [88,89]. Both these methods represent a viable first step analysis, especially because of short protocols and limited amounts of samples required. The latter criterion can be crucial in protein–RNA interaction studies, since some RBPs can be very difficult to isolate, tagged RNA synthesis can be expensive and advanced methods for more reliable and accurate determination of molecular binding characteristics generally have specific sample requirements and higher operation costs. Despite these advantages, since EMSA and filter binding experiments are performed within conditions very distinct from the in vivo ones (polyacrylamide gel and nitrocellulose filter, respectively), more reliable kinetic data may be obtained by other techniques. Examples of such methods are bio-layer interferometry (BLI), multi-channel surface plasmon resonance (SPR), microscale thermophoresis (MST), the employment of fluorescence resonance energy transfer (FRET), and the most recent switchSENSE (Table 2, [86–89, 95]). Some of these require the labelling of one of the interactors with fluorescent dyes, as in the case of MST and FRET, or sample immobilisation on biosensors, as needed for SPR and BLI experiments. These requisites can impose certain structural restraints and thus compromise the binding and the results obtained. Comparing how the same RBP binds to RNA molecules that differ only in one nucleotide can help the identification of the RNA portion physically interacting with the protein. However, the determination of the exact nucleotides involved in the binding can be highly challenging without structural studies. Once the RNA targets have been precisely characterised, selective mutations on the protein region thought to be responsible for the binding may be useful to determine which amino acids are directly responsible for the interaction. The same approach could be extended to mutating selected nucleotides or regions of the RNA molecules for further characterisation of the binding spots without obtaining an atomic resolution structure.

Kinetic studies are essential for determining the binding propensity and validating different binding partners [10]. However, experiments conducted in vitro with isolated components often do not allow studies under physiologically relevant conditions. Isolated proteins can be sensitive to higher temperatures, especially under prolonged experiments that far exceed their in-cell lifespans. Since chemical equilibriums are temperature-dependent, in vivo kinetics may thus significantly differ from those measured in vitro. The thermodynamics of binding can also reveal the energetic landscape of protein–RNA interactions [100–102]. It is therefore very important to measure the thermodynamics parameters directly, when experimentally possible, or to calculate them. Isothermal titration calorimetry (ITC) allows direct measurements of the equilibrium constant, stoichiometry and reaction enthalpy (ΔH) at a given temperature [99]. Using the van't Hoff equation, these can then be used to determine the temperature dependency of the equilibrium constant. SPR and MST also allow calculation of thermodynamic data based on the stable temperature of the measurement cell, while BLI is considered less reliable.

Kinetic and thermodynamic data obtained from interaction studies provide a good numerical description of the binding and, through the use of RNA-centric methods, a library of RNA sequences with high binding propensities [103]. However, protein and RNA sequences alone may not be sufficient to fully characterise their binding specificity. This is especially true in the case of non-canonical RBPs that do not contain a consensus RNA-recognition sequence [104]. RNA structure has been shown to be the driver behind most non-canonical RBP binding events, with highly structured RNAs having bigger protein interactomes [105,106]. The approaches to determine the structure of separate components have been extensively reviewed [107,108]. Among the most pertinent techniques to determine the macromolecular structure of a complex, there are X-ray crystallography, nuclear magnetic resonance (NMR) and cryogenic electron microscopy (cryo-EM) (Table 3). However, the definition of structural details by these methodologies can be challenging [107]. A large number of RBPs contain intrinsically disordered regions and tend to form macromolecular condensates, making their crystallisation impractical; NMR, that bypasses the need for crystals, is dependent on the molecular weight of the complexes that can make relaxation times slow and signal-to-noise ratio poor; obtaining results at atomic resolution with cryo-EM remains difficult [108]. These limitations make a strong case for the employment of complementary structural biology techniques [107]. A particularly promising development is the advance of small angle scattering and computational methods for data analysis, in particular small angle neutron scattering (SANS) [109]. Selective deuteration enables a higher degree of contrast between binding partners and can therefore provide a rough position for each component of the complex [110]. Data obtained from the above mentioned methods can be used as structural restraints for molecular dynamics simulations and data-driven docking [111,112] and can be integrated together in a hybrid multi-level approach for studying RNA-protein complexes, thus completing the full circle of integrated methodologies [75,113,119,120].

• Our understanding of many physiological and pathological phenomena cannot be exempted from an in-depth knowledge of protein–RNA interactions underlying them.

• Such comprehension, which goes from the identification of targets RNAs to binding modes characterization, requires a multidisciplinary approach involving biochemistry, molecular biology, bioinformatics and physics.

• As new and more powerful high-throughput methods and predictors are being developed, an integrated and productive usage of both approaches becomes more and more feasible. For instance, predictive tools such as catRAPID, which is general enough to be applied to any protein–RNA pair, could also be employed to improve the specificity of omics studies.

The authors declare that there are no competing interests associated with the manuscript.

The research leading to this work has been supported by European Research Council (RIBOMYLOME_309545 and ASTRA_855923), by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754490 and under projects IASIS_727658 and INFORE_825080, by the Spanish Ministry of Economy and Competitiveness BFU2017-86970-P as well as the collaboration with Peter St. George-Hyslop financed by the Wellcome Trust.

A.C., J.R. and A.V. drafted the sections of the review pertinent to their expertise. E.Z. and G.G.T. guided, supervised and coordinated the work. All authors contributed to writing and editing the manuscript and to the original production of all figures and tables provided.

The authors would like to thank all members of Tartaglia lab.

BLI

biolayer interferometry

CIMS

crosslinking-induced mutation sites

CITS

crosslinking-induced truncation sites

CLIP

cross-linking followed by immunoprecipitation

CLIP-Seq

CLIP followed by followed by high-throughput sequencing

eCLIP

enhanced CLIP

FRET

Förster resonance energy transfer

HITS-CLIP

high-throughput sequencing CLIP

iCLIP

individual nucleotide-resolution CLIP

IP

immunoprecipitation

MST

microscale thermophoresis

NMR

nuclear magnetic resonance

PAR-CLIP

photoactivatable ribonucleoside CLIP

RBP

RNA-binding protein

RIP

RNA immunoprecipitation

RIP-Seq

RIP followed by high-throughput sequencing

SPR

surface plasmon resonance

SVM

support vector machine

UMI

unique molecule identifier