We assume that an RBP binds an oligomer over a binding site s of length Lw and that the likelihood of binding taking place, according to Boltzmann’s law, goes as ∝exp(∑i=1LwEisi)≡eE(s) , where si is the nucleotide at position i , so si∈{A,C,G,T} . Therefore, each element of the position weight matrix (PWM) M can be identified with miα≡exp(Eiα) , with columns being normalized as ∑αmiα=1∀i=1,…,Lw .

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(E0) . Combining these two possibilities, we arrive at the probability of a site s being bound by the RBP

where the 1 in the denominator represents the (constant) chance of the RBP being unbound. Assuming that the binding of RBPs to oligomers is not saturated, i.e. c(eE(s)+(LS−Lw+1)eE0)≪1 , we can linearize (2.1)

Consequently, the chance of an RBP being bound somewhere on a longer oligomer S with LS≥Lw is

where eE(S)≡∑s∈SeE(s) is the sum over all possible Lw -mers s in S . The probability to observe each read S in the pool of oligomers that were selected by the interaction with the RBP is

with D being the data set of reads, IP being a binary variable indicating if the read is immuno-percipitated or not, and fS=P(S) a frequency prior that corrects for the fact that the pool of oligomers has a non-uniform nucleotide composition. Note that, due to the linearization in c , P(S|bound,c,M,E0) is independent of the protein concentration c , which cancels out as an overall prefactor in both numerator and denominator. Equation (2.4) is essentially a formulation of Bayes’ theorem, with P(S|bound,c,M,E0) being the likelihood of having a read S being immuno-percipitated, P(S)=fS being the prior probability of the read S in the input, and having an overall normalization in the denominator.

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

where nS is the number of copies of read S in the library.

The log-likelihood of binding derived in the previous section can be related to the dissociation constant KD known from the formalism of chemical reactions as follows. The concentration [RBP]bound of RBP bound to a multivalent oligomer with concentration [S] is given by¹¹

where KD is the dissociation constant and n is the number of binding sites on the oligomer. We make the simplification of considering binding sites to be in principle identical and their maximum number be given by the total number of configurations in which the RBP can contact the oligomer, i.e. LS−Lw+1 .

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(S|M,E0)∝[RBP]bound , we can re-express the log-likelihood of the oligomer pool from (2.5) using the formalism of (2.6):

where we used that the concentration [S] of any given oligomer S is proportional to fS , and C is a proportionality constant independent of the binding specificity. Under the assumption that the binding is typically of low affinity, i.e. 1KD[RBP]+1≈[RBP]KD , we arrive at

Hence, we can rank interactions by the dissociation constants relative to some reference PWM Mref and E0ref

As expected, the result corrects for the library size-dependence of P(D|M,E0) by dividing by the total foreground reads ∑S∈DnS , and for differences in binding domain or read length by the last term.

For the reference we can assume a random oligomer of length Lwref for which there is no unspecific binding ( E0ref→−∞ ), which would bring (2.5) in the form

Substituting (2.10) and (2.5) (for the latter using the parameter values given by the optimization procedure) into (2.9) 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 Lwref=5 , a typical size for an RBP binding site. Relative KD ‘s also provide a measure to rank different binding specificities, even of different lengths, for a given RBP.

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. 1). The number of possible 20 mers is 420≈1012 , much larger than the library sizes of ∼107 . 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. 12. 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 fS . 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 order d=4 , which seemed to give a good tradeoff between the accuracy of background read frequency prediction, size of Bind’n Seq libraries and available computational resources. The predicted frequency priors fS – and kmer-frequencies, accordingly – need to be normalized such that their sum over the background frequency pool B satisfies ∑S∈BfS=1 .

Having constructed our model with the final expression (2.5), as well as the background frequencies fS described in the subsection above, the remaining task is to identify the PWM and E0 that maximize the log-likelihood (c.f. (2.5)). To this end, we rely on the expectation maximization algorithm.^13,14 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.²

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 motif, but also the slight T bias at the preceding position 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.

The log-values of relative KD ‘s are in the range 0.48 - 2.39 , meaning that the two motifs shown in Figure 3(a) and (b) differ by ≈100 fold in their KD . The unspecific binding term in the run that yielded the motif from panel (a) was E0=−15.57 , while the one from the panel (b) run was E0=−2.49 . The overall log-likelihoods of the dataset were −8.51×109 vs. −8.77×109 , 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).

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.¹

