On the limits of inferring biophysical parameters of RBP-RNA interactions from in vitro RNA Bind’n Seq data

1. Introduction

RNA-binding proteins (RBPs) interact with RNAs at every step of their life cycle. Their modular structure, usually an assortment of RNA-binding domains, underlies their ability to interact with both RNAs and proteins, and couple various layers of gene expression.¹ While $\sim 2500$ RBPs are currently known, most remain to be functionally characterized. A first step in this process is to determine the interaction partners and the sequence/structure specificity of the RBP. Many RBPs recognize their targets in a sequence-specific manner, although the accessibility of binding sites within the targets also plays a role.² The sequence specificity is usually represented in a position weight matrix (PWM), which specifies the probability of finding each of the four nucleotides at each position in the RBP binding site. This is an obvious simplification, as dependencies between positions in the binding site likely occur. However, training more complex models requires substantially more data, which are often not available. Moreover, the improvement in binding site predictability by more complex models is modest, at least in the case of other nucleic acid binding proteins, transcription factors.³ With the realization that the presence of a canonical RNA-binding domains is not necessary for the ability of a protein to bind RNAs⁴ came a pressing need to determine the determinants of RNA-RBPs interactions and the sequence/structure specificity of the proteins newly found to interact with RNAs.

The past two decades have seen the development and broad application of experimental methods for RBP target identification. They include in vivo high-throughput approaches such as HITS-CLIP, PAR-CLIP, iCLIP and eCLIP (reviewed in Ref. 5), and more recently-developed in vitro approaches such as RNA Bind’n Seq.⁶ While the CLIP methods rely on the sequencing of RNAs that interact and can therefore be crosslinked to RBPs in vivo, RNA Bind’n Seq relies on the affinity-dependent interaction of RBPs with random RNAs in vitro. The oligonucleotides whose interaction with the RBP (or domain thereof) of interest are computationally analyzed to identify short sequence motifs that are enriched in the affinity-selected pool of RNAs. So far, analyses of such data involved the identification of enriched k-mers (short oligonucleotide sequences of a specified length, $k$), and then a greedy alignment procedure yielded PWM representations of the RBP binding motifs. This left open the question of whether the derived PWMs accurately predicted the interaction energies of RBPs with their binding sites. In contrast, the aim of our work was to develop a biophysics-anchored method to directly infer the PWMs from RNA Bind’n Seq data. Our paper is organized as follows: Sec. 2 explains how we derive our thermodynamical model. We comment on the practical implementation of this model in Sec. 3, where we also explain how we account for sequence composition biases in the pool of oligomers. Results for different RBPs are presented in Sec. 4, where we also comment on the accuracy of the results obtained from this type of data for different RBPs. Concluding remarks are given in Sec. 5, and a list of publicly available data sets that we analyzed is provided in Sec. 6.

2. Model

Our model is an adaptation of a Bayesian, thermodynamic model that was constructed to infer di-nucleotide weight tensors from SELEX data.⁷ In the following, we derive the log-likelihood of Bind’N Seq data given the PWM for the RBP of interest, which will be inferred by expectation-maximization as described in Sec. 3.

We assume that an RBP binds an oligomer over a binding site $s$ of length $L_w$ and that the likelihood of the binding taking place, according to Boltzmann’s law, goes as $\propto exp\left({\sum }_{i=1}^{L_w}{E}_i^{s_i}\right)\equiv {\mathrm{e}}^{E\left(s\right)}$, where $s_i$ is the nucleotide at position $i$, so $s_i\in \left\{\mathrm{A},\mathrm{C},\mathrm{G},\mathrm{T}\right\}$. Therefore, each element of the position weight matrix (PWM) $M$ can be identified with $m_i^{\alpha }\equiv exp\left(E_i^{\alpha }\right)$, with their columns being normalized as ${\sum }_{\alpha }{m}_i^{\alpha }=1\forall i=1,\dots ,L_w$.

Additionally, we account for the fact that there are genuinely two different ways of binding, sequence-specific binding as described by the PWM, and unspecific binding to RNAs with a probability $exp\left(E_0\right)$. Combining these two possibilities, we arrive at the probability for an RBP binding to a site $s$

