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 ∝ exp ( ∑ i = 1 L w E i s i ) ≡ e E ( s ) , where s i is the nucleotide at position i , so s i ∈ { A , C , G , T } . Therefore, each element of the position weight matrix (PWM) M can be identified with m i α ≡ exp ( E i α ) , with columns being normalized as ∑ α m i α = 1 ∀ i = 1 , … , 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 ( E 0 ) . Combining these two possibilities, we arrive at the probability of a site s being bound by the RBP P ( bound | s , c , M , E 0 ) = c ( e E ( s ) + e E 0 ) 1 + c ( e E ( s ) + e E 0 ) , (2.1) 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 ( e E ( s ) + e E 0 0 ) ≪ 1 , we can linearize (2.1) P ( bound | s , c , M , E 0 ) ≈ c ( e E ( s ) + e E 0 ) . (2.2)

Consequently, the chance of an RBP being bound somewhere on a longer oligomer S with L S ≥ L w is P ( bound | S , c , M , E 0 ) = ∑ s ∈ S P ( bound | s , c , M , E 0 ) ≈ c ( e E ( S ) + ( L S − L w + 1 ) e E 0 ) , (2.3) where e E ( S ) ≡ ∑ s ∈ S e E ( s ) is the sum over all possible L w -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 P ( IP | S , c , M , E 0 ) = f S P ( bound | S , c , M , E 0 ) ∑ σ ∈ D f σ P ( bound | σ , c , M , E 0 ) ≈ f S ( e E ( S ) + ( L S − L w + 1 ) e E 0 ) ∑ σ ∈ D f σ ( e E ( σ ) + ( L σ − L w + 1 ) e E 0 ) , (2.4) with D being the data set of reads, IP being a binary variable indicating if the read is immuno-percipitated or not, and f S = 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 ( IP | S , M , E 0 ) is independent of the protein concentration c , which cancels out as an overall prefactor in both numerator and denominator. (2.4) is essentially a formulation of Bayes’ theorem, with P ( IP | S , M , E 0 ) being the likelihood of having a read S being immuno-percipitated, P ( S ) = f S being the likelihood 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 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 log P ( D | M , E 0 ) ≈ ∑ S ∈ D n S log [ f S ( e E ( S ) + ( L S − L w + 1 ) e E 0 ) ∑ σ ∈ D f σ ( e E ( σ ) + ( L σ − L w + 1 ) e E 0 ) ] , (2.5) where n S 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 K D 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 ⁹ [ RBP ] bound = n [ S ] K D [ RBP ] + 1 , (2.6) where K D 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. L S − L w + 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 , E 0 ) ∝ [ RBP ] bound , we can re-express (2.5) in terms of (2.6) as log P ( D | M , E 0 ) = C + ∑ S ∈ D n S log ( n [ S ] K D [ RBP ] + 1 ) = C + ∑ S ∈ D n S log ( f S ( L S − L w + 1 ) K D [ RBP ] + 1 ) , (2.7) where we used that the concentration [ S ] of any given oligomer S is proportional to f S , C is a proportionality constant independent of the binding specificity. Under the assumption that the binding is typically of low affinity, i.e. 1 K D [ RBP ] + 1 ≈ [ RBP ] K D , we arrive at log P ( D | M , E 0 ) ≈ C ′ + ∑ S ∈ D n S log ( [ RBP ] K D f S ( L S − L w + 1 ) ) . (2.8)

Hence, we can rank interactions by the dissociation constants relative to some reference PWM M ref and E 0 ref log ( K D K D ref ) = 1 ∑ S ∈ D n S [ log P ( D | M ref , E 0 ref ) − log P ( D | M , E 0 ) + ∑ S ∈ D n S log ( L S − L w + 1 L S − L w ref + 1 ) ] . (2.9)

As expected, the result corrects for the library size-dependence of P ( D | M , E 0 ) by dividing by the total foreground reads ∑ S ∈ D n S , and for differences in binding domain or read length by the last term. Assuming a random oligomer of length L w ref for which there is no unspecific binding ( E 0 ref → − ∞ ) as reference allows us to bring log P ( D | M rel ) into the simple closed form log P ( D | M ref , E 0 ) = ∑ S ∈ D n S log [ f S ( L S − L w ref + 1 ) ∑ σ ∈ D f σ ( L σ − L w ref + 1 ) ] . (2.10)

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 ref = 5 , a typical size for an RBP binding site. Relative K D ’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. 10). The number of possible 20mers is 4 20 ≈ 10 12 , much larger than the library sizes of ∼ 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 ∑ S ∈ B f S = 1 .

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: 1.