Relative log- KD values ranged from −0.76 to 1.63 , the unspecific binding term was E0=−6.10 for Figure 4(a) and E0=−2.61 for Figure 4(b), and the log-likelihoods of the data were −3.80×109 and −7.91×109 , respectively.

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 ( KD , E0 , and logP(D|M,E0) ) 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×109 to −1.43×109 .

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- KD values for T-rich motifs are ≥−0.76 , with −10.53≥E0≳−7.5 and logP(D|M,E0)≤−2.81×109 . 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 logKDrel>−0.33 and E0>−5.49 , with logP(D|M,E0)≤−4.04×109 . The motif with the highest dissociation constant obtained from random initializations is shown in Figure 6(c), has logKDrel=0.98 , and corresponds to E0=−0.64 and logP(D|M,E0)=−5.68×109 .

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- KD ‘s range from 0.71 to 1.94 , unspecific binding strength from E0=−3.46 to E0=−2.24 , and the log-likelihood of the data from −1.95×109 to −2.07×109 .

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 logKDrel=−0.76 , E0=−8.54 , and logP(D|M,E0)=−2.09×109 . Mildly subdominant motifs are A-rich and less repetitive (Figure 8(b), (c)) with logKDrel=−0.76 , E0=−7.67 , and logP(D|M,E0)=−2.09×109 , as well as logKDrel=−1.01 , E0=−8.54 , and logP(D|M,E0)=−2.09×109 .

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) ( logKDrel=−0.81 ) with E0=−13.26 and logP(D|M,E0)=−2.99×109 . The lowest affinity motif (Figure 9(b)) that random PWM-initialization yielded had logKDrel=1.62 , E0=−1.49 , logP(D|M,E0)=−6.07×109 .

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 ∼2×107 reads of LS=40 .

Since 16 runs with completely random PWM intialization for Lw=5,6,…,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. The relative log-dissociation constants were 1.67≳logKDrel≳2.17 , unspecific binding was characterized by −4.0≳E0≳−0.75 , and the likelihood of the daat was between logP(D|M,E0)−8.41×109 and −8.48×109 .

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≤logKDrel≤0.91 in log- KD , unspecific binding energy −4.18≤E0≤−0.55 , and data likelihoods −2.2×109≥logP(D|M,E0)≥−4.50×109 .

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≥E0≥−2.48 , where Figure 12(b) corresponds to the highest affinity motif (Figure 12(a)) for which E0=−3.17 . The two motifs shown in Figure 12 come from a range of dissociation constants 0.42≤logKDrel≤1.95 and data likelihoods −2.23×109≥logP(D|M,E0)≥−2.30×109 . Direct initialization with the kmer-enrichment consensus ATAAT did not lead to any convergent and specific result.

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 Lw>5 in this data set, whether we used random initialization or shorter motif-guided initialization.

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 E0=−2.39 , logP(D|M,E0)=−3.30×109 , and logKDrel=1.59 . The penta(T) is also part of a higher affinity motif (see Figure 13(c)) with E0=−5.90 , logP(D|M,E0)=2.92×109 , and logKDrel=0.00 . More complex motifs with higher affinity, but less specificity are displayed in Figure 13(a) with E0=−2.05 , logP(D|M,E0)=−2.56×109 , and logKDrel=−0.10 and in Figure 13(b) with E0=−1.32 , logP(D|M,E0)=−2.56×109 , and logKDrel=−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 ≈50 -fold, which is quite substantial.

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- KD is relatively high, logKDrel=0.11 , as is the contribution of unspecific binding, E0=−0.67 . logP(D|M,E0)=−7.44×109 .

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 KD features E0=−3.02 , logKDrel=1.28 , and logP(D|M,E0)=−1.15×109 . The motif corresponding to the highest KD had E0=−3.86 , logKDrel=2.14 , and logP(D|M,E0)=−1.12×109 . Thus, compared to other RBPs (c.f. Figure 14), the specificity of PTBP3 remains to be better defined.

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 ( ∼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 ≠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.