Consequently, the chance of an RBP being bound somewhere on a longer oligomer $S$ with $L_S\ge L_w$ is

Eventually, the logarithmic likelihood of our library of oligomers $D$ reads

3. Implementation

Our goal is to optimize the parameters in (2.5) such that they maximize the likelihood of our library to be realized in the present way. As a side note, it is equivalent to optimize $P\left(D\right)$, or its logarithm because the logarithm strictly increases as its argument increases and decreases as its argument decreases. Since the libraries are typically quite big it is beneficial for us to maximize the logarithm in order to keep the overall numbers under control. While the library $D$, the copy number of a read $n_S$, the read and binding site length $L_S$ and $L_w$, and – with some limitations – the frequency priors $f_S$ are given from our data, the position-specific binding encoded by the PWM and the position-unspecific binding ${\mathrm{e}}^{E_0}$ have to be obtained during the optimization process. Eventually, we want to obtain the the PWM, whereas ${\mathrm{e}}^{E_0}$ represents a hidden parameter which will be inferred via the expectation-maximization procedure. In principle, this would also apply to the concentration $c$ but none of our final expressions depend on $c$ any more due to the linearization. Before diving into the details of the EM procedure’s implementation we would like to comment on how to infer the frequency priors $f_S$.

4. Results

In analyzing Bind’n Seq datasets for various RBPs, we found that only a small subset of random initializations deliver a convergent EM process. The larger $L_w$ is the larger is the space of possible initializations, therefore it becomes increasingly unlikely to accidentally hit a region of initialization which converges. This could be compensated for by increasing the number of runs by a factor of $4$ for each additional position in the PWM. To avoid that, one can use the knowledge of previous runs, done with shorter PWMs and initialize the longer PWM from the shorter PWM, filling up the additional entries with randomly initialized values to check if the shorter PWM is part of a longer PWM. We carry on with this procedure until the EM algorithm does not find any minimum any more amongst 12 different random initializations of $E_0$ and the random initialization of the PWM at the additional positions. Sometimes, neither a fully random intialization nor an initialization of the PWM “guided” by prior knowledge lead the EM algorithm converging to a minimum in log-likelihood (2.5). The relative efficiency of the algorithm finding true maxima is displayed in Figure 1 for the RBPs discussed in the following section. We consider an outcome of the algorithm to be a “true” maximum if the posterior log-likelihood is larger than the initial one and the algorithm is not stuck in a region where $E_0$ is large, meaning that the unspecific binding dominates. The maximization algorithm is eventually terminated by a limit of 200 iterations. For readability, we list the used Bind’n Seq data files from Ref. 8 in Table 1 in Sec. 6.

Figure 1. Summary of fraction of convergent outcomes for different investigated RBPs and binding site length $L_{\mathrm{w}}$.

RBPs and binding site lengths with no convergent maxima of the log-likelihood (as described in the corresponding subsections, e.g. subsection 4.6) were discarded. While $E_0$ is always initialized with a negative random number, the PWM can be initialized either “guided” by already obtained shorter motifs or literature motifs (c.f. corresponding subsection of Sec. 4), or with a fully “random” PWM.

4.1 Benchmark: PWM of length 6 for RBFOX2

To benchmark our method, we started our evaluation with RBFOX2, a key regulator of alternative splicing¹³ that was extensively studied with a variety of methods (e.g. Ref. 14). The RBFOX2 Bind’n Seq dataset⁸ consists in nine libraries at nine different protein concentrations and two protein-free control libraries, all containing reads of 50 nucleotides (nts) in length, including the adaptor. RBFOX2 is widely used to benchmark computational analysis methods (c.f. Ref. 15) and thus the corresponding dataset was carefully generated, to include multiple, high-quality libraries. Established techniques like kmer-enrichment analysis and the streaming-kmer-algorithm (SKA) predict a consensus 6mer TGCATG as the most prominent motif followed by other GCATG-containing 6mers.¹⁵ Our results in Figure 2 reproduce the predicted TGCATG 6mer as a part of the motif Figure 2(a). Moreover, we find the subdominant PWM Figure 2(b) which has a quite substantial overlap with Figure 2(a) in the first four positions. We therefore consider an important real world data test of our code to be passed.

Figure 2. Findings of our model for PWMs of $L_w=6$.

(a) The consensus motif that features a higher $logP\left(D\right)$ than (b). However, (b) is also present.

4.2 Other PWMs found for RBFOX2

The big overlap of the motifs Figure 2 suggests also searching for longer motifs, which may subsume the shorter ones. Indeed, the motif shown in Figure 3(d) contains both 6mers. Along with that we find local minima in the probability landscape, i.e. PWMs of $L_w\le 9$ (see Figure 3). All motifs have the consensus TGCATG in common. For RBFOX2, we found no evidence for our model to converge beyond $L_w=9$. The posterior probability $log\widehat{P}\left(D\right)$ – (2.5) at the optimized parameters – serves as a measure to compare and rank different motifs at equal $L_w$. The Bayesian Information Criterion (BIC)¹⁶ estimates the information content of every obtained local minimum,

(4.1)

$$\mathrm{BIC}=klogn-2log\widehat{P}\left(D\right),$$

with the number of degrees of freedom $k=\left(4-1\right)L_w+1$ (four nucleotides minus one for the normalization per position, one extra for $E_0$), and the number of data points $n={\sum }_S{n}_S\left(L_S-L_w+1\right)$, i.e. the number of possible binding sites in the entire foreground pool. Closely related is the Akaike Information Criterion (AIC)¹⁷

(4.2)

$$\mathrm{AIC}=2k-2log\widehat{P}\left(D\right),$$

which is a bit less susceptible to overfitting. Both criteria rank the longer PWMs as having the higher information content, while we rather expect to find an optimum in information content with respect to $L_w$. Therefore, we compare posterior probabilities only among equal $L_w$ in the following and leave the search for a comparison criterion among different $L_w$ for future work. Having found that the model retrieves the expected motif for a well-studied RBP, we sought to further explore its performance for others.

Figure 3. Non-consensus PWMs of different $L_w$ for RBFOX2.

4.3 CELF1

CELF1 is an RBP of the CUG-binding CELF family.¹⁸ CELF1 participates in multiple steps of post-transcriptional processing of RNAs, including splicing, translation and decay,¹⁹ and requires UGU motifs for high-affinity interaction with RNAs.²⁰ The corresponding Bind’n Seq dataset⁸ consists of libraries generated for seven different RBP concentrations, each containing $\sim 2\times {10}^7$ reads of $L_S=40$. Since 40 runs with completely random PWM intialization for $L_w=3,4,5$ did not yield any local optima of the probability landscape we decided to test whether the biased initialization of the PWM with the known motif ($\mathit{UGU}$/$\mathit{TGT}$), which was found as as enriched 3-mer in RNA Bind’n Seq⁸ enables the recovery of longer motifs. Indeed, our procedure yielded multiple extended versions of TGT Figure 4 with $L_w$ up to $8$. The context in which the reduced motif occurs is A/T-rich, the presence of an A immediately upstream indicates that CELF1 could recognize the AUG start codon. This could be interesting in light of CELF1 being a translational regulator of epithelial-mesenchymal transition via the binding of both cap-binding EIF4E and the poly(A)-binding protein.²¹

Figure 4. PWMs of different $L_w$ for CELF1.

4.4 HNRNPD

Within the class of heterogeneous ribonucleoproteins (hnRNPs), hnRNPD (also known as AUF1) is a well-known A/U-rich element RNA binding protein with important role in RNA decay.²² HNRNPD has been reported to bind clusters of AUUUA elements.²² The ENCODE-database⁸ lists AUAAU as another possible binding site for hnRNPD. While entirely random initializations do not deliver any convergent runs, we recover both AUAAU, and, with a smaller binding log-likelihood, AUUUA, as binding sites. Building on this shorter motifs enables the discovery of UAAAU-containg longer motifs that can be extended up to $L_w=14$ (see Figure 5), the highest length for which we found convergent results. We did find PWMs with $L_w=7,\dots ,13$ which we omitted in Figure 5 since they are parts of the two $L_w=14$ PWMs.

Figure 5. PWMs of different $L_w$ for HNRNPD.

The benchmark cases $L_w=5$ are show in (a) and (b), whereas (c) and (d) show the two motifs of $L_w=14$, the longest motifs that our algorithm found.

4.6 Other RBPs

There are other proteins covered in the Bind’n Seq data⁸ whose specificity was studied before. For example, we analyzed the data corresponding to MBNL1,²⁴ hnRNPL,²⁵ FUS,²⁶ TAF15.²⁷ For these, our model did not deliver any convergent results, even if the PWM was directly intialized with the expected consensus motif. This indicates that the enrichment did not work equally well for all the RBPs studied with the Bind’n Seq method. Interestingly, expected motifs were identified for these proteins with another method, the so-called kmer-enhancement that underlies most of the consensus motifs reported in the ENCODE database.⁸ Kmer enhancements are computed by counting the number of occurrences of every possible kmer in the foreground samples (RBP concentration $\ne 0$) and in the background samples (RBP concentration $=0$), and finally normalizing the foreground abundances by the background to extract the respective enhancement. The higher the enhancement of a given kmer, the higher the likelihood of it being bound by the RBP used in the experiment is thought to be. We computed these enhancements for 6mers, as done in the ENCODE studies. The results, shown in Figure 6 indicate that only RBFOX2 has a few highly enhanced 6mers with a clear hierarchy of enrichment, while all other investigated RBPs show a much flatter hierarchy of motif enhancements. An analysis of the Levenshtein distance of these motifs showed no clear difference in the pattern of distances among the leading motifs across the investigated RBPs. This suggests that these motifs correspond to many local minima of comparable depth, which precludes our algorithm finding clear PWMs representing the binding sites. Conversely, it becomes unclear whether the specificity of these RBPs would be well represented as weight matrices, or whether another model, for e.g. clusters of short, degenerate motifs may better represent the specificity of these RBPs.

Figure 6. Logarithmic (base eenhancement (counts foreground normalized by counts background) of all $4^6\approx 4000$ possible 6mers of all investigated RBPs, ranked by enhancement.

The top most enhanced motifs are RBFOX2: TGCATG, FUS: GCGCGC, hnRNPL: CACACA, MBNL1: GCTGCT, TAF15: GGGGGG, followed by variants thereof. All except the RBFOX2 motif are repeats of shorter oligomers.

5. Conclusion

We constructed a thermodynamical model that can be used to infer characteristic position weight matrices for the binding domains of RNA binding proteins from data obtained from affinity-based enrichment of oligonucleotides. Since we directly model the RBP-binding specificity as PWMs, our method bypasses arbitrary choices in the alignment of k-mers found to be individually enriched in the data. We evaluate our model on data in the public domain⁸ using expectation-maximization. For the benchmark case of RBFOX2, where very high-quality data is available, our model reproduces the known binding motif TGCATG where the first position features an almost uniform superposition of T and A in our result. Subdominantly, we find another PWM of $L_w=6$ as well as longer PWMs for RBFOX2. Unfortunately, our principled model does not robustly recover the binding motifs of other RBPs. For a few, e.g. CELF1, HNRNPD, HNRNPK, we can still recover the motifs as well as some longer variants, if the search starts from a PWM closely matching the expected motif. However, for most of the data sets our model did not deliver any prediction. Rather, other motifs, e.g. poly(A), often show higher enrichment in the data than the expected motifs. This indicates that experimental details that do not have to do with the affinity of the RBP for oligomers affect the frequency of oligomer capture in the data, complicating its analysis and raising questions about the biophysical realism of the motifs derived from the data. The motifs of the RBPs for which we did not recover a PWM tend to be more degenerate than those of RBPs for which some motif emerged. They consist of repeated occurrences of mono, di or trinucleotides. It is likely that for these motifs, it is crucially important to construct an appropriate background model. How to best do this remains to be determined in future work. Of note, while crosslinking and immunoprecipitation data is available for the proteins studied here, PWMs with enriched binding sites were also not recovered.²⁸ Thus, it will be interesting to explore models that allow more flexible spacing of RBP contact points on RNAs in the future.