Initialize E 0 and the PWM elements m i α with respectively well-defined real numbers, i.e. E 0 ∈ ( − ∞ , 0 ] and ∑ α m i α = 1 ∀ i = 1 , … , L w . This can either be done in an entirely unbiased way or by pre-determining some motifs and initializing PWMs with values reflecting those motifs.

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

The log-values of relative K D ’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 K 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 × 10 9 vs. − 8.77 × 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).

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- K 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 × 10 9 and − 7.91 × 10 9 , 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 ( K D , E 0 , and log P ( D | 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 × 10 9 to − 1.43 × 10 9 .

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 D values for T-rich motifs are ≥ − 0.76 , with − 10.53 ≥ E 0 ≳ − 7.5 and log P ( D | M , E 0 ) ≤ − 2.81 × 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 D rel > − 0.33 and E 0 > − 5.49 , with log P ( D | M , E 0 ) ≤ − 4.04 × 10 9 . The motif with the highest dissociation constant obtained from random initializations is shown in Figure 6(C), has log K D rel = 0.98 , and corresponds to E 0 = − 0.64 and log P ( D | M , E 0 ) = − 5.68 × 10 9 .

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 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 × 10 9 to − 2.07 × 10 9 .

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 D rel = − 0.76 , E 0 = − 8.54 , and log P ( D | M , E 0 ) = − 2.09 × 10 9 . Mildly subdominant motifs are A-rich and less repetitive ( Figure 8(B), (C)) with log K D rel = − 0.76 , E 0 = − 7.67 , and log P ( D | M , E 0 ) = − 2.09 × 10 9 , as well as log K D rel = − 1.01 , E 0 = − 8.54 , and log P ( D | M , E 0 ) = − 2.09 × 10 9 .

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 D rel = − 0.81 ) with E 0 = − 13.26 and log P ( D | M , E 0 ) = − 2.99 × 10 9 . The lowest affinity motif ( Figure 9(B)) that random PWM-initialization yielded had log K D rel = 1.62 , E 0 = − 1.49 , log P ( D | M , E 0 ) = − 6.07 × 10 9 .

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 × 10 7 reads of L S = 40 .

Since 16 runs with completely random PWM intialization for L w = 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(A). The relative log-dissociation constants were 1.67 ≳ log K D rel ≳ 2.17 , unspecific binding was characterized by − 4.0 ≳ E 0 ≳ − 0.75 , and the likelihood of the daat was between log P ( D | M , E 0 ) − 8.41 × 10 9 and − 8.48 × 10 9 .

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 ≤ log K D rel ≤ 0.91 in log- K D , unspecific binding energy − 4.18 ≤ E 0 ≤ − 0.55 , and data likelihoods − 2.2 × 10 9 ≥ log P ( D | M , E 0 ) ≥ − 4.50 × 10 9 .

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 ≥ E 0 ≥ − 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 ≤ log K D rel ≤ 1.95 and data likelihoods − 2.23 × 10 9 ≥ log P ( D | M , E 0 ) ≥ − 2.30 × 10 9 . 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 L w > 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 E 0 = − 2.39 , log P ( D | M , E 0 ) = − 3.30 × 10 9 , and log K D rel = 1.59 . The penta(T) is also part of a higher affinity motif (see Figure 13(C)) with E 0 = − 5.90 , log P ( D | M , E 0 ) = 2.92 × 10 9 , and log K D rel = 0.00 . More complex motifs with higher affinity, but less specificity are displayed in Figure 13(A) with E 0 = − 2.05 , log P ( D | M , E 0 ) = − 2.56 × 10 9 , and log K D rel = − 0.10 and in Figure 13(B) with E 0 = − 1.32 , log P ( D | M , E 0 ) = − 2.56 × 10 9 , and log K D 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 ≈ 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- K D is relatively high, log K D rel = 0.11 , as is the contribution of unspecific binding, E 0 = − 0.67 . log P ( D | M , E 0 ) = − 7.44 × 10 9 .

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 D features E 0 = − 3.02 , log K D rel = 1.28 , and log P ( D | M , E 0 ) = − 1.15 × 10 9 . The motif corresponding to the highest K D had E 0 = − 3.86 , log K D rel = 2.14 , and log P ( D | M , E 0 ) = − 1.12 × 10 9 . 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.