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

2.1 Derivation of the data likelihood

We assume that an RBP binds an oligomer over a binding site $s$ of length $L_w$ and that the likelihood of 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 \{\mathrm{A},\mathrm{C},\mathrm{G},\mathrm{T}\}$. 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 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 of a site $s$ being bound by the RBP

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

Eventually, the log-likelihood of the entire library of oligomers $D$ explained by the specific binding to the RBP described by the PWM $M$ of length $L_w$ as well as unspecific binding described by parameter $E_0$ is given by

2.2 Relationship to the dissociation constant K_D

The log-likelihood of binding derived in the previous section can be related to the dissociation constant $K_{\mathrm{D}}$ known from the formalism of chemical reactions as follows. The concentration ${\left[\mathrm{RBP}\right]}_{\text{bound}}$ of RBP bound to a multivalent oligomer with concentration $\left[S\right]$ is given by⁹

Further assuming that the probability to observe individual oligomers in the sequencing pool is proportional to their relative abundance in complex with the RBP (which can only hold when at most one RBP molecule is bound to an individual oligomer), i.e. $P\left(S\right|M,E_0)\propto {\left[\mathrm{RBP}\right]}_{\text{bound}}$, we can re-express (2.5) in terms of (2.6) as

Hence, we can rank interactions by the dissociation constants relative to some reference PWM $M^{\mathrm{ref}}$ and $E_0^{\mathrm{ref}}$

As expected, the result corrects for the library size-dependence of $P\left(D\right|M,E_0)$ by dividing by the total foreground reads ${\sum }_{S\in D}{n}_S$, and for differences in binding domain or read length by the last term. Assuming a random oligomer of length $L_w^{\mathrm{ref}}$ for which there is no unspecific binding ($E_0^{\mathrm{ref}}\to -\infty$) as reference allows us to bring $logP\left(D\right|M^{\mathrm{rel}})$ into the simple closed form

Combining (2.9) with (2.5), the output of the optimization procedure, and the reference (2.10) allows us to compute the logarithm of the dissociation constant of the RBP-RNA binding for a representative RBP binding site relative to a random oligomer of length $L_w^{\mathrm{ref}}=5$, a typical size for an RBP binding site. Relative $K_{\mathrm{D}}$’s also provide a measure to rank different binding specificities, even of different lengths, for a given RBP.

3.1 Construction of the frequency priors f_S from a Markov model

RNA Bind’n Seq data does not only comprise libraries of reads pulled down by specific RBPs at non-vanishing RBP concentrations, but also libraries of input reads. The oligonucleotides that were used for RBP affinity-based selection were short, typically 20 nucleotides in length (c.f. Ref. 10). The number of possible 20mers is $4^{20}\approx {10}^{12}$, much larger than the library sizes of $\sim {10}^7$. Thus, even in the absence of selection ($c=0$), the expected overlap of two libraries is extremely small.

To preserve the statistical power of the foreground pool, i.e. use all the reads detected in the foreground sample in the analysis, even though they were not represented in the background sample, we have to predict the frequency of foreground reads under the assumption of no selection for binding the RBP. A commonly used approach for this type of problem is to train a Markov model from the background pool and construct the expected frequency of each read in the foreground from the trained model, just as in Ref. 11. For an completely unbiased process of oligomer synthesis, capture and sequencing the degree $d$ of the Markov model would be $0$, i.e. each base would be equally likely to occur at any position in the oligomer, and all 20mers would have the same prior frequency of occurrence $f_S$. However, various types of biases can lead to some sequences, with specific composition of short nucleotide motifs, being present in the data more often than others. In principle, the higher the degree of the Markov model, the larger the sequence context that can be resolved. However, an over-parametrization of the Markov-model needs to be avoided. As binding specificities of RBPs tend to be short and our main aim is to appropriately detect enrichment of motifs on this length scale, we used a Markov model of $d=4$, which seemed to give a good tradeoff between the accuracy of the background read frequency prediction, size of Bind’n Seq libraries, and available compuational resources. The predicted frequency priors $f_S$ – and kmer-frequencies, accordingly – need to be normalized such that their sum over the background frequency pool $B$ satisfies ${\sum }_{S\in B}{f}_S=1$.

3.2 Inferring PWMs by expectation maximization

Having constructed our model with the final expression (2.5), as well as the background frequencies $f_S$ described in the subsection above, the remaining task is to identify the PWM and $E_0$ maximize the log-likelihood (2.5). To this end, we rely on the expectation maximization algorithm.^12,13 Provided that only some of our model parameters can be directly inferred from the data, the algorithm optimizes the”hidden” parameters to maximize (2.5). The expectation-maximization procedure (EM) can be divided into the following steps:

Our code is written in C++ and python and is publicly available.¹⁵

4.1 RBFOX2

RBFOX2 is a key regulator of alternative splicing¹⁶ that was extensively studied with a variety of methods (e.g. Ref. 17). The RBFOX2 Bind’n Seq dataset¹⁰ consists in nine libraries obtained using 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. 18) 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 3 reproduce the importance of the GCATG morif, but also the slight T bias at the preceding position, see Figure 3(A). Moreover, we find the subdominant PWM Figure 3(B) which shares a CATG core with the highest affinity motif, but over-emphasizes the A bias downstream. These results demonstrate that our algorithm identifies the known consensus for RBFOX2.

Figure 3. Lowest (A) and highest (B) K_D PWMs for RBFOX2 obtained from our model with random initialization.

The highest affinity motif contains the known consensus GCATG, flanked by positions with low A/T bias. Both high and low affinity motifs share a CATG core.

The log-values of relative $K_{\mathrm{D}}$’s are in the range $0.48$ - $2.39$, meaning that the two motifs shown in Figure 3(A) and (B) differ by $\approx 100$ fold in their $K_{\mathrm{D}}$. The unspecific binding term in the run that yielded the motif from panel (A) was $E_0=-15.57$, while the one from the panel (B) run was $E_0=-2.49$. The overall log-likelihoods of the dataset were $-8.51\times {10}^9$ vs. $-8.77\times {10}^9$, respectively. Thus, the higher affinity motif corresponds to a higher peak in the likelihood landscape, and accounts for a more specific interaction. The binding specificity of RBFOX2 in Figure 3(A) is similar to that of the closely related RBFOX3 (c.f. Figure 2).

4.2 NOVA1

The neuro-oncological ventral antigen 1 (NOVA1) is a neuron-specific RBP¹⁹ found by SELEX experiments to bind a triple CAT-repeat.²⁰ This is confirmed by the motif identified by our algorithm (Figure 4). In contrast, the streaming kmer-enrichment analysis of RNA Bind’n Seq data in ENCODE predicted a single CAT.¹⁰

Figure 4. Lowest (A) and highest (B) K_D NOVA1 motifs obtained by our model with random initialization.

(A) shows the characteristic triple CAT-repeat.

Relative log-$K_{\mathrm{D}}$ values ranged from $-0.76$ to $1.63$, the unspecific binding term was $E_0=-6.10$ for Figure 4(A) and $E_0=-2.61$ for Figure 4(B), and the log-likelihoods of the data were $-3.80\times {10}^9$ and $-7.91\times {10}^9$, respectively.

4.3 PUM1

Pumilio homolog 1 is a well-studied member of the PUF famility of proteins, which binds to 3’UTRs of a variety of mRNAs and regulates their translation.²¹ While immunoprecipitation (IP) experiments found its binding specificity to be TGTAHATA,²² SKA analysis of RNA Bind’n Seq yielded TGTAH without further context.¹⁰ In contrast, we do find the full motif in the Bind’n Seq data (Figure 5). All obtained motifs are related to each other via shifts and show only gradual differences at single positions, which also reflects in characteristic numbers ($K_{\mathrm{D}}$, $E_0$, and $logP\left(D\right|M,E_0)$) being distributed only over a narrow range. The relative log-dissociation constants for Figure 5 are $0.36$, $0.36$, and $0.67$, where the former two differ by only $0.0004$. The corresponding unspecific binding strengths are $-3.47$, $-1.21$, and $-2.58$. The log likelihoods of the data range from $-1.39\times {10}^9$ to $-1.43\times {10}^9$.

Figure 5. Lowest (A), intermediate (B), and highest (C) K_D PWMs for PUM1 obtained from our model with random initialization.

The intermediate motif (B) represents the consensus.²²

4.4 PUF60

A representative from the same family of RBPs is the Poly(U)-binding splicing factor PUF60,²³ which regulates both transcription and translation.²⁴ Our algorithm finds the expected motif (poly(T)) as having the highest affinity Figure 6(A). The relative log-$K_{\mathrm{D}}$ values for T-rich motifs are $\ge -0.76$, with $-10.53\ge E_0\gtrsim -7.5$ and $logP\left(D\right|M,E_0)\le -2.81\times {10}^9$. Different variations (see Figure 6(B)) of the poly(T) motif are found, having in common the alternating pattern of very polarized – less polarized T (the latter blended with C). Poly(A)-containing motifs (c.f. Figure 6(C)) also appear at $log{K}_{\mathrm{D}}^{\mathrm{rel}}>-0.33$ and $E_0>-5.49$, with $logP\left(D\right|M,E_0)\le -4.04\times {10}^9$. The motif with the highest dissociation constant obtained from random initializations is shown in Figure 6(C), has $log{K}_{\mathrm{D}}^{\mathrm{rel}}=0.98$, and corresponds to $E_0=-0.64$ and $logP\left(D\right|M,E_0)=-5.68\times {10}^9$.

Figure 6. Lowest (A), intermediate (B) and highest (C) K_D PWMs for PUF60 obtained from our model with random initialization.

4.5 ELAVL4

The ELAV-like protein 4 (ELAVL4) is an RBP that is exclusively expressed in neurons.²⁵ It binds A/T-rich elements according to genome-wide IP experiments,^26,27 a pattern that we also observe in PWMs from our model (see Figure 7). Relative log-$K_{\mathrm{D}}$’s range from $0.71$ to $1.94$, unspecific binding strength from $E_0=-3.46$ to $E_0=-2.24$, and the log-likelihood of the data from $-1.95\times {10}^9$ to $-2.07\times {10}^9$.

Figure 7. Lowest (A), intermediate (B), and highest (C) K_D PWMs for ELVAL4 obtained from our model with random initialization.

4.6 HNRNPC

Turning towards another large familiy of RBPs, the heterogeneous nuclear ribonucleoproteins (hnRNPs), we highlight hnRNPC,²⁸ a protein that is involved pre-mRNA processing.²⁹ It is known to bind to poly(U) motifs, typically pentamers, found in both SKA analysis of RNA Bind’n Seq data¹⁰ and in electrophoretic mobility shift assays.³⁰ Our model does recover the dominant poly(T) motif at $log{K}_{\mathrm{D}}^{\mathrm{rel}}=-0.76$, $E_0=-8.54$, and $logP\left(D\right|M,E_0)=-2.09\times {10}^9$. Mildly subdominant motifs are A-rich and less repetitive (Figure 8(B), (C)) with $log{K}_{\mathrm{D}}^{\mathrm{rel}}=-0.76$, $E_0=-7.67$, and $logP\left(D\right|M,E_0)=-2.09\times {10}^9$, as well as $log{K}_{\mathrm{D}}^{\mathrm{rel}}=-1.01$, $E_0=-8.54$, and $logP\left(D\right|M,E_0)=-2.09\times {10}^9$.

Figure 8. Lowest (A), intermediate (B), and highest (C) K_D PWMs for hnRNPC obtained from our model with random initialization.

Figure 9. Lowest (A) and highest (B) K_D PWMs for PCBP1 obtained from our model with random initialization.

4.7 PCBP1

As a last example of RBPs of which we predict a motif that agrees with prior data, we highlight the poly (rC)-binding protein 1 (PCBP1).³¹ Its binding specificity consists of a few very polarized C’s linked by A/T-enriched sequence.³² This pattern is clearly reproduced by our highest affinity motif in Figure (9)(A) ($log{K}_{\mathrm{D}}^{\mathrm{rel}}=-0.81$) with $E_0=-13.26$ and $logP\left(D\right|M,E_0)=-2.99\times {10}^9$. The lowest affinity motif (Figure 9(B)) that random PWM-initialization yielded had $log{K}_{\mathrm{D}}^{\mathrm{rel}}=1.62$, $E_0=-1.49$, $logP\left(D\right|M,E_0)=-6.07\times {10}^9$.

4.8 CELF1

CELF1 is an RBP of the CUG-binding CELF family,^33,34 that participates in multiple steps of post-transcriptional processing of RNAs, including splicing, translation and decay.³⁵ CELF1 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 $16$ runs with completely random PWM intialization for $L_w=5,6,\dots \mathrm{,14,15}$ did not yield any local optimum of the probability landscape we decided to test whether the biased initialization of the PWM with the known motif (TGT), which was found as as enriched 3-mer in RNA Bind’n Seq¹⁰ enables the recovery of longer motifs. However, our model predicted only poly(A) stretches as convergent optima of the probability landscape, see Figure 10(A). The relative log-dissociation constants were $1.67\gtrsim log{K}_{\mathrm{D}}^{\mathrm{rel}}\gtrsim 2.17$, unspecific binding was characterized by $-4.0\gtrsim E_0\gtrsim -0.75$, and the likelihood of the daat was between $logP\left(D\right|M,E_0)-8.41\times {10}^9$ and $-8.48\times {10}^9$.

4.9 EWSR1

EWSR1, or Ewing’s sarcoma protein, forms fusions with a number of other proteins and serves as a transcriptional activator in human solid tumors like Ewing’s sarcoma and malignant melanoma.³⁷ Information about its binding specificity is scarce, while SKA analysis of RNA Bind’n Seq data predicts G-rich elements.¹⁰ Our investigation lead to different varieties of poly(A) motifs (Figure 11(A) and (B)), which vary in relative dissociation constants $-0.79\le log{K}_{\mathrm{D}}^{\mathrm{rel}}\le 0.91$ in log-$K_{\mathrm{D}}$, unspecific binding energy $-4.18\le E_0\le -0.55$, and data likelihoods $-2.2\times {10}^9\ge logP\left(D\right|M,E_0)\ge -4.50\times {10}^9$.

Figure 10. Lowest (A) and highest (B) K_D PWMs for CELF1 predicted by our model from random initialization.

Enrichment analysis suggests a prominent $\mathit{TGT}$ as consensus, which our model cannot recover, even when initializing with $\mathit{\text{NTGTN}}$.

Figure 11. Lowest (A) and highest (B) K_D PWMs for EWSR1 predicted by our model from random initialization.

4.10 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. We find poly(A) stretches of different lengths as convergent and specific binding motifs (see Figure 12). While A is the dominating nucleotide at every position, it is followed consistently by T, which fits the reported A/U-rich binding domains in the literature. We find unspecific binding terms in the range $-7.56\ge E_0\ge -2.48$, where Figure 12(B) corresponds to the highest affinity motif (Figure 12(A)) for which $E_0=-3.17$. The two motifs shown in Figure 12 come from a range of dissociation constants $0.42\le log{K}_{\mathrm{D}}^{\mathrm{rel}}\le 1.95$ and data likelihoods $-2.23\times {10}^9\ge logP\left(D\right|M,E_0)\ge -2.30\times {10}^9$. Direct initialization with the kmer-enrichment consensus ATAAT did not lead to any convergent and specific result.

Figure 12. Lowest (A) and highest (B) K_D PWMs for hnRNPD predicted by our model from random initialization.

Enrichment analysis suggests a quite unpolar A/T-rich binding domain. Both motifs, (A) and (B), show A as the dominant nucleotide, followed by T at each position.

4.11 HNRNPK

We were also interested in determining whether we can recover G/C-rich binding motifs from the data and therefore applied the model to heterogeneous nuclear ribonucleoprotein K (hnRNPK), a member of the poly(C) binding family of proteins.³⁹ We could only recover one of the two consensus motifs reported in the ENCODE analysis of these data (GCCCA, from SKA)¹⁰ when specifically initializing the PWM with this motif. The second reported motif, with the CACGC consensus, could not be found by our algorithm even when the PWM was initialized with the motif itself and even when sequences containing the first motif were eliminated, indicating that this motif does not correspond to a local maximum of the likelihood function. We did not find any PWMs of $L_w>5$ in this data set, whether we used random initialization or shorter motif-guided initialization.

4.12 CPEB1

The cytoplasmic polyadenylation element-binding protein 1 (CPEB1)⁴⁰ serves as a translational regulator by binding specific U-rich sequences in 3’UTRs inducing cytoplasmic adenylation. The streaming kmer algorithm predicts a poly(T) motif from RNA Bind’n Seq data, which is what we recover as the lowest affinity motif in Figure 13(D) with $E_0=-2.39$, $logP\left(D\right|M,E_0)=-3.30\times {10}^9$, and $log{K}_{\mathrm{D}}^{\mathrm{rel}}=1.59$. The penta(T) is also part of a higher affinity motif (see Figure 13(C)) with $E_0=-5.90$, $logP\left(D\right|M,E_0)=2.92\times {10}^9$, and $log{K}_{\mathrm{D}}^{\mathrm{rel}}=0.00$. More complex motifs with higher affinity, but less specificity are displayed in Figure 13(A) with $E_0=-2.05$, $logP\left(D\right|M,E_0)=-2.56\times {10}^9$, and $log{K}_{\mathrm{D}}^{\mathrm{rel}}=-0.10$ and in Figure 13(B) with $E_0=-1.32$, $logP\left(D\right|M,E_0)=-2.56\times {10}^9$, and $log{K}_{\mathrm{D}}^{\mathrm{rel}}=-0.10$. While the difference in affinity between the first three motifs shown in Figure 13(A-C) is negligible, while the penta(T) motif alone (Figure 13(D) differs from the highest affinity one by $\approx 50$-fold, which is quite substantial.

Figure 13. Lowest (A), intermediate (B) and (C), and highest (D) K_D PWMs for CPEB1 predicted by our model from random initialization.

Although the consensus binding motif is recovered in (C) and (D), we predict (A) and (B) to be stronger binding specificities.

4.13 EIF4H

As a representative of the family of eukaryotic translation initiation factors, we highlight the eukaryotic translation initiation factor 4H (EIF4H).⁴¹ SKA enrichment studies of RNA Bind’n Seq data predict a poly(G) motif for this protein,¹⁰ but we did not find corroborating data obtained with other approaches. Our model identified only one convergent and specific motif shown in Figure 14, which is pyrimidine-rich. The relative log-$K_{\mathrm{D}}$ is relatively high, $log{K}_{\mathrm{D}}^{\mathrm{rel}}=0.11$, as is the contribution of unspecific binding, $E_0=-0.67$. $logP\left(D\right|M,E_0)=-7.44\times {10}^9$.

Figure 14. Single convergent and specific PWM for EIF4H predicted by our model from random initialization.

4.14 PTBP3

Polypyrimidine tract-binding protein 3 (PTBP3) is involved in the regulation of numerous steps of protein production, i.e. splicing, alternative 3’ end processing, mRNA stability and RNA localization.⁴² Data on PTBP3 binding specificity is lacking, and enriched kmer was reported by the SKA from RNA Bind’n Seq data.¹⁰ Our model predicts different variations of poly(A), see Figure 15. The predicted binding specificities with the lowest $K_{\mathrm{D}}$ features $E_0=-3.02$, $log{K}_{\mathrm{D}}^{\mathrm{rel}}=1.28$, and $logP\left(D\right|M,E_0)=-1.15\times {10}^9$. The motif corresponding to the highest $K_{\mathrm{D}}$ had $E_0=-3.86$, $log{K}_{\mathrm{D}}^{\mathrm{rel}}=2.14$, and $logP\left(D\right|M,E_0)=-1.12\times {10}^9$. Thus, compared to other RBPs (c.f. Figure 14), the specificity of PTBP3 remains to be better defined.

Figure 15. Lowest (A) and highest (B) K_D PWMs for PTBP3 predicted by our model from random initialization.

4.15 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.⁴⁶

Random and consensus intialization did result in convergent and specific binding specificities for FUS, however, these were poly(A) motifs which do not agree with the GGUG consensus from SELEX experiments.⁴⁷ For the other threee proteins our model did not deliver any convergent results, even when the PWM was directly initialized 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 or that our method does not identify well motifs that are very short ($\sim 2$nts, observed in the SKA enrichment analysis). In this analysis, kmer enrichments 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 enrichment. The higher the enrichment 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 enrichments for 6mers, as done in the ENCODE studies. The results, shown in Figure 16 indicate that only RBFOX2 has a very clear hierarchy of top enriched 6mers, while the other investigated RBPs show a much flatter hierarchy of motif enrichments. 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 16. Logarithmic (base $e$) enrichment (count in foreground reads normalized by count in background) for all $4^6=4096$ possible 6mers in all datasets (corresponding to different RBPs), ranked by enrichment.

The top most enriched motifs are for RBFOX2: TGCATG, FUS: GCGCGC, hnRNPL: CACACA, MBNL1: GCTGCT, TAF15: GGGGGG. All except the RBFOX2 motif are repeats of shorter oligomers.