Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

De Novo Prediction of PTBP1 Binding and Splicing Targets Reveals Unexpected Features of Its RNA Recognition and Function

De Novo Prediction of PTBP1 Binding and Splicing Targets Reveals Unexpected Features of Its RNA... The splicing regulator Polypyrimidine Tract Binding Protein (PTBP1) has four RNA binding domains that each binds a short pyrimidine element, allowing recognition of diverse pyrimidine-rich sequences. This variation makes it difficult to evaluate PTBP1 binding to particular sites based on sequence alone and thus to identify target RNAs. Conversely, transcriptome-wide binding assays such as CLIP identify many in vivo targets, but do not provide a quantitative assessment of binding and are informative only for the cells where the analysis is performed. A general method of predicting PTBP1 binding and possible targets in any cell type is needed. We developed computational models that predict the binding and splicing targets of PTBP1. A Hidden Markov Model (HMM), trained on CLIP-seq data, was used to score probable PTBP1 binding sites. Scores from this model are highly correlated (r =20.9) with experimentally determined dissociation constants. Notably, we find that the protein is not strictly pyrimidine specific, as interspersed Guanosine residues are well tolerated within PTBP1 binding sites. This model identifies many previously unrecognized PTBP1 binding sites, and can score PTBP1 binding across the transcriptome in the absence of CLIP data. Using this model to examine the placement of PTBP1 binding sites in controlling splicing, we trained a multinomial logistic model on sets of PTBP1 regulated and unregulated exons. Applying this model to rank exons across the mouse transcriptome identifies known PTBP1 targets and many new exons that were confirmed as PTBP1-repressed by RT-PCR and RNA-seq after PTBP1 depletion. We find that PTBP1 dependent exons are diverse in structure and do not all fit previous descriptions of the placement of PTBP1 binding sites. Our study uncovers new features of RNA recognition and splicing regulation by PTBP1. This approach can be applied to other multi-RRM domain proteins to assess binding site degeneracy and multifactorial splicing regulation. Citation: Han A, Stoilov P, Linares AJ, Zhou Y, Fu X-D, et al. (2014) De Novo Prediction of PTBP1 Binding and Splicing Targets Reveals Unexpected Features of Its RNA Recognition and Function. PLoS Comput Biol 10(1): e1003442. doi:10.1371/journal.pcbi.1003442 Editor: Roderic Guigo, Center for Genomic Regulation, Spain Received April 1, 2013; Accepted November 23, 2013; Published January 30, 2014 Copyright:  2014 Han et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: This work was supported by NIH Grant R01 GM49662 and an associated ARRA supplement to DLB. DLB is an Investigator of the Howard Hughes Medical Institute. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. * E-mail: [email protected] not currently feasible, in part due to our incomplete understanding Introduction of RNA recognition by the splicing regulators and their Alternative splicing of pre-mRNA commonly determines the mechanisms of action. protein output of mammalian genes, with most genes generating Whole-transcriptome crosslinking methods for individual pro- multiple mRNA and protein products [1]. A typical alternative teins in vivo are allowing the identification of large numbers of exon is affected by multiple pre-mRNA binding proteins that may protein/RNA interaction sites [8,9,10,11]. These data can be either enhance or repress splicing [2]. The expression and activity overlapped with functional data on splicing to identify possible of these splicing regulatory proteins can vary with development, direct target exons for particular proteins [12,13,14,15]. However, cell type, or cellular stimulus [3]. This complex combinatorial there are limitations in the interpretation of these data. Cross- regulation can be seen in the conserved sequences within and linking efficiency can vary between different proteins and between surrounding alternative exons, which generally contain the individual binding sites, making it difficult to relate the crosslinking binding sites for many different regulators. These sequences make signal to the actual binding affinity. These signals are also up what is sometimes called the splicing code as they determine dependent on the expression of the bound RNA, and since these where and when the exon is spliced into an mRNA [4,5,6,7]. Such data are generated one tissue or cell type at a time it is not always a code should allow the development of models that predict exon feasible to extend the results from one setting to a new cell type or regulation based solely on the RNA binding affinity of the many point in development. It would be extremely useful to be able to regulatory proteins and their other interactions. However, this is scan for binding affinity across the complete transcriptome and to PLOS Computational Biology | www.ploscompbiol.org 1 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets elements for each domain are not known. Moreover, the flexible Author Summary linkers separating some of the RRM domains and the requirement A key step in the regulation of mammalian genes is the for a gap between elements simultaneously bound to domains splicing of the messenger RNA precursor to produce a three and four allow for substantial degeneracy in PTBP1 binding mature mRNA that can be translated into a particular sites. This degeneracy and the lack of understanding of the protein needed by the cell. Through the process of sequence features that contribute to binding affinity have made it alternative splicing, mRNAs encoding different proteins difficult to identify PTBP1 binding sites based on sequence alone, can be derived from the same primary gene transcript. The and to assess which sequences surrounding an exon might regulation of this process plays essential roles in the contribute to PTBP1 regulation. development of differentiated tissues and is mediated by Experiments with model substrates indicate that a single high special pre-mRNA binding proteins. To understand how affinity PTBP1 binding site placed upstream of an exon, or within these proteins control gene expression, one must charac- it, can repress splicing [27]. However, strong repression of an terize what they recognize in RNA and identify these efficiently spliced exon requires an additional binding site either binding sites across the genome in order to predict their within the exon or downstream from an exon with an upstream targets. Models that allow this prediction are essential to high affinity site [17,27,28]. PTBP1 is also known to enhance the understanding developmental regulatory programs and splicing of certain exons [13,19,20]. The properties of these exons their perturbation by disease causing mutations. In this and how they differ from those that are repressed by PTBP1 are study, we use statistical methods to build models of RNA unclear, with different studies coming to different conclusions recognition by the important splicing regulator PTBP1 and [13,19]. An analysis of CLIP data in HeLa cells found that PTBP1 then apply these models to predict PTBP1 regulation of new gene transcripts. We show that PTBP1 has different sites near the adjacent constitutive exons could enhance the specificity for RNA than was previously recognized and inclusion of an alternative exon between them [13]. In contrast, that its target exons are more diverse than was known examination of exons whose splicing was reduced by double before. There are many similar splicing regulators in knockdown of Ptbp1 and Ptbp2 found that they frequently had mammalian cells, and these analyses provide a general binding sites immediately downstream [19], whereas splicing framework for the computational analysis of their RNA repression often involved upstream binding sites: a pattern binding and target identification. observed for other splicing regulators. These results are not mutually exclusive. It is possible that the two groups examined predict exon targets in tissues that have not yet been subjected to different subsets of the many exons regulated by PTBP1, and that experimental analysis. the protein may show additional patterns of protein binding Splicing regulatory proteins commonly contain multiple RRM or adjacent to its target exons. other RNA binding domains, with each domain recognizing a short In this study we sought to understand the sequence features that element of a few nucleotides [2,16]. Subtle variation in the optimal determine RNA binding by PTBP1 and to examine how they are binding element of each domain and flexible peptide linkers combined in exons that are targeted by the protein. We first between them allow for significant degeneracy within high affinity developed a statistical model of PTBP1 binding sites that identifies binding sites. Although the short sequence motifs that are common new features of RNA recognition by the protein. This binding to a set of binding sites are readily identified, these likely constitute model was then applied to the assessment of exon regulation by only a portion of a full high affinity site. To rank binding sites and PTBP1 across the transcriptome. assess their finer structures, we need an approach to search for clusters of these short motifs and to score for binding affinity. Results The Polypyrimidine tract binding protein 1 (PTBP1) is a widely G containing triplets contribute to PTBP1 binding studied splicing regulatory protein [17,18]. PTBP1 is known to repress the splicing of a large number of exons by binding in their To examine the interactions of PTBP1 across many binding adjacent introns or within the exons themselves. PTBP1 is down sites, we used a set of PTBP1-bound sequences identified by regulated in differentiating neurons and muscle cells to allow crosslinking immunoprecipitation (CLIP) [13]. PTBP1 has four inclusion of PTBP1 repressed exons during development of these RRMs separated by linker peptides, with each RRM recognizing a tissues [19,20,21]. In neurons the loss of PTBP1 is accompanied by pyrimidine triplet. In previous studies we found that a minimal the up-regulation of the homologous protein PTBP2 [17,20,22]. high affinity binding site for the protein extended across 25 to 30 PTBP2 has similar binding properties to PTBP1 and represses nucleotides, about the average size of the CLIP clusters (29 nt) some of the same exons [23]. Other exons are more sensitive to [27]. Given the triplet recognition and the need for spacers PTBP1 than PTBP2 and are induced to splice when PTBP2 between the direct RRM contacts, it is unlikely that every replaces PTBP1 in early neurons [24]. nucleotide within a CLIP cluster makes a direct base-specific contact with the protein or otherwise contributes to binding PTBP1 contains four RRM domains that recognize short pyrimidine elements [25]. Flexible linkers separate RRM domains affinity. This information about direct binding is hidden in the one and two, and domains two and three. RRM domains three examination of a CLIP tag, but should affect the triplet frequencies and four interact through a hydrophobic interface that position within the entire set of tags. We designed a two-state Hidden their RNA binding surfaces on opposite faces of the two-domain Markov Model (HMM) based on triplets to assess whether triplets structure. This orientation requires that the RNA elements would segregate into two states and whether these two states interacting with the structure be separated by an RNA loop differed in their PTBP1 binding or non-binding potential. The [26]. The structure of each of the PTBP1 RRM domains has been 48,604 CLIP clusters from the human transcriptome were solved in complex with the hexanucleotide, CUCUCU [25]. extracted and used to train the HMM (Figure 1A) [29,30]. This These structures show each domain binding a nucleotide triplet training defined two states showing distinctly different triplet with some additional contacts, and making similar base specific distributions (Figure 1B). Pleasingly, all of the pyrimidine triplets interactions with CU or UC dinucleotides. Other sequences can segregated into State 1. We called this state the PTBP1 binding likely make different base specific contacts, and the optimal state, as we confirm below. We found that 20 triplets have higher PLOS Computational Biology | www.ploscompbiol.org 2 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 1. PTBP1 binding model. A. Scheme of the PTBP1 binding model. The two-state HMM model was trained on PTBP1 bound RNA sequences (48,604 clusters) from published PTBP1-CLIP experiments. Triplets from these CLIP clusters were predictive of two states, with all of the pyrimidine triplets preferred by State 1. The diagram presents the structure of the PTBP1 HMM (Hidden Markov Model) and its trained transition probabilities. B. The probabilities that triplets are seen states 1 or 2 (emission probabilities) are plotted in black and gray bars, respectively. Asterisks indicate G containing pyrimidine triplets. doi:10.1371/journal.pcbi.1003442.g001 probabilities to be seen in the PTBP1 binding state. All triplets of its RRM domains can presumably make specific contacts with containing only pyrimidines were included in this 20-triplet set G residues. On the other hand, all A containing triplets have (Figure 1B), with the top-scoring triplet UCU showing the modest positive emission probabilities for state 2 and are likely to alternating C and U nucleotides seen in many characterized be either neutral or to inhibit PTBP1 binding. PTBP1 binding sites. We next tested the HMM scoring, which strongly weights the Interestingly, multiple triplets containing G residues are also triplets from state 1 over state 2, for prediction of PTBP1 binding. preferred in State1 (Figure 1B). These triplets often contain U We performed cross validation experiments on the Hela CLIP residues as the other nucleotides. Some of these triplets, such as dataset. A background dataset was generated using ten randomly UGU, have output (emission) probabilities in State 1 that are picked sequences from each gene identified as containing a CLIP similar to pyrimidine triplets, presumably also making them cluster. Applying the model to this data set gave us a distribution of predictive of PTBP1 binding. In contrast, triplets containing A scores that was compared to scores generated by subsets of the residues, even if the other two nucleotides are pyrimidines, were all CLIP clusters removed from the training set prior to training. As preferred by the non-PTBP1 binding State 2. These results shown in Figure S1, sequences from subsets of the CLIP clusters indicate that PTBP1 is not strictly pyrimidine specific. At least one scored significantly higher than background. PLOS Computational Biology | www.ploscompbiol.org 3 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets We also tested our model on an independent iCLIP dataset this background model generated scores that predicted affinity from human embryonic stem cells (ESC) (Figure S2). Unlike reasonably well. However, it did generate negative scores for a standard CLIP, iCLIP tags define the probable crosslink site as couple of probes that are shown to bind (data not shown). Thus, the uniform model gave the most accurate scoring of the being the 59 terminus of the tag. We used a Viterbi algorithm to predict the most probable state path predicted by the PTBP1 background models we tested. HMM model for each iCLIP tag. Defining triplets from the State1 The data demonstrate that HMM scoring based on triplet (PTBP1 binding) and triplets from State 2 (nonbinding), we found frequencies can accurately predict the observed binding affinities that the frequency of predicted binding triplets is highly enriched across a wide range of Kd values (from ,250 nM to 1 nM). Probe in the iCLIP cluster regions and peaks precisely at the crosslink 6 yields a z-score of 0.82 and binds with a Kd of 257 nM, whereas site. This indicates that State 1 probability is highly associated with probe 10 scores 2.74 in the model and binds with a Kd of 73 nM PTBP1 crosslinking in vivo. (Figure 2B). These sequences include G containing triplets that To more quantitatively assess the relationship between the contribute to the binding scores. This method allows any sequence to now be quantitatively assessed for possible PTBP1 binding, HMM score and RNA binding, we applied the trained model to a set of 100,000 random 69 nucleotide sequences. This length allows which was not previously possible by simply looking for clusters of a limited number of motifs. This HMM based approach should be for one hexanucleotide binding site for each of the four RRMs applicable to the prediction of binding sites and affinity for other with 15 nucleotide gaps, the minimum gap required for multi-domain RNA binding proteins. simultaneous binding by RRMs 3 and 4 [25,26]. The scores are calculated as a log-odds ratio of the probabilities of the sequence having been generated by the HMM over a background model Placement of PTBP1 binding sites adjacent to target that assigns equal probability to all triplets. The random sequences exons generated a distribution of scores that was used to normalize the With our new method of defining PTBP1 binding sites, we next binding scores, with the average score for random sequence set to examined PTBP1 target exons for the location of predicted PTBP1 zero, and the z-score defined as the deviation from the average as binding. In part, we wanted to reassess two previous studies that shown in Figure S3A [29]. Thus a sequence with a z-score of 2.74 came to differing conclusions regarding the placement of PTBP1 is 2.74 standard deviations from the average (empirical sites adjacent to its target exons. One group mapped PTBP1 CLIP p-value = 0.005), and is predicted to be a significantly stronger clusters adjacent to a limited number of PTBP1 repressed and binder than the average sequence (500 of the 100,000 random enhanced exons [13]. This study described PTBP1 repressed sequences have scores equal or greater than this sequence). A exons as enriched for binding sites both upstream and down- negative z-score is predicted to bind less well than the average stream, as has been seen in studies of individual exons. They did sequence. We isolated thirteen sequences from the mouse not observe PTBP1 CLIP clusters within repressed exons, even transcriptome that exhibited a range of scores from 22.62 to + though such exons have been described [17,32,33]. The PTBP1 4.40 (Figure 2A). These were transcribed in vitro and subjected to enhanced exons they examined showed a trend in PTBP1 binding electrophoretic mobility shift assay to measure binding to near the flanking constitutive exons. A second study examined recombinant PTBP1 (Figure 2B; Figure S3B). Sequences yielding exons showing altered splicing on splicing-sensitive microarrays negative scores all failed to bind PTBP1 within the protein after Ptbp1/Ptbp2 double knockdown [19]. CLIP clusters derived concentration range tested, with the exception of probe 4, which from the first study were mapped to these exons. The authors bound weakly, below the level that would allow measurement of found CLIP cluster enrichment upstream and within PTBP1/ an affinity constant. Positive scoring sequences all yielded PTBP1 PTBP2 repressed exons. In contrast to the previous study, they bound complexes that were assayable by gel shift to derive found that PTBP1/PTBP2 enhanced exons showed enrichment apparent binding affinities. The apparent Kds of these RNAs for CLIP tags in the downstream region. This pattern of binding showed a very strong negative correlation with their binding score site placement relative to repressed and enhanced exons has been from the model (Pearson correlation coefficient =2 0.9), where a observed for several other splicing regulatory proteins [14,34]. higher score predicts a lower Kd and hence a higher affinity In our study, we defined four groups of exons from a set of (Figure 2A). Thus, the scoring system performed very well in exons previously assessed for splicing after Ptbp1 knockdown predicting PTBP1 binding affinity. [20,35]. These included 68 PTBP1-repressed exons whose splicing Two sequences (probes 9 and 11) showed variable binding that increases after Ptbp1 knockdown, 37 PTBP1-enhanced exons shifted their Kd’s slightly off the fitted curve relating z-score to Kd. whose splicing decreases after knockdown, 69 control exons that These may have secondary structures that reduce binding affinity are not affected by Ptbp1 depletion but are known to be thus increase their apparent Kd. To look at this, we examined the alternatively spliced (PTBP1-non regulated), and 1,000 constitu- predicted structure of each probe using the RNA fold program tive exons. We determined the density of predicted PTBP1 [31]. Probes 9 and 11 did not show an overall free energy of binding states within a 24-nucleotide window sliding along the folding substantially lower than other RNAs. However, it is exon region. We also examined the sequence encompassing the difficult to rule out that they contain a local structure that adjacent constitutive exons (Figure 3A). As expected, the non- sequesters some key feature for PTBP1 recognition. regulated control and constitutive exon sets did not exhibit high In addition to the background model using uniform triplet probabilities of PTBP1 binding except in the polypyrimidine tract frequencies, we also tested control sequence sets using different of the 39 splice site. On the other hand, the introns upstream of nucleotide frequencies (Figure S4). Control sets that maintain the PTBP1 repressed exons show enrichment of potential PTBP1 mono or dinucleotide frequencies of the PTBP1 CLIP tags while binding sites starting from 250 nucleotides upstream of the exon. shuffling the triplet frequencies did not perform well. This is not Relative to the control exons, exons repressed by PTBP1 also surprising because these sequences are highly skewed in nucleotide exhibited substantial enrichment of PTBP1 binding sites within the content and the shuffling does not change the triplet frequencies exon itself and within the first 100 nucleotides of the downstream dramatically. We also tested a background model based on intron. The repressed exons thus exhibit binding site placement random sequences selected from genes containing PTBP1 CLIP that combines the findings of the two previous studies [13,19]. The clusters (ten sequences from each gene). Like the random dataset, PTBP1-enhanced exon set also shows enrichment of PTBP1 PLOS Computational Biology | www.ploscompbiol.org 4 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 2. Validation of the PTBP1 binding model. A. To validate binding scores, thirteen RNAs with various PTBP1 binding scores were transcribed in vitro and subjected to binding assay. Apparent Kd’s (dissociation constant) were highly negatively correlated with PTBP1 binding scores (Pearson correlation =20.9). B. Four RNA sequences with predicted PTBP1 binding scores (Full data binding data in Figure S3). Potential PTBP1 binding sites are underlined and in bold. Experimental binding affinities were assessed by electrophoretic mobility shift of RNA by PTBP1 and compared with prediction scores. Apparent dissociation constants (Kd) were defined as the concentration at which half the protein was bound to RNA. doi:10.1371/journal.pcbi.1003442.g002 binding sites within the downstream intron relative to control (Figure 3A). Similar to what was seen in the previous study by exons, although the distribution of binding sites across this region Llorian, we found little enrichment of PTBP1 sites within was different between the repressed and enhanced exon sets enhanced exons [19]. There is a limited enrichment adjacent to PLOS Computational Biology | www.ploscompbiol.org 5 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 3. Sequence characteristics of PTBP1-dependent alternatively spliced exons. A. An RNA map shows enrichment of predicted PTBP1 binding sites near PTBP1-dependent exons. The Y-axis plots average density of predicted PTBP1 binding states within a 24 nt window; the length of overlap between two adjacent windows was 8 nt. B. To assess PTBP1 binding signatures of individual exons, known PTBP1 regulated exons were clustered by their PTBP1 binding score profiles and visualized as heat maps. These heat maps indicate wide variation in the positions of PTBP1 binding sites between individual exons. C. Four sequence features including the PTBP1 binding scores and 39 splice site strength show statistically significant differences between regulated and control exon groups (one-tailed Student’s t-tests). doi:10.1371/journal.pcbi.1003442.g003 PLOS Computational Biology | www.ploscompbiol.org 6 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets the exons flanking enhanced exons. Interestingly however, we find to consider many different factors. Nevertheless, models based on some PTBP1 enhanced exons that have PTBP1 binding sites single factors will be useful for understanding the relative upstream of the exon. These were not seen in either previous contributions of individual proteins to patterns of splicing study. Our results are generally consistent with the known regulation. Such models will be easier to interpret regarding the placement of PTBP1 binding sites in PTBP1 target exons and contributions of individual factors to individual exons than more imply that rules correlating the position of PTBP1 binding to its complex models. Moreover in the longer term, models developed effect on a target exon are not as strict as seen for some other for different individual factors can be combined to make more splicing regulators. The mechanisms proposed from previous maps accurate predictions. To assess how well one might model splicing of PTBP1 binding do not appear to be generalizable to all PTBP1 regulation by a single factor, we examined whether the strength targets [13,19,27]. and placement of predicted PTBP1 binding sites could be used to Binding maps for PTBP1 and other splicing regulators show the predict new PTBP1 dependent exons. We plotted the scores for a averages of multiple exons. Since the data indicated a high level of variety of sequence features against the percent of exons exhibiting variability in binding site placement between individual exons, we that score that also exhibit PTBP1 dependent exon repression wanted to visualize target exons relative to each other. To display (Figure S5). These plots produced distinct sigmoidal curves where binding signals for individual exons we created heat maps of the most exons regulated by PTBP1 were found above or below a binding scores upstream, within, and downstream of each exon in particular score. This strongly suggests that a logistic regression the PTBP1 target set (Figure 3B). This display makes clear that the model incorporating each of these scores will be predictive of location of PTBP1 binding sites within its known target exons is PTBP1 repression. variable. We found that 60% of PTBP1 repressed exons are We developed a multinomial logistic regression model and predicted to have strong binding sites within the upstream intron. trained it on three classes of regulated exons (Figure 4A) [38]. The Most of these exons also have strong binding sites within either the training set included PTBP1 repressed exons, PTBP1 enhanced exon or the downstream intron, patterns that were observed exons, and non-regulated exons. Each exon in each class was previously [13,19,27]. However, other patterns of binding site scored for the four features found to correlate with PTBP1 placement are also seen, suggesting PTBP1 dependent exons are regulation (x through x ), including the 39 splice site strength, and 1 4 following multiple rules. Some repressed exons score highly for the PTBP1 binding scores for each of three regions: the 250 PTBP1 binding only within the exon or in both the exon and the nucleotides upstream of the exon, the exon itself, and the 100 downstream intron. About half of PTBP1 enhanced exons have nucleotides downstream of the exon. These intron lengths strong PTBP1 binding sites downstream (Figure 3B). These can encompass the regions of binding site enrichment for PTBP1 co-occur with upstream intron-binding sites, but rarely with exon dependent exons (Figure 3). binding sites. Interestingly, there are exons enhanced by PTBP1 The PTBP1-enhanced exons are fewer in number and show with strong upstream binding in the absence of other sites. These more limited enrichment of PTBP1 binding sites than PTBP1- data demonstrate the heterogeneity in the position of PTBP1 repressed exons making the prediction for these exons less accurate. binding sites for its target exons. This heterogeneity needs to be We first tested models that considered just PTBP1-repressed exons considered for predicting PTBP1 dependent regulation. relative to control groups. However, we found that including the PTBP1 repressed exons exhibited significantly higher average enhanced exons as a separate training group improved the binding scores in both the upstream intron and in the exon itself, prediction of repressed exons, even though enhanced exons than either the control group of alternative exons or the PTBP1 themselves are not as easily identified (data not shown). enhanced exons (Figure 3C). The average binding scores in the The trained model yielded values for the b coefficients that downstream introns were higher for both the PTBP1-repressed weight the different features contributing to the regulation. As and PTBP1-enhanced exons than the control group (Figure 3C), expected the upstream binding score was weighted most heavily in although not at the same statistical significance. The variability of predicting PTBP1 repression (Table S1), although binding scores in binding site placement within the smaller group of PTBP1- all three regions contributed to the score for PTBP1 repression. In enhanced exons presumably contributes to the weaker statistical contrast, we found that only the downstream binding score was correlation of binding scores with positive regulation. significantly associated with PTBP1 enhancement. The upstream We also compared the three exon sets for other features that score generated a b coefficient close to zero making it essentially might contribute to their ability to be regulated by PTBP1, neutral in the prediction of enhanced exons. The exon binding score including exon length, flanking intron length, and 59 and 39 splice was subject to a negative b coefficient, indicating that exon binding site strength. Most of these features were not statistically different reduces the probability of PTBP1 enhancement. Using these b among the three-exon groups. However, both PTBP1 enhanced coefficients, the trained models for repression or enhancement each and PTBP1 repressed exons were found to carry significantly yield a value of the g-function (logit) for an exon (x) given by the log weaker 39 splice sites than the control exon set, as measured by the of the ratio of the probability of repression or enhancement over the Analyzer Splice Tool (Figure 3C) [36,37]. probability that the exon is not regulated. From this, the probability These results indicate that PTBP1-repressed exons, and perhaps that an exon is repressed by PTBP1 can be determined from the two PTBP1-enhanced exons, exhibit an ensemble of sequence features g-values as shown in Figure 4A. that define them as PTBP1 regulated and that should allow their We assessed the multinomial logistic regression model by identification by sequence alone. recursively retraining on exon sets with one exon left out and then scoring the missing exon. This leave-one-out cross validation Prediction of PTBP1 repressed exons enabled assessment of the overall performance of the model [38] Alternative exons are generally regulated by multiple factors (Figure S6). The PTBP1 dependent exon repression logit showed that act both positively and negatively on their ability to be spliced. good prediction, with an area under the curve (AUC) value of Thus, an exon controlled by a regulator in one context might not 0.72, substantially greater than random guessing (AUC = 0.5). As be affected by it under other conditions where counteracting expected, the enhanced exon logit was not as accurate as the factors are present, or required cofactors are absent. This means repression logit (AUC = 0.57), although it was better than random that the most accurate predictions of splicing regulation will need (Figure S6A). Using these data, we assessed the sensitivity and PLOS Computational Biology | www.ploscompbiol.org 7 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 4. Scheme of the PTBP1 splicing regulation model and its application to an exon in Ptbp3.A. The PTBP1 splicing regulation model was trained on known PTBP1-regulated and non-regulated exons and used to predict new PTBP1-dependent exons. Prediction results were compared to changes in exon inclusion (PSI) measured by RT-PCR and RNA-seq. An exon from Ptbp3 is presented as a prediction example. From intron and exon sequences, PTBP1 binding scores and 39 splice site strength were calculated and fed into the regulation model. B. The model predicts exon 2 of Ptbp3 as repressed by PTBP1 with high probability (0.89). Ptbp1 knockdown in mouse neuroblastoma cells (N2A) confirmed de-repression of the exon (from PSI = 45 to PSI = 70). doi:10.1371/journal.pcbi.1003442.g004 specificity across the range of scores to define a decision threshold numbers of repressed exons. We sought to choose a threshold that for exon repression scores (Figure S6B). Increasing the threshold gave a low false positive rate over one that yielded more regulated increases the specificity by eliminating many false positives, but exons. We found that above a threshold score of 0.65 the false decreases the sensitivity of the model in identifying maximum positive rate was 10% or lower (Figure S6B). PLOS Computational Biology | www.ploscompbiol.org 8 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Applying the model to 4494 alternative cassette exons from (2.3%) with probability scores above 0.65 for being repressed by UCSC genome browser database, we found 243 exons (5.4%) that PTBP1. Among other activities, these exons were enriched in yielded a PTBP1 repression probability score greater than 0.65 and genes that function in calcium ion transport, cytoskeletal which were not in the training set. The 50 top-scoring cassette exons organization, intracellular transport, and synaptic transmission, are listed in Table 1. These included two exons that were reported all functions affected by previously known PTB targets (Table S2). previously to be PTBP1 targets. An exon of Gabrg2 yields a To assess splicing of this large set of predicted PTB targets, we probability score of 0.92. Although we could not confirm its used RNA-seq to generate a large dataset of exons that change after repression in N2A cells because of low expression of the transcript, Ptbp1 knockdown. RNA from control and PTBP1-depleted N2A the orthologous exon in rat is a well-characterized PTBP1 cells was subjected to high density short read sequencing on the repression target [39]. Exon 2 of Ptbp3 (Rod1), another known Illumina HiSeq platform using a strand specific, paired end protocol PTBP1 target [40], yielded a repression probability score of 0.89 [41]. Exons whose inclusion changed between the two samples were and was confirmed by RT/PCR to show increased inclusion after identified by alignment to an exon database and quantification of Ptbp1 knockdown (Figure 4B). We performed additional RT-PCR exon inclusion using the SpliceTrap program [42]. After filtering for validation in triplicate on a series of high and low scoring exons from read coverage and removing the training set, we identified 573 transcripts expressed in N2A cells (Figure 5 & Figures S7, S8 and alternative exons whose splicing was assayable in N2A cells. These S9). Seven of ten exons scoring above 0.65 were de-repressed after exons exhibit changes in percent exon inclusion (delta PSI) ranging Ptbp1 knockdown in N2A cells, yielding a validation rate of 70%. from 229% to 62% upon PTBP1 depletion. The exons were The actual false positive rate is difficult to estimate because exons binned by their PTBP1 repression probability scores and plotted for with high repression scores that are not affected by Ptbp1 depletion their change in PSI (Figure 6). The average changes in splicing were in N2A cells might be regulated by PTBP1 in other cells. An significantly correlated with the repression probability. Exons indication that this might be occurring is that the average inclusion scoring below 0.5 distributed around zero change in PSI, but above level (or percent spliced in value, PSI) of the putative false positives is this score the average exon inclusion is altered by PTBP1 depletion. significantly higher than the confirmed true positives in N2A cells, Most notably, exons with a repression probability score above 0.65 indicating that they will be less prone to change upon Ptbp1 exhibited significantly larger changes in splicing than exons with depletion and be more difficult to validate (Figure S8B). Thus, the lower scores. Exons with intermediate scores and hence weaker true positive rate may be greater than 70%. Importantly, the high binding sites show smaller changes in splicing than high scoring validation rate for exons scoring above 0.65 indicates that the exons. Setting a threshold of a 5% change in PSI as validation, 22 of binding model and the regulation model based upon it can identify 33 exons (67%) that scored above 0.65 for PTBP1 regulation were many new PTBP1 targets that were not previously known (Table1). confirmed as PTBP1 repression targets in N2A cells. At least some High scoring exons might also fail to be validated because of of the other 11 exons are presumably PTBP1 targets in other cells. regulation by other proteins. Knockdown of Ptbp1 induces To test the model in another cell type, we examined exons expression of its close homolog Ptbp2, which targets some of the reported to change after Ptbp1 knockdown in mouse C2C12 same exons [20] (Figure S7). To test whether PTBP2 was also myoblasts, as measured on splicing sensitive microarrays [43]. Very targeting the predicted PTBP1 repressed exons, we knocked down similar to what was observed in N2A cells, we found that exons with Ptbp2 or both Ptbp1 and Ptbp2 expression in N2A cells and re- high repression probabilities showed significant de-repression upon assayed the exons in triplicate (Figures S10, S11 & S8A). Although the Ptbp1 knockdown compared to exons with low repression some exons showed greater inclusion in the double knockdown probabilities (Figure S12). Of 29 exons assayed on the arrays with a compared to depletion of Ptbp1 alone, this did not validate any repression probability above 0.65, 19 exons were confirmed as additional predicted PTBP1 repressed exons. We did identify some PTBP1 repressed on the array (q-value,0.05), yielding a validation high and low scoring exons showing more complex regulation by rate of 66%. Thus the model performed very similarly in C2C12 the two PTB proteins (Figure S10 & S11). and N2A cells. Among the 11 high scoring exons identified as We also examined a set of low scoring exons (probability score# unchanged after PTBP1 knockdown in N2A cells only 3 were 0.2) by RT-PCR after Ptbp1 and/or Ptbp2 depletion (Figure 5B assayed on the array and expressed in C2C12 cells. These again and Figure S11). All of these exons (8 of 8) failed to respond to the showed high inclusion in C2C12 prior to knockdown and so were loss of PTBP1 and are likely true negatives. Thus, PTBP1 difficult to assay for derepression. Thus, it is difficult to use the repression scores above 0.65 and below 0.2 were highly predictive C2C12 data to draw conclusions about the false positive rate. for regulation and its absence, respectively. As expected, interme- The logistical model gives us a new tool for studying the diate scores were less consistent in their predictive value (Figure regulation of alternative splicing. Using it, we can now scan genomic S9). Some exons in the intermediate scoring group were affected sequence to score exons for PTBP1 regulation. Applying the model by PTB proteins and will be interesting to assess further. genomewide, the PTBP1 repression probability scores were The prediction of PTBP1-repressed exons was improved by integrated into the UCSC genome browser. These data, displayed treating PTBP1-enhanced exons as a separate class, but the with the RNAseq data from N2A cells are available at our website probability scores for PTBP1 enhancement did not consistently (http://www.mimg.ucla.edu/faculty/black/ptbatweb/). A novel identify new PTBP1 target exons (data not shown). This is likely in PTBP1 repressed exon in the Kcnq2 gene is shown in Figure 6B. part due to the smaller number of exons in the training set and The logistic model thus allows the assessment of any exon across the their heterogeneity, with some possibly being indirect targets. transcriptome for likely PTBP1 regulation. These predictions will likely improve with training on larger numbers of PTBP1 enhanced exons as they are identified. Discussion However, it is possible that simply the presence of the PTBP1 binding site is not sufficient for predicting PTBP1 enhancement New features of PTBP1 binding sites and that binding sites for other factors will need to be considered. We have developed two computational models, one that allows We next tested the model on a genomewide scale, by applying it accurate prediction of PTBP1 binding sites and another that to a set of 168,111 mouse internal exons and ranking them by their predicts likelihood of PTBP1 repression of exons across the probability of PTBP1 repression. This analysis yielded 3824 exons transcriptome. These models uncovered several new features of PLOS Computational Biology | www.ploscompbiol.org 9 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Table 1. PTBP1 repressed exons identified by the splicing model. PTB Binding Scores Gene Name Gene Description mm9 coordinates 39 p(Repressed) Upstream Downstream Splice site Intron Intron (250 nt) Exon (100 nt) Strength Pax6 paired box gene 6 chr2:105523985–105524115(+) 8.35 21.49 20.80 0.27 0.99 Mbd5 methyl-CpG binding chr2:49134101–49135303(+) 6.46 20.93 2.27 20.32 0.98 domain protein 5 Arhgap24 Rho GTPase activating chr5:102981145–102981338(+) 6.47 0.10 20.05 0.46 0.97 protein 24 Tle1 transducin-like enhancer chr4:71819247–71819451(2) 4.71 0.05 21.26 22.56 0.94 of split 1 Acsl6 acyl-CoA synthetase chr11:54150438–54150515(+) 4.16 1.40 20.03 20.82 0.94 long-chain family Ryr1 ryanodine receptor 1, chr7:29829938–29829955(2) 4.78 20.08 20.88 21.71 0.94 skeletal muscle Ankhd1 ankyrin repeat and KH chr18:36784163–36784921(+) 4.37 0.03 1.44 20.64 0.93 domain containing 1 Slc39a14 solute carrier family 39 chr14:70713408–70713577(2) 3.51 1.16 21.27 23.20 0.92 (zinc transporter) Gabrg2 gamma-aminobutyric chr11:41727472–41727495(2) 1.95 2.77 0.56 23.79 0.92 acid (GABA) A receptor Itga7 integrin alpha 7 chr10:128378878–128378997(+) 4.14 0.29 1.13 20.25 0.92 Iqsec2 IQ motif and Sec7 chrX:148615540–148615635(+) 4.88 0.68 20.13 0.71 0.91 domain 2 Smarca2 SWI/SNF related, chr19:26825612–26825646(+) 3.94 20.09 1.23 20.86 0.91 matrix associated, actin dependent regulator of chromatin Zfand3 zinc finger, AN1-type chr17:30197755–30197795(+) 4.17 2.34 0.90 1.80 0.91 domain 3 Agap2 ArfGAP with GTPase chr10:126527198–126527257(+) 3.57 0.06 20.53 23.08 0.90 domain, ankyrin repeat and PH domain 2 Ttn Titin chr2:76723554–76723832(2) 2.93 1.06 1.19 21.83 0.90 Ptbp3 ROD1 regulator chr4:59559021–59559054(2) 3.80 0.57 1.66 0.73 0.89 of differentiation 1 (S. pombe) Mapk8 mitogen-activated chr14:34203859–34203930(2) 2.35 1.17 1.01 23.48 0.89 protein kinase 8 Snap91 synaptosomal-associated chr9:86693373–86693534(2) 2.60 1.89 20.35 22.17 0.88 protein 91 Fmnl1 formin-like 1 chr11:103059449–103059547(+) 3.93 20.60 21.36 22.70 0.88 Phldb1 pleckstrin homology-like chr9:44509029–44509169(2) 3.20 0.57 1.05 20.66 0.87 domain, family B 2310035C23Rik RIKEN cDNA 2310035C23 chr1:107637012–107637094(+) 2.03 1.76 0.85 22.46 0.87 gene Arnt aryl hydrocarbon receptor chr3:95270715–95270759(+) 3.48 20.36 2.53 0.09 0.87 nuclear translocator Smyd2 SET and MYND domain chr1:191723697–191723807(2) 3.33 0.33 20.28 21.64 0.86 containing 2 Ap2a1 adaptor protein complex chr7:52158832–52158897(2) 3.35 20.12 20.89 22.69 0.86 AP-2, alpha 1 subunit Klra killer cell lectin-like chr6:130329011–130329100(2) 2.82 3.18 0.62 0.77 0.86 receptor, subfamily A Spag9 sperm associated chr11: 93942054–93942068(+) 0.99 3.01 1.62 23.03 0.86 antigen 9 Col4a3bp collagen, type IV, chr13:97386949–97387026(+) 2.81 1.23 0.74 20.93 0.86 alpha 3 binding protein Garnl3 GTPase activating chr2:32941395–32941464(2) 4.15 0.38 0.36 0.94 0.86 RANGAP domain-like 3 PLOS Computational Biology | www.ploscompbiol.org 10 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Table 1. Cont. PTB Binding Scores Gene Name Gene Description mm9 coordinates 39 p(Repressed) Upstream Downstream Splice site Intron Intron (250 nt) Exon (100 nt) Strength Dennd1a DENN/MADD domain chr2:37982049–37982168(2) 3.37 0.80 1.35 0.59 0.86 containing 1A Ms4a7 membrane-spanning chr19:11400297–11400353(2) 2.79 2.35 0.37 20.05 0.86 4-domains, subfamily A BC030307 cDNA sequence BC030307 chr10:86169981–86170089(+) 2.75 20.16 1.95 22.40 0.85 Phactr1 phosphatase and actin chr13:43154940–43155146(+) 2.73 1.25 0.50 20.97 0.85 regulator 1 R3hdm2 R3H domain containing 2 chr10:126902187–126902240(+) 1.66 1.98 1.96 21.57 0.84 Cdc14b CDC14 cell division chr13:64306579–64306725(2) 1.42 2.75 2.44 20.58 0.84 cycle 14B Ubqln1 ubiquilin 1 chr13:58282183–58282266(2) 2.88 0.98 20.06 21.17 0.84 Ttn Titin chr2:76739898–76740179(2) 2.63 20.07 1.49 22.38 0.84 Stx3 syntaxin 3 chr19:11857290–11857400(2) 3.00 21.12 2.26 23.62 0.84 Slc8a3 solute (sodium/calcium) chr12: 82310340–82310458(2) 1.84 1.25 1.76 22.25 0.84 carrier family 8 Zfp62 zinc finger protein 62 chr11:49028057–49028156(+) 3.27 1.98 20.51 0.19 0.83 Dlg1 discs, large homolog 1 chr16:31771843–31771941(+) 1.53 1.98 1.87 21.65 0.83 (Drosophila) Nrxn2 neurexin II chr19:6463824–6463847(+) 3.35 21.37 1.33 22.26 0.83 Klra7 killer cell lectin-like chr6:130179953–130180042(2) 2.68 2.11 20.39 20.63 0.83 receptor, subfamily A Picalm phosphatidylinositol chr7:97330729–97330878(+) 1.15 2.15 2.30 22.37 0.83 binding clathrin assembly Acad8 acyl-Coenzyme A chr9:26798168–26798277(2) 2.61 0.88 20.31 21.86 0.83 dehydrogenase family Epn1 epsin 1 chr7:5033620–5033723(+) 3.92 0.65 0.07 1.06 0.82 Grip1 glutamate receptor chr10:119422530–119422685(+) 2.66 20.74 2.61 23.13 0.82 interacting protein 1 Csmd3 CUB and Sushi multiple chr15:47587514–47587627(2) 2.42 2.16 0.20 20.45 0.82 domains 3 Lrrfip1 leucine rich repeat (in FLII) chr1:92990137–92990214(+) 2.02 0.40 3.44 22.21 0.82 interacting protein 1 Srsf11 serine/arginine-rich chr3:157703405–157703586(+) 1.09 2.05 20.56 24.75 0.82 splicing factor 11 Tmem209 transmembrane chr6:30441087–30441184(2) 3.82 0.16 0.10 0.64 0.82 protein 209 The 50 highest scoring exons predicted to be repressed by PTBP1 based on sequence alone. doi:10.1371/journal.pcbi.1003442.t001 RNA recognition by PTBP1 and the properties of its target exons. The base-specific contacts that PTBP1 makes with Guanosine The PTBP1 binding model was based on triplets following the are not yet clear. Recent studies of RNA recognition by SRSF2 structures of the PTBP1 RRM domains, whose sequence specific (SC35) protein have shown that the element GGAG can be contacts are each primarily to three nucleotides. We find that the recognized by the same RRM as CCAG by flipping the initial two set of triplets that increase the probability of binding includes the G nucleotides to the syn conformation [45]. It will be very expected pyrimidine motifs, particularly those with alternating interesting to investigate whether a similar anti to syn switch cytosines and uridines. However, many triplets with guanosine occurs in RNA bound by PTBP1, when C residues are replaced residues also increase binding probability. In contrast, adenosine with G. residues have a negative effect on binding. Thus, RNA recognition Previous characterizations of PTBP1 binding sites have focused by PTBP1 is not solely dependent on pyrimidine nucleotides. The on finding enriched short motifs within populations of bound recognition of G residues by PTB was unexpected, although some RNAs or regulated exon sequences [13,44,46,47,48]. These previously characterized PTB binding sites did contain G residues methods generally identify elements whose short length will allow [13,44]. With this model, we can now predict PTBP1 binding interaction with only one RRM domain. Searching for new affinity to any site in the transcriptome. binding sites comprised of clusters of these short elements can PLOS Computational Biology | www.ploscompbiol.org 11 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets PLOS Computational Biology | www.ploscompbiol.org 12 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 5. Validation of novel PTBP1-repressed exons by RT-PCR. A. Candidate PTBP1-repressed exons with probability greater than 0.65 were validated by RT-PCR following Ptbp1 knockdown. Data shown are averages 6 standard error of PSI (Percent of Spliced In) from biological triplicates. Statistical analysis was performed using paired one-tailed Student’s t-test (p-values,0.01**, ,0.05*). B. Exons with low PTBP1 repression probabilities (#0.2) were also validated by RT-PCR following Ptbp1 knockdown in biological triplicates. doi:10.1371/journal.pcbi.1003442.g005 identify higher affinity sites but does not consider all elements or to achieve accurate predictions from the model. However, treating rank them. Crosslinking-immunoprecipitation experiments allow these as a separate exon class improves the prediction of PTBP1 large numbers of binding regions to be identified. However, not all repression. We find that for probability scores above 0.65 the the sequence within a CLIP tag will be contacting the protein and model is strongly predictive of PTBP1 repression. Applying this it is difficult to relate CLIP signals to binding affinity. The HMM criterion across the transcriptome, we identified hundreds of new allowed the individual assessment of different short elements PTBP1 target exons. within the CLIP clusters, showing that they segregated into two Alternative exons are generally regulated by multiple proteins states. The ranking of the triplets for their contributions to one of acting in combination, and a particular exon will often be subject these states yielded a model where complex clusters of short to both positive and negative regulation by antagonistic factors. elements could be assessed for binding and yielded accurate For a model based on one factor, these other proteins will predictions of binding affinity. Many RNA binding proteins are confound predictions. Exons with high PTBP1 binding scores may similar to PTBP1 in having multiple domains that may each make be counteracted by antagonistic factors in some cell types. different base specific contacts with RNA. The widespread Alternatively, synergistic factors may allow an exon with a generation of CLIP-seq datasets will allow the modeling of RNA relatively weak binding site to still recruit PTBP1. Thus, a model recognition by almost any protein based on a large number of based on one factor will be limited in its predictive power. In this known binding sites. study, our intent was to measure the effect of PTBP1 binding alone Using the same modeling approach, we also developed a before considering the contributions of other factors. The logistic binding model for PTBP2 (neuronal PTB) using a published modeling allowed the contributions of different binding site PTBP2 CLIP dataset [49]. PTBP2 is about 70% identical to placements to PTBP1 regulation to be measured. PTBP1 in sequence, and has only two amino acid changes among Several studies have used Bayesian models to dissect the the residues making direct contact with RNA [17]. We found that regulatory properties of exons [7,54]. These models can generate the binding models for two PTB proteins were also nearly identical accurate predictions by incorporating a wide variety of sequence, indicating that the two proteins are likely to differ more in their expression and conservation data. However, because so many protein/protein interactions than in their RNA binding sites (Data disparate variables are incorporated, it can be difficult to draw not shown). mechanistic conclusions from these models regarding any one protein. For example, the presence of high pyrimidine density upstream from the branch point can be predictive of exons Defining PTBP1 target exons showing neuronal specific inclusion [7,55]. This is presumably in Several PTBP1 target exons have been analyzed in detail part due to many neuronal exons being regulated by PTBP1 and [17,50]. These exons vary in the placement and action of their PTBP2. However, a subset of these exons may be regulated by PTBP1 binding sites. It is common for PTBP1-repressed exons to other factors with pyrimidine rich binding sites. In the long term, it have a binding site upstream, often encompassing the branch will be most accurate to develop predictive binding models for point of the 39 splice site [39]. Exons can also be repressed by each protein, similar to the PTBP1 model here, and then to PTBP1 binding within the exon [19,32,33]. Other exons contain incorporate each of these binding models into a larger network downstream binding sites that are needed in conjunction with an model. Such an approach will allow the analysis of the many upstream site to achieve splicing repression [51,52,53]. Although overlapping regulatory programs controlled by RNA binding acting as a repressor for most of its targets, PTBP1 also activates proteins. the splicing of a group of exons. There have been divergent reports about placement of PTBP1 binding sites needed to mediate PTBP1 enhancement of splicing. The PTBP1 binding model Materials and Methods allowed us to examine PTBP1 binding site placement across a Hidden Markov Model for PTBP1 binding affinity large set of known PTBP1 target exons. Nearly all exons had prediction predicted high affinity PTBP1 binding sites nearby. We found that A Hidden Markov Model (HMM) was designed and trained by more than half of PTBP1 repressed exons have high affinity binding sites upstream, and a fraction of PTBP1 enhanced exons an expectation–maximization (EM) method (Baum-Welch algo- have high affinity sites downstream. These exons fit with recent rithm) using published PTBP1 CLIP data [13,29,30]. In total, results on several other splicing regulators where the placement of 48,604 PTBP1-CLIP cluster sequences were used to train model the binding site determines the direction of the regulatory effect parameters. During the training step, multiple initial values were [12,14,34]. However, for PTBP1 these rules are not so clear. Some tested to avoid a local maximum problem. Trained parameters PTBP1 repressed exons have their strongest predicted binding site included emission probabilities for nucleotide triplets, initial downstream or within the exon. These results indicate that there probabilities and transition probabilities between states [29,30]. are fundamental differences between the mechanisms of PTBP1 The trained model was used to score RNA sequences. The raw mediated splicing regulation, and those governing regulation by PTBP1 binding score is defined as a log-odds ratio that compares certain other splicing factors. the score of a sequence from the HMM over the score from a To quantify the predictive value of the PTBP1 binding scores background model. Since CLIP experiments do not have an inherent corresponding negative dataset, we generated computa- for PTBP1 repression, we built a logistic model for PTBP1 regulation. For exons repressed by PTBP1, binding scores for the tional negative datasets and tested different background models upstream, downstream and exon sequences all contribute to the (Figure S4). We found that a background model that values all probability of repression. Exons enhanced by PTBP1 were too few triplets equally yielded the most accurate binding scores [29]. Raw PLOS Computational Biology | www.ploscompbiol.org 13 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 6. Large-scale validation of novel PTBP1-repressed exons by RNA-seq. A. Validation of the PTBP1 splicing model using RNA-seq. After Ptbp1 knockdown, we performed RNA-seq experiments and estimated changes in PSI (Percent of Spliced In) for 573 cassette exons. The graph shows average delta PSI values for exons, grouped by their probabilities to be repressed by PTBP1. The number of exons in the corresponding probability bin is given by n. P-values were calculated from one-tailed Student’s t-test. B. A genome browser screenshot of a novel PTBP1-regulated exon: exon 2 of the Kcnq2 gene. For whole internal mouse exons, we created custom genome browser tracks to visualize the PTBP1 splicing model and mapped RNA seq reads. doi:10.1371/journal.pcbi.1003442.g006 PLOS Computational Biology | www.ploscompbiol.org 14 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets scores were further normalized and converted to z-scores. For the with canonical splice sites (GU-AG). An exon was classified as 69 mer RNA sequences used in binding assays, scores were PTBP1 repressed or enhanced when 1) the inclusion level (PSI) of normalized by 100,000 random sequences with same length its minor isoform was greater than 5% in both the control and (Figure S3). This yielded very accurate predictions of binding knock-down samples and 2) the inclusion level of its minor isoform affinity (Figure 2). was changed by 30% or more in the Ptbp1 knock down condition When considering binding scores in genomic sequence, exons compared to the control sample. Next, we collected sequence features for each exon and its flanking exons. The features and upstream or downstream intron regions have different base compositions and will yield different average binding scores. Thus, included PTBP1 binding scores, 59 and 39 splice site strengths, exon/intron lengths, and word frequencies. The PTBP1 binding to score binding sites adjacent to possible regulated exons, it is scores were calculated from the PTBP1 binding model described more informative to score sites relative to equivalent sequence above. The strength of splice sites was calculated by the splice-site regions. From the annotated mouse genome, we retrieved 168,111 analyzer tool [37]. Using a mouse whole internal exon set, we internal exons and their flanking introns as separate sequence sets normalized features and fed them into the model. The PTBP1 using a python library, Pygr. We scored log odds of these sequences splicing model is based on a multinomial logistic regression with the trained model. Since the lengths and base compositions of framework using the following steps: 1) selection of initial variables intronic and exonic sequences are different, and binding scores with a moderate level of association (p-value from t-test,0.25), 2) automatically increase with length (Figure S13) [29], we grouped removal of outlier exons, 3) stepwise variable selection [38]. We sequences by their location and sequences in each group were scored mouse internal exons with the trained PTBP1 splicing sorted according to length into bins of 1000 sequences each. The model and validated candidate exons with RT-PCR and RNA-seq average score and standard deviation were determined for each experiments. Exons from the training set were excluded from the bin. These values were used to transform the raw scores into validation. z-scores for each upstream intron, downstream intron, and exon sequence. We localized the PTBP1 binding sites along each RNA sequence using the Viterbi algorithm [29,30]. Validation of exon candidates by RT-PCR and RNA-seq To test alternative splicing events for candidate exons, we assayed exon inclusion levels in cells following Ptbp1, Ptbp2, and Validation of PTBP1 binding model scores by binding double Ptbp1 & Ptbp2 knock down. The knockdown experiment assay was performed as described previously with minor modification To test predicted PTBP1 binding scores, we selected thirteen [20]. Mouse neuroblastoma (N2A) cells were cultured in DMEM mouse exon/intron RNA sequences (69 nucleotides) exhibiting a with 10% FBS and 2 mM L-glutamine. At 70 to 80% confluency, range scores. In the selection, other sequence features such as cells were trypsinized and suspended in the growth medium. secondary structure were not considered. Target RNAs were DNA–Lipofectamine 2k (Invitrogen) complexes were prepared transcribed in vitro from dsDNA using T7 RNA polymerase and and mixed with cells in a tube according to manufacturer’s subjected to an electrophoretic mobility shift assay (EMSA). instructions. Tubes were incubated for 5 h with mixing every half During the transcription, radioactive a-32P UTP was incorporat- hour. Then cells were centrifuged and cultured in plates for 3 d. ed into RNA to visualize the probes. The RNA probes were then Proteins and RNA was extracted from collected cells. Protein denatured for 2 min at 85uC and cooled down on ice immediately samples were subjected to fluorescence immunoblotting to to reduce secondary structure formation. Binding assays were monitor knockdown efficiency of Ptbp1 and Ptbp2. Total RNA carried out as previously described with some modifications [27]. was collected using Trizol (Invitrogen) according to the manufac- Specifically, each gel mobility shift reaction (10 mL) contained the turer’s instructions. The RNA was further treated with DNase I to indicated amounts of recombinant human PTBP1 in 6 mLDG avoid DNA contamination. For RT-PCR (Reverse Transcription- buffer (20 mM Hepes-KOH ph 7.9, 20% glycerol, 80 mM PCR) assays, the RNA was reverse transcribed to cDNA with potassium glutamate, 0.2 mM EDTA, 0.2 mM PMSF), 1 mL random hexamers using SuperScript enzyme (Invitrogen) follow- 22 mM MgCl2, 1 mL 0.5 mg/ml tRNA, 0.5 mL RNase inhibitor ing the manufacturer’s instructions. PCR reactions were per- (20 unit, RNaseOut from invitrogen), 0.5 mL DEPC treated H2O, formed to assay alternative splicing of particular target exons. and 1 mL 100 nM RNA probe. At first, all reaction components First, forward and reverse PCR primers were designed for the excluding RNase inhibitor, tRNA, and RNA probes were mixed flanking exons using PRIMER3 program [56]. To label PCR and incubated for 8 min at 30uC. Then RNase inhibitor and products, a 59 fluorescent-labeled universal primer (59-FAM- tRNA were added and mixed. RNA probe was then added and CGTCGCCGTCCAGCTCGACCAG-39) was added to the PCR the reaction was incubated for an additional 15 min. The reaction and a universal priming site was introduced to the 59 end reactions were put on ice for 5 min and mixed with 1.2 mL of the forward primer (59-CGTCGCCGTCCAGCTCGACCAG- glycerol loading dye (30% glycerol). They were separated on 8% Forward Primer-39). Each PCR reaction (15 mL) was carried out native polyacrylamide gels with 25 mM Tris-Gly running buffer in with 1.5 picomole of the forward primer and 6.75 picomole of the a cold room. Gels were dried and exposed to a phosphor screen. reverse and universal primers [57]. PCR amplification proceeded Then images were scanned using Typhoon 9410 and quantified with an initial denaturation at 94uC for 4 m followed by 24 cycles using ImageQuant TL program (GE Lifesciences). The apparent of 94uC for 30 s, at a melting temperature of the reverse primer for Kd values were estimated by fitting the data to non-linear curves 45 s, and 72uC for 45 s, with a final extension step at 72uC for using Prism software. 10 m. The samples were mixed with 26 formamide buffer (Formamide with 1 mM EDTA pH 8.0) and denatured at 95uC Logistic regression model for PTBP1 dependent exon for 5 min. Then samples were chilled on ice and run on 8% prediction denaturing polyacrylamide gels. Gels were directly scanned by An exon training set was compiled from previous microarray Typhoon and quantified by ImageQuant program. and RT-PCR experiments [20,35]. The training set was composed RNA-seq libraries were constructed following standard proto- with 68 PTBP1 repressed, 37 PTBP1 enhanced, and 69 non- cols (Illumina TruSeq RNA Sample Prep Kit). To make strand- PTBP1 regulated simple cassette exons. We only considered exons specific libraries, we added two extra steps to the protocol [41]. PLOS Computational Biology | www.ploscompbiol.org 15 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets After first strand cDNA synthesis, remaining dNTPs were slightly improves the linear fit (20.91 vs. 20.95) for some strong removed by a size selection on beads (AMPure XP). Second- binders. However, it wrongly predicted some binders as non- strand cDNA was synthesized with a dNTP mix containing dUTP binders, which reduced the rank correlation (20.95 to 20.90). instead of dTTP. The reaction contained samples eluted in 50 ml The two shuffled models did not perform well. resuspension buffer, 2 ml56 FS buffer, 1 ml 50 mM MgCl2, 1 ml (TIF) 100 mM DTT, 2 ml 10 mM dUTP nucleotides mix, 15 ml Second Figure S5 Correlations of particular sequence features Strand Buffer (Invitrogen), 0.5 ml E.coli DNA Ligase (10 U/ with PTBP1 repression. For PTBP1-repressed and PTBP1 ml;NEB), 0.5 ml RNase H (2 U/ml;Invitrogen), 2 ml DNA E.coli non-regulated exons, we calculated scores for sequence features Polymerase I (10 U/ml;NEB). The reaction was incubated for 2 h and determined the fraction of PTBP1-repressed exons in each at 16uC. After sequencing adaptors were ligated, 1 ml USER score bin. Shown are the graphs for the 39 splice site score, and the (Uracil-Specific Excision Reagent enzyme; NEB) was added to PTBP1 binding scores in the upstream intron, the exon, and the reactions to degrade the second strand cDNA. The samples were downstream intron plotted against the percent of exons within the incubated for 15 min at 37uC and the reaction were inactivated at score bin that are PTBP1-repressed. 94uC for 5 min. The samples were put in ice and then subjected to (TIF) PCR amplification. Average size of inserts was about 225 bp and the libraries were subjected to 100 bp paired-end sequencing Figure S6 Performance of PTBP1-dependent splicing (Illumina HiSeq2000 platform). Using SpliceTrap [42], 60–65% models. A. Receiver Operating Characteristic Curves in a leave- of reads were mapped to exon duos or trios. In total, 180M one-out cross validation for each logit: exon repression (left) and (179,511,116) and 145M (145,334,711) paired end reads were exon enhancement (right). B. Sensitivity and specificity plotted used to infer exon inclusion ratios in the control and Ptbp1 across the whole threshold range. Sensitivity is defined as the knockdown conditions, respectively. The data have been deposited percent of true repressed exons that are correctly predicted as in NCBI’s Gene Expression Omnibus [58] and are accessible repressed at the corresponding threshold. Specificity is defined as through GEO Series accession number GSE45119. the percent of actual non-repressed exons that are correctly predicted as non-repressed at the corresponding threshold. Supporting Information (TIF) Figure S1 Five-fold cross-validation of the PTBP1 Figure S7 ShRNA mediated depletion of PTBP1 and binding model using Hela CLIP clusters. To test the two PTBP2. Duplicate immunoblots after shRNA knockdown of state model of binding and non-binding triplets, we divided the PTBP1, PTBP2 or both proteins. Note that depletion of PTBP1 CLIP-data to the five subsets. In each plot, four sub sets were used induces expression of PTBP2 as observed previously. Numbers in training the model and one subset was subjected to scoring. We above each lane indicate the fluorescence intensity for PTBP1 or then compared scores from the CLIP-subset sequences to random PTBP2 relative to the control lane. sequences picked from same genic regions that contained the (TIF) CLIP clusters. As shown, sequences from CLIP-subset generated Figure S8 Characteristics of false positive exons. A. significantly higher scores than random. The results indicate that Exon inclusion was measured for three false positives exons after the triplets identified by the HMM as predictive of state 1 are Ptbp1 knockdown (left), or Ptbp2 knockdown, and Ptbp1 & Ptbp2 predictive of PTBP1 CLIP sites and thus of protein binding. double knockdown (right). P-values were calculated from biolog- (TIF) ical triplicates using paired one-tailed t-tests. B. False positive Figure S2 The density of PTBP1 binding triplets exons exhibit higher PSI values prior to PTBP1 depletion. Box predicted by the HMM peaks at the aligned crosslink plot of exon inclusion for twenty-nine false positive exons showing sites from Human ESC iCLIP clusters. little change in splicing by RNA-seq after PTBP1 depletion (delta (TIF) PSI,5%) that score with high probability to be repressed (.0.55). Exon inclusion levels prior to PTBP1 depletion are compared to Figure S3 PTBP1 binding model scores and validation. thirty-six true positive exons. A. Summary statistics and the distribution of raw and normalized (TIF) PTBP1 binding scores for 100,000 random sequences. B. Electrophoretic mobility shift assay of RNAs with various PTBP1 Figure S9 Two exons with intermediate scores for binding scores. RNAs were transcribed in vitro, incubated with PTBP1 repression show complex responses to PTBP1 increasing concentrations of purified PTBP1 (0 to 200 nM), and and PTBP2 depletion. Exon inclusion was measured after the bound and unbound RNA separated on native gels. Arrows Ptbp1 depletion (left), or after Ptbp2 and Ptbp1/Ptbp2 double indicate RNA-protein complexes. The fraction of PTBP1-bound depletion (right). P-values were calculated from biological RNA is plotted below for each RNA. triplicates using paired one-tailed t-tests. (TIF) (TIF) Figure S4 Evaluation of different background models Figure S10 PTBP2 dependence of predicted PTBP1 for scoring PTBP1 binding. Four background models were target exons. RT-PCR of high probability PTBP1 exon targets evaluated. The uniform distribution model assumes equal following Ptbp2 knockdown or Ptbp2/Ptbp1 double knockdown. frequencies of triplets. The PTBP1 target gene set model used Relative band intensities of the gels in triplicate on the right are random sequences from genes containing PTBP1 CLIP clusters. plotted on the left to show the average delta PSI 6 SE (Percent The two shuffled models used shuffled CLIP cluster sequences Spliced In). P-values were calculated from paired one-tailed t-tests maintaining mono or di nucleotide ratios. We calculated PTBP1 with PSI values in control samples. binding scores based on each background model and compared (TIF) the scores to the measured dissociation constants. Based on the rank correlation, the uniform distribution model worked best. The Figure S11 PTBP2 dependence of predicted non-PTBP1 PTBP1 target gene set model showed comparable performance. It target exons. RT-PCR of exon with low probabilities for PTBP1 PLOS Computational Biology | www.ploscompbiol.org 16 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets repression (#0.2). Ptbp2 knockdown and Ptbp1/Ptbp2 double knock Table S1 Trained PTBP1 splicing regulation model. The table presents a summary of the multinomial logistic regression down with data analysis as in Figure S10. (TIF) model for PTBP1 splicing regulation, including estimated coefficients and their statistics. Figure S12 Boxplot of PSI (Percent of Spliced In) values (TIF) estimated from splicing sensitive microarray data for Table S2 Enriched gene ontology categories for novel exons expressed in mouse C2C12 myoblasts. Exons with PTBP1-repressed exons. The table lists ontology entries higher (.0.65) and lower (,0.2) repression probabilities are enriched in genes with predicted PTBP1-repressed exons (prob- compared. Exons used in the original training set were excluded ability score of exon repression .0.65). Whole mouse internal from the plot. The number of exons in the corresponding exons were used as the control set, and p-value cut off was 0.05. probability bin is given by n. The p-value was calculated from a Gene ontology analysis was performed using the GOTM web one-tailed Student’s t-test. server. (TIF) (TIF) Figure S13 Distribution of PTBP1 binding scores of exons and introns before and after normalization. Acknowledgments 168,111 exons and their flanking introns from the set of annotated We thank Manny Ares, Matteo Pellegrini, Yi Xing, and Christopher Lee mouse internal exons were subjected to scoring and normalization. for thoughtful suggestions on the experiments, the analysis and manuscript, Raw PTBP1 Binding scores are affected by sequence length and and all members of the Black laboratory for helpful discussion. We thank base composition. To account differences in these features Shalini Sharma and Batoul Amir-Ahmady for advice on the EMSA between introns and exons in normalized scores, we grouped the experiments, and Paul Boutz for advice on the PTBP1 knock down exons and their upstream and downstream introns separately. The experiments. We especially thank Shalini Sharma for the recombinant sequences in each group were sorted according to length into bins PTBP1 proteins. of 1,000 sequences each. The average scores and standard deviations were determined for each bin. These values were used Author Contributions to transform the raw scores into z-scores for each sequence per Conceived and designed the experiments: AH DLB. Performed the bin. experiments: AH PS AJL. Analyzed the data: AH DLB. Contributed (TIF) reagents/materials/analysis tools: YZ XDF. Wrote the paper: AH DLB. References 1. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, et al. (2008) Alternative 19. Llorian M, Schwartz S, Clark TA, Hollander D, Tan LY, et al. (2010) Position- isoform regulation in human tissue transcriptomes. Nature 456: 470–476. dependent alternative splicing activity revealed by global profiling of alternative 2. Black DL (2003) Mechanisms of alternative pre-messenger RNA splicing. splicing events regulated by PTB. Nature structural & molecular biology 17: Annual review of biochemistry 72: 291–336. 1114–1123. 20. Boutz PL, Stoilov P, Li Q, Lin CH, Chawla G, et al. (2007) A post- 3. Gabut M, Chaudhry S, Blencowe BJ (2008) SnapShot: The splicing regulatory transcriptional regulatory switch in polypyrimidine tract-binding proteins machinery. Cell 133: 192 e191. reprograms alternative splicing in developing neurons. Genes & development 4. Matlin AJ, Clark F, Smith CWJ (2005) Understanding alternative splicing: 21: 1636–1652. Towards a cellular code. Nature Reviews Molecular Cell Biology 6: 386–398. 5. Blencowe BJ (2006) Alternative splicing: New insights from global analyses. Cell 21. Boutz PL, Chawla G, Stoilov P, Black DL (2007) MicroRNAs regulate the 126: 37–47. expression of the alternative splicing factor nPTB during muscle development. Genes & development 21: 71–84. 6. Wang Z, Burge CB (2008) Splicing regulation: from a parts list of regulatory 22. Makeyev EV, Zhang J, Carrasco MA, Maniatis T (2007) The MicroRNA miR- elements to an integrated splicing code. RNA 14: 802–813. 124 promotes neuronal differentiation by triggering brain-specific alternative 7. Barash Y, Calarco JA, Gao W, Pan Q, Wang X, et al. (2010) Deciphering the pre-mRNA splicing. Mol Cell 27: 435–448. splicing code. Nature 465: 53–59. 8. Ule J, Jensen K, Mele A, Darnell RB (2005) CLIP: a method for identifying 23. Zheng S, Gray EE, Chawla G, Porse BT, O’Dell TJ, et al. (2012) PSD-95 is protein-RNA interaction sites in living cells. Methods 37: 376–386. post-transcriptionally repressed during early neural development by PTBP1 and PTBP2. Nature Neuroscience 15: 381–U202. 9. Darnell RB (2010) HITS-CLIP: panoramic views of protein-RNA regulation in 24. Tang ZZ, Sharma S, Zheng SK, Chawla G, Nikolic J, et al. (2011) Regulation of living cells. Wiley interdisciplinary reviews RNA 1: 266–286. the Mutually Exclusive Exons 8a and 8 in the CaV1.2 Calcium Channel 10. Hafner M, Landthaler M, Burger L, Khorshid M, Hausser J, et al. (2010) PAR- Transcript by Polypyrimidine Tract-binding Protein. Journal of Biological CliP–a method to identify transcriptome-wide the binding sites of RNA binding Chemistry 286: 10007–10016. proteins. Journal of visualized experiments : JoVE. 25. Oberstrass FC, Auweter SD, Erat M, Hargous Y, Henning A, et al. (2005) 11. Konig J, Zarnack K, Rot G, Curk T, Kayikci M, et al. (2010) iCLIP reveals the function of hnRNP particles in splicing at individual nucleotide resolution. Structure of PTB bound to RNA: specific binding and implications for splicing Nature structural & molecular biology 17: 909–915. regulation. Science 309: 2054–2057. 26. Lamichhane R, Daubner GM, Thomas-Crusells J, Auweter SD, Manatschal C, 12. Wang ET, Cody NA, Jog S, Biancolella M, Wang TT, et al. (2012) et al. (2010) RNA looping by PTB: Evidence using FRET and NMR Transcriptome-wide regulation of pre-mRNA splicing and mRNA localization spectroscopy for a role in splicing repression. Proc Natl Acad Sci U S A 107: by muscleblind proteins. Cell 150: 710–724. 4105–4110. 13. Xue Y, Zhou Y, Wu T, Zhu T, Ji X, et al. (2009) Genome-wide analysis of PTB- 27. Amir-Ahmady B, Boutz PL, Markovtsov V, Phillips ML, Black DL (2005) Exon RNA interactions reveals a strategy used by the general splicing repressor to repression by polypyrimidine tract binding protein. RNA 11: 699–716. modulate exon inclusion or skipping. Molecular cell 36: 996–1006. 14. Yeo GW, Coufal NG, Liang TY, Peng GE, Fu XD, et al. (2009) An RNA code 28. Kafasla P, Mickleburgh I, Llorian M, Coelho M, Gooding C, et al. (2012) for the FOX2 splicing regulator revealed by mapping RNA-protein interactions Defining the roles and interactions of PTB. Biochemical Society transactions 40: in stem cells. Nature structural & molecular biology 16: 130–137. 815–820. 29. Durbin R (1998) Biological sequence analysis : probabalistic models of proteins and 15. Licatalosi DD, Mele A, Fak JJ, Ule J, Kayikci M, et al. (2008) HITS-CLIP yields nucleic acids. Cambridge, UK New York: Cambridge University Press. xi, 356 p. p. genome-wide insights into brain alternative RNA processing. Nature 456: 464– U422. 30. Rabiner LR (1989) A Tutorial on Hidden Markov-Models and Selected 16. Auweter SD, Oberstrass FC, Allain FHT (2006) Sequence-specific binding of Applications in Speech Recognition. Proceedings of the Ieee 77: 257–286. single-stranded RNA: is there a code for recognition? Nucleic Acids Research 31. Hofacker IL (2003) Vienna RNA secondary structure server. Nucleic acids 34: 4943–4959. research 31: 3429–3431. 32. Izquierdo JM, Majos N, Bonnal S, Martinez C, Castelo R, et al. (2005) 17. Keppetipola N, Sharma S, Li Q, Black DL (2012) Neuronal regulation of pre- Regulation of Fas alternative splicing by antagonistic effects of TIA-1 and PTB mRNA splicing by polypyrimidine tract binding proteins, PTBP1 and PTBP2. on exon definition. Molecular cell 19: 475–484. Critical reviews in biochemistry and molecular biology 47: 360–378. 18. Spellman R, Smith CWJ (2006) Novel modes of splicing repression by PTB. 33. Shen HH, Kan JLC, Ghigna C, Biamonti G, Green MR (2004) A single Trends in Biochemical Sciences 31: 73–76. polypyrimidine tract binding protein (PTB) binding site mediates splicing PLOS Computational Biology | www.ploscompbiol.org 17 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets inhibition at mouse IgM exons M1 and M2. Rna-a Publication of the Rna 46. Ray D, Kazan H, Chan ET, Pena Castillo L, Chaudhry S, et al. (2009) Rapid Society 10: 787–794. and systematic analysis of the RNA recognition specificities of RNA-binding proteins. Nature biotechnology 27: 667–670. 34. Ule J, Stefani G, Mele A, Ruggiu M, Wang X, et al. (2006) An RNA map 47. Singh R, Valcarcel J, Green MR (1995) Distinct binding specificities and predicting Nova-dependent splicing regulation. Nature 444: 580–586. functions of higher eukaryotic polypyrimidine tract-binding proteins. Science 35. Xing Y, Stoilov P, Kapur K, Han A, Jiang H, et al. (2008) MADS: a new and 268: 1173–1176. improved method for analysis of differential alternative splicing by exon-tiling 48. Perez I, Lin CH, McAfee JG, Patton JG (1997) Mutation of PTB binding sites microarrays. RNA 14: 1470–1479. causes misregulation of alternative 39 splice site selection in vivo. RNA 3: 764– 36. Carmel I, Tal S, Vig I, Ast G (2004) Comparative analysis detects dependencies among the 59 splice-site positions. RNA 10: 828–840. 49. Licatalosi DD, Yano M, Fak JJ, Mele A, Grabinski SE, et al. (2012) Ptbp2 37. Shapiro MB, Senapathy P (1987) RNA splice junctions of different classes of represses adult-specific splicing to regulate the generation of neuronal precursors eukaryotes: sequence statistics and functional implications in gene expression. in the embryonic brain. Genes & development 26: 1626–1642. Nucleic acids research 15: 7155–7174. 50. Wagner EJ, Garcia-Blanco MA (2001) Polypyrimidine tract binding protein 38. Hosmer DW, Lemeshow S (2000) Applied logistic regression. New York: Wiley. antagonizes exon definition. Molecular and cellular biology 21: 3281–3288. xii, 373 p. p. 51. Carstens RP, Wagner EJ, Garcia-Blanco MA (2000) An intronic splicing silencer 39. Ashiya M, Grabowski PJ (1997) A neuron-specific splicing switch mediated by an causes skipping of the IIIb exon of fibroblast growth factor receptor 2 through array of pre-mRNA repressor sites: evidence of a regulatory role for the involvement of polypyrimidine tract binding protein. Molecular and cellular polypyrimidine tract binding protein and a brain-specific PTB counterpart. biology 20: 7388–7400. RNA 3: 996–1015. 52. Chan RC, Black DL (1995) Conserved intron elements repress splicing of a 40. Spellman R, Llorian M, Smith CWJ (2007) Crossregulation and functional neuron-specific c-src exon in vitro. Molecular and cellular biology 15: 6377– redundancy between the splicing regulator PTB and its paralogs nPTB and ROD1. Molecular cell 27: 420–434. 53. Gooding C, Roberts GC, Smith CW (1998) Role of an inhibitory pyrimidine 41. Parkhomchuk D, Borodina T, Amstislavskiy V, Banaru M, Hallen L, et al. element and polypyrimidine tract binding protein in repression of a regulated (2009) Transcriptome analysis by strand-specific sequencing of complementary alpha-tropomyosin exon. RNA 4: 85–100. DNA. Nucleic acids research 37: e123. 54. Zhang C, Frias MA, Mele A, Ruggiu M, Eom T, et al. (2010) Integrative 42. Wu J, Akerman M, Sun S, McCombie WR, Krainer AR, et al. (2011) modeling defines the Nova splicing-regulatory network and its combinatorial SpliceTrap: a method to quantify alternative splicing under single cellular controls. Science 329: 439–443. conditions. Bioinformatics 27: 3010–3016. 55. Castle JC, Zhang C, Shah JK, Kulkarni AV, Kalsotra A, et al. (2008) Expression 43. Hall MP, Nagel RJ, Fagg WS, Shiue L, Cline MS, et al. (2013) Quaking and of 24,426 human alternative splicing events and predicted cis regulation in 48 PTB control overlapping splicing regulatory networks during muscle cell tissues and cell lines. Nature genetics 40: 1416–1425. differentiation. RNA 19: 627–638. 56. Rozen S, Skaletsky H (2000) Primer3 on the WWW for general users and for 44. Reid DC, Chang BL, Gunderson SI, Alpert L, Thompson WA, et al. (2009) biologist programmers. Methods in molecular biology 132: 365–386. Next-generation SELEX identifies sequence and structural determinants of 57. Schuelke M (2000) An economic method for the fluorescent labeling of PCR splicing factor binding in human pre-mRNA sequence. RNA 15: 2385–2397. fragments. Nature biotechnology 18: 233–234. 45. Daubner GM, Clery A, Jayne S, Stevenin J, Allain FH (2012) A syn-anti 58. Edgar R, Domrachev M, Lash AE (2002) Gene Expression Omnibus: NCBI conformational difference allows SRSF2 to recognize guanines and cytosines gene expression and hybridization array data repository. Nucleic acids research equally well. The EMBO journal 31: 162–174. 30: 207–210. PLOS Computational Biology | www.ploscompbiol.org 18 January 2014 | Volume 10 | Issue 1 | e1003442 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png PLoS Computational Biology Public Library of Science (PLoS) Journal

De Novo Prediction of PTBP1 Binding and Splicing Targets Reveals Unexpected Features of Its RNA Recognition and Function

PLoS Computational Biology , Volume 10 (1) – Jan 30, 2014
Download PDF
18 pages

Loading next page...
 
/lp/public-library-of-science-plos-journal/de-novo-prediction-of-ptbp1-binding-and-splicing-targets-reveals-W1JJ0v0cpg

References (116)

Publisher
Public Library of Science (PLoS) Journal
Copyright
Copyright: © 2014 Han et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: This work was supported by NIH Grant R01 GM49662 and an associated ARRA supplement to DLB. DLB is an Investigator of the Howard Hughes Medical Institute. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.
ISSN
1553-734X
eISSN
1553-7358
DOI
10.1371/journal.pcbi.1003442
Publisher site
See Article on Publisher Site

Abstract

The splicing regulator Polypyrimidine Tract Binding Protein (PTBP1) has four RNA binding domains that each binds a short pyrimidine element, allowing recognition of diverse pyrimidine-rich sequences. This variation makes it difficult to evaluate PTBP1 binding to particular sites based on sequence alone and thus to identify target RNAs. Conversely, transcriptome-wide binding assays such as CLIP identify many in vivo targets, but do not provide a quantitative assessment of binding and are informative only for the cells where the analysis is performed. A general method of predicting PTBP1 binding and possible targets in any cell type is needed. We developed computational models that predict the binding and splicing targets of PTBP1. A Hidden Markov Model (HMM), trained on CLIP-seq data, was used to score probable PTBP1 binding sites. Scores from this model are highly correlated (r =20.9) with experimentally determined dissociation constants. Notably, we find that the protein is not strictly pyrimidine specific, as interspersed Guanosine residues are well tolerated within PTBP1 binding sites. This model identifies many previously unrecognized PTBP1 binding sites, and can score PTBP1 binding across the transcriptome in the absence of CLIP data. Using this model to examine the placement of PTBP1 binding sites in controlling splicing, we trained a multinomial logistic model on sets of PTBP1 regulated and unregulated exons. Applying this model to rank exons across the mouse transcriptome identifies known PTBP1 targets and many new exons that were confirmed as PTBP1-repressed by RT-PCR and RNA-seq after PTBP1 depletion. We find that PTBP1 dependent exons are diverse in structure and do not all fit previous descriptions of the placement of PTBP1 binding sites. Our study uncovers new features of RNA recognition and splicing regulation by PTBP1. This approach can be applied to other multi-RRM domain proteins to assess binding site degeneracy and multifactorial splicing regulation. Citation: Han A, Stoilov P, Linares AJ, Zhou Y, Fu X-D, et al. (2014) De Novo Prediction of PTBP1 Binding and Splicing Targets Reveals Unexpected Features of Its RNA Recognition and Function. PLoS Comput Biol 10(1): e1003442. doi:10.1371/journal.pcbi.1003442 Editor: Roderic Guigo, Center for Genomic Regulation, Spain Received April 1, 2013; Accepted November 23, 2013; Published January 30, 2014 Copyright:  2014 Han et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: This work was supported by NIH Grant R01 GM49662 and an associated ARRA supplement to DLB. DLB is an Investigator of the Howard Hughes Medical Institute. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. * E-mail: [email protected] not currently feasible, in part due to our incomplete understanding Introduction of RNA recognition by the splicing regulators and their Alternative splicing of pre-mRNA commonly determines the mechanisms of action. protein output of mammalian genes, with most genes generating Whole-transcriptome crosslinking methods for individual pro- multiple mRNA and protein products [1]. A typical alternative teins in vivo are allowing the identification of large numbers of exon is affected by multiple pre-mRNA binding proteins that may protein/RNA interaction sites [8,9,10,11]. These data can be either enhance or repress splicing [2]. The expression and activity overlapped with functional data on splicing to identify possible of these splicing regulatory proteins can vary with development, direct target exons for particular proteins [12,13,14,15]. However, cell type, or cellular stimulus [3]. This complex combinatorial there are limitations in the interpretation of these data. Cross- regulation can be seen in the conserved sequences within and linking efficiency can vary between different proteins and between surrounding alternative exons, which generally contain the individual binding sites, making it difficult to relate the crosslinking binding sites for many different regulators. These sequences make signal to the actual binding affinity. These signals are also up what is sometimes called the splicing code as they determine dependent on the expression of the bound RNA, and since these where and when the exon is spliced into an mRNA [4,5,6,7]. Such data are generated one tissue or cell type at a time it is not always a code should allow the development of models that predict exon feasible to extend the results from one setting to a new cell type or regulation based solely on the RNA binding affinity of the many point in development. It would be extremely useful to be able to regulatory proteins and their other interactions. However, this is scan for binding affinity across the complete transcriptome and to PLOS Computational Biology | www.ploscompbiol.org 1 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets elements for each domain are not known. Moreover, the flexible Author Summary linkers separating some of the RRM domains and the requirement A key step in the regulation of mammalian genes is the for a gap between elements simultaneously bound to domains splicing of the messenger RNA precursor to produce a three and four allow for substantial degeneracy in PTBP1 binding mature mRNA that can be translated into a particular sites. This degeneracy and the lack of understanding of the protein needed by the cell. Through the process of sequence features that contribute to binding affinity have made it alternative splicing, mRNAs encoding different proteins difficult to identify PTBP1 binding sites based on sequence alone, can be derived from the same primary gene transcript. The and to assess which sequences surrounding an exon might regulation of this process plays essential roles in the contribute to PTBP1 regulation. development of differentiated tissues and is mediated by Experiments with model substrates indicate that a single high special pre-mRNA binding proteins. To understand how affinity PTBP1 binding site placed upstream of an exon, or within these proteins control gene expression, one must charac- it, can repress splicing [27]. However, strong repression of an terize what they recognize in RNA and identify these efficiently spliced exon requires an additional binding site either binding sites across the genome in order to predict their within the exon or downstream from an exon with an upstream targets. Models that allow this prediction are essential to high affinity site [17,27,28]. PTBP1 is also known to enhance the understanding developmental regulatory programs and splicing of certain exons [13,19,20]. The properties of these exons their perturbation by disease causing mutations. In this and how they differ from those that are repressed by PTBP1 are study, we use statistical methods to build models of RNA unclear, with different studies coming to different conclusions recognition by the important splicing regulator PTBP1 and [13,19]. An analysis of CLIP data in HeLa cells found that PTBP1 then apply these models to predict PTBP1 regulation of new gene transcripts. We show that PTBP1 has different sites near the adjacent constitutive exons could enhance the specificity for RNA than was previously recognized and inclusion of an alternative exon between them [13]. In contrast, that its target exons are more diverse than was known examination of exons whose splicing was reduced by double before. There are many similar splicing regulators in knockdown of Ptbp1 and Ptbp2 found that they frequently had mammalian cells, and these analyses provide a general binding sites immediately downstream [19], whereas splicing framework for the computational analysis of their RNA repression often involved upstream binding sites: a pattern binding and target identification. observed for other splicing regulators. These results are not mutually exclusive. It is possible that the two groups examined predict exon targets in tissues that have not yet been subjected to different subsets of the many exons regulated by PTBP1, and that experimental analysis. the protein may show additional patterns of protein binding Splicing regulatory proteins commonly contain multiple RRM or adjacent to its target exons. other RNA binding domains, with each domain recognizing a short In this study we sought to understand the sequence features that element of a few nucleotides [2,16]. Subtle variation in the optimal determine RNA binding by PTBP1 and to examine how they are binding element of each domain and flexible peptide linkers combined in exons that are targeted by the protein. We first between them allow for significant degeneracy within high affinity developed a statistical model of PTBP1 binding sites that identifies binding sites. Although the short sequence motifs that are common new features of RNA recognition by the protein. This binding to a set of binding sites are readily identified, these likely constitute model was then applied to the assessment of exon regulation by only a portion of a full high affinity site. To rank binding sites and PTBP1 across the transcriptome. assess their finer structures, we need an approach to search for clusters of these short motifs and to score for binding affinity. Results The Polypyrimidine tract binding protein 1 (PTBP1) is a widely G containing triplets contribute to PTBP1 binding studied splicing regulatory protein [17,18]. PTBP1 is known to repress the splicing of a large number of exons by binding in their To examine the interactions of PTBP1 across many binding adjacent introns or within the exons themselves. PTBP1 is down sites, we used a set of PTBP1-bound sequences identified by regulated in differentiating neurons and muscle cells to allow crosslinking immunoprecipitation (CLIP) [13]. PTBP1 has four inclusion of PTBP1 repressed exons during development of these RRMs separated by linker peptides, with each RRM recognizing a tissues [19,20,21]. In neurons the loss of PTBP1 is accompanied by pyrimidine triplet. In previous studies we found that a minimal the up-regulation of the homologous protein PTBP2 [17,20,22]. high affinity binding site for the protein extended across 25 to 30 PTBP2 has similar binding properties to PTBP1 and represses nucleotides, about the average size of the CLIP clusters (29 nt) some of the same exons [23]. Other exons are more sensitive to [27]. Given the triplet recognition and the need for spacers PTBP1 than PTBP2 and are induced to splice when PTBP2 between the direct RRM contacts, it is unlikely that every replaces PTBP1 in early neurons [24]. nucleotide within a CLIP cluster makes a direct base-specific contact with the protein or otherwise contributes to binding PTBP1 contains four RRM domains that recognize short pyrimidine elements [25]. Flexible linkers separate RRM domains affinity. This information about direct binding is hidden in the one and two, and domains two and three. RRM domains three examination of a CLIP tag, but should affect the triplet frequencies and four interact through a hydrophobic interface that position within the entire set of tags. We designed a two-state Hidden their RNA binding surfaces on opposite faces of the two-domain Markov Model (HMM) based on triplets to assess whether triplets structure. This orientation requires that the RNA elements would segregate into two states and whether these two states interacting with the structure be separated by an RNA loop differed in their PTBP1 binding or non-binding potential. The [26]. The structure of each of the PTBP1 RRM domains has been 48,604 CLIP clusters from the human transcriptome were solved in complex with the hexanucleotide, CUCUCU [25]. extracted and used to train the HMM (Figure 1A) [29,30]. This These structures show each domain binding a nucleotide triplet training defined two states showing distinctly different triplet with some additional contacts, and making similar base specific distributions (Figure 1B). Pleasingly, all of the pyrimidine triplets interactions with CU or UC dinucleotides. Other sequences can segregated into State 1. We called this state the PTBP1 binding likely make different base specific contacts, and the optimal state, as we confirm below. We found that 20 triplets have higher PLOS Computational Biology | www.ploscompbiol.org 2 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 1. PTBP1 binding model. A. Scheme of the PTBP1 binding model. The two-state HMM model was trained on PTBP1 bound RNA sequences (48,604 clusters) from published PTBP1-CLIP experiments. Triplets from these CLIP clusters were predictive of two states, with all of the pyrimidine triplets preferred by State 1. The diagram presents the structure of the PTBP1 HMM (Hidden Markov Model) and its trained transition probabilities. B. The probabilities that triplets are seen states 1 or 2 (emission probabilities) are plotted in black and gray bars, respectively. Asterisks indicate G containing pyrimidine triplets. doi:10.1371/journal.pcbi.1003442.g001 probabilities to be seen in the PTBP1 binding state. All triplets of its RRM domains can presumably make specific contacts with containing only pyrimidines were included in this 20-triplet set G residues. On the other hand, all A containing triplets have (Figure 1B), with the top-scoring triplet UCU showing the modest positive emission probabilities for state 2 and are likely to alternating C and U nucleotides seen in many characterized be either neutral or to inhibit PTBP1 binding. PTBP1 binding sites. We next tested the HMM scoring, which strongly weights the Interestingly, multiple triplets containing G residues are also triplets from state 1 over state 2, for prediction of PTBP1 binding. preferred in State1 (Figure 1B). These triplets often contain U We performed cross validation experiments on the Hela CLIP residues as the other nucleotides. Some of these triplets, such as dataset. A background dataset was generated using ten randomly UGU, have output (emission) probabilities in State 1 that are picked sequences from each gene identified as containing a CLIP similar to pyrimidine triplets, presumably also making them cluster. Applying the model to this data set gave us a distribution of predictive of PTBP1 binding. In contrast, triplets containing A scores that was compared to scores generated by subsets of the residues, even if the other two nucleotides are pyrimidines, were all CLIP clusters removed from the training set prior to training. As preferred by the non-PTBP1 binding State 2. These results shown in Figure S1, sequences from subsets of the CLIP clusters indicate that PTBP1 is not strictly pyrimidine specific. At least one scored significantly higher than background. PLOS Computational Biology | www.ploscompbiol.org 3 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets We also tested our model on an independent iCLIP dataset this background model generated scores that predicted affinity from human embryonic stem cells (ESC) (Figure S2). Unlike reasonably well. However, it did generate negative scores for a standard CLIP, iCLIP tags define the probable crosslink site as couple of probes that are shown to bind (data not shown). Thus, the uniform model gave the most accurate scoring of the being the 59 terminus of the tag. We used a Viterbi algorithm to predict the most probable state path predicted by the PTBP1 background models we tested. HMM model for each iCLIP tag. Defining triplets from the State1 The data demonstrate that HMM scoring based on triplet (PTBP1 binding) and triplets from State 2 (nonbinding), we found frequencies can accurately predict the observed binding affinities that the frequency of predicted binding triplets is highly enriched across a wide range of Kd values (from ,250 nM to 1 nM). Probe in the iCLIP cluster regions and peaks precisely at the crosslink 6 yields a z-score of 0.82 and binds with a Kd of 257 nM, whereas site. This indicates that State 1 probability is highly associated with probe 10 scores 2.74 in the model and binds with a Kd of 73 nM PTBP1 crosslinking in vivo. (Figure 2B). These sequences include G containing triplets that To more quantitatively assess the relationship between the contribute to the binding scores. This method allows any sequence to now be quantitatively assessed for possible PTBP1 binding, HMM score and RNA binding, we applied the trained model to a set of 100,000 random 69 nucleotide sequences. This length allows which was not previously possible by simply looking for clusters of a limited number of motifs. This HMM based approach should be for one hexanucleotide binding site for each of the four RRMs applicable to the prediction of binding sites and affinity for other with 15 nucleotide gaps, the minimum gap required for multi-domain RNA binding proteins. simultaneous binding by RRMs 3 and 4 [25,26]. The scores are calculated as a log-odds ratio of the probabilities of the sequence having been generated by the HMM over a background model Placement of PTBP1 binding sites adjacent to target that assigns equal probability to all triplets. The random sequences exons generated a distribution of scores that was used to normalize the With our new method of defining PTBP1 binding sites, we next binding scores, with the average score for random sequence set to examined PTBP1 target exons for the location of predicted PTBP1 zero, and the z-score defined as the deviation from the average as binding. In part, we wanted to reassess two previous studies that shown in Figure S3A [29]. Thus a sequence with a z-score of 2.74 came to differing conclusions regarding the placement of PTBP1 is 2.74 standard deviations from the average (empirical sites adjacent to its target exons. One group mapped PTBP1 CLIP p-value = 0.005), and is predicted to be a significantly stronger clusters adjacent to a limited number of PTBP1 repressed and binder than the average sequence (500 of the 100,000 random enhanced exons [13]. This study described PTBP1 repressed sequences have scores equal or greater than this sequence). A exons as enriched for binding sites both upstream and down- negative z-score is predicted to bind less well than the average stream, as has been seen in studies of individual exons. They did sequence. We isolated thirteen sequences from the mouse not observe PTBP1 CLIP clusters within repressed exons, even transcriptome that exhibited a range of scores from 22.62 to + though such exons have been described [17,32,33]. The PTBP1 4.40 (Figure 2A). These were transcribed in vitro and subjected to enhanced exons they examined showed a trend in PTBP1 binding electrophoretic mobility shift assay to measure binding to near the flanking constitutive exons. A second study examined recombinant PTBP1 (Figure 2B; Figure S3B). Sequences yielding exons showing altered splicing on splicing-sensitive microarrays negative scores all failed to bind PTBP1 within the protein after Ptbp1/Ptbp2 double knockdown [19]. CLIP clusters derived concentration range tested, with the exception of probe 4, which from the first study were mapped to these exons. The authors bound weakly, below the level that would allow measurement of found CLIP cluster enrichment upstream and within PTBP1/ an affinity constant. Positive scoring sequences all yielded PTBP1 PTBP2 repressed exons. In contrast to the previous study, they bound complexes that were assayable by gel shift to derive found that PTBP1/PTBP2 enhanced exons showed enrichment apparent binding affinities. The apparent Kds of these RNAs for CLIP tags in the downstream region. This pattern of binding showed a very strong negative correlation with their binding score site placement relative to repressed and enhanced exons has been from the model (Pearson correlation coefficient =2 0.9), where a observed for several other splicing regulatory proteins [14,34]. higher score predicts a lower Kd and hence a higher affinity In our study, we defined four groups of exons from a set of (Figure 2A). Thus, the scoring system performed very well in exons previously assessed for splicing after Ptbp1 knockdown predicting PTBP1 binding affinity. [20,35]. These included 68 PTBP1-repressed exons whose splicing Two sequences (probes 9 and 11) showed variable binding that increases after Ptbp1 knockdown, 37 PTBP1-enhanced exons shifted their Kd’s slightly off the fitted curve relating z-score to Kd. whose splicing decreases after knockdown, 69 control exons that These may have secondary structures that reduce binding affinity are not affected by Ptbp1 depletion but are known to be thus increase their apparent Kd. To look at this, we examined the alternatively spliced (PTBP1-non regulated), and 1,000 constitu- predicted structure of each probe using the RNA fold program tive exons. We determined the density of predicted PTBP1 [31]. Probes 9 and 11 did not show an overall free energy of binding states within a 24-nucleotide window sliding along the folding substantially lower than other RNAs. However, it is exon region. We also examined the sequence encompassing the difficult to rule out that they contain a local structure that adjacent constitutive exons (Figure 3A). As expected, the non- sequesters some key feature for PTBP1 recognition. regulated control and constitutive exon sets did not exhibit high In addition to the background model using uniform triplet probabilities of PTBP1 binding except in the polypyrimidine tract frequencies, we also tested control sequence sets using different of the 39 splice site. On the other hand, the introns upstream of nucleotide frequencies (Figure S4). Control sets that maintain the PTBP1 repressed exons show enrichment of potential PTBP1 mono or dinucleotide frequencies of the PTBP1 CLIP tags while binding sites starting from 250 nucleotides upstream of the exon. shuffling the triplet frequencies did not perform well. This is not Relative to the control exons, exons repressed by PTBP1 also surprising because these sequences are highly skewed in nucleotide exhibited substantial enrichment of PTBP1 binding sites within the content and the shuffling does not change the triplet frequencies exon itself and within the first 100 nucleotides of the downstream dramatically. We also tested a background model based on intron. The repressed exons thus exhibit binding site placement random sequences selected from genes containing PTBP1 CLIP that combines the findings of the two previous studies [13,19]. The clusters (ten sequences from each gene). Like the random dataset, PTBP1-enhanced exon set also shows enrichment of PTBP1 PLOS Computational Biology | www.ploscompbiol.org 4 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 2. Validation of the PTBP1 binding model. A. To validate binding scores, thirteen RNAs with various PTBP1 binding scores were transcribed in vitro and subjected to binding assay. Apparent Kd’s (dissociation constant) were highly negatively correlated with PTBP1 binding scores (Pearson correlation =20.9). B. Four RNA sequences with predicted PTBP1 binding scores (Full data binding data in Figure S3). Potential PTBP1 binding sites are underlined and in bold. Experimental binding affinities were assessed by electrophoretic mobility shift of RNA by PTBP1 and compared with prediction scores. Apparent dissociation constants (Kd) were defined as the concentration at which half the protein was bound to RNA. doi:10.1371/journal.pcbi.1003442.g002 binding sites within the downstream intron relative to control (Figure 3A). Similar to what was seen in the previous study by exons, although the distribution of binding sites across this region Llorian, we found little enrichment of PTBP1 sites within was different between the repressed and enhanced exon sets enhanced exons [19]. There is a limited enrichment adjacent to PLOS Computational Biology | www.ploscompbiol.org 5 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 3. Sequence characteristics of PTBP1-dependent alternatively spliced exons. A. An RNA map shows enrichment of predicted PTBP1 binding sites near PTBP1-dependent exons. The Y-axis plots average density of predicted PTBP1 binding states within a 24 nt window; the length of overlap between two adjacent windows was 8 nt. B. To assess PTBP1 binding signatures of individual exons, known PTBP1 regulated exons were clustered by their PTBP1 binding score profiles and visualized as heat maps. These heat maps indicate wide variation in the positions of PTBP1 binding sites between individual exons. C. Four sequence features including the PTBP1 binding scores and 39 splice site strength show statistically significant differences between regulated and control exon groups (one-tailed Student’s t-tests). doi:10.1371/journal.pcbi.1003442.g003 PLOS Computational Biology | www.ploscompbiol.org 6 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets the exons flanking enhanced exons. Interestingly however, we find to consider many different factors. Nevertheless, models based on some PTBP1 enhanced exons that have PTBP1 binding sites single factors will be useful for understanding the relative upstream of the exon. These were not seen in either previous contributions of individual proteins to patterns of splicing study. Our results are generally consistent with the known regulation. Such models will be easier to interpret regarding the placement of PTBP1 binding sites in PTBP1 target exons and contributions of individual factors to individual exons than more imply that rules correlating the position of PTBP1 binding to its complex models. Moreover in the longer term, models developed effect on a target exon are not as strict as seen for some other for different individual factors can be combined to make more splicing regulators. The mechanisms proposed from previous maps accurate predictions. To assess how well one might model splicing of PTBP1 binding do not appear to be generalizable to all PTBP1 regulation by a single factor, we examined whether the strength targets [13,19,27]. and placement of predicted PTBP1 binding sites could be used to Binding maps for PTBP1 and other splicing regulators show the predict new PTBP1 dependent exons. We plotted the scores for a averages of multiple exons. Since the data indicated a high level of variety of sequence features against the percent of exons exhibiting variability in binding site placement between individual exons, we that score that also exhibit PTBP1 dependent exon repression wanted to visualize target exons relative to each other. To display (Figure S5). These plots produced distinct sigmoidal curves where binding signals for individual exons we created heat maps of the most exons regulated by PTBP1 were found above or below a binding scores upstream, within, and downstream of each exon in particular score. This strongly suggests that a logistic regression the PTBP1 target set (Figure 3B). This display makes clear that the model incorporating each of these scores will be predictive of location of PTBP1 binding sites within its known target exons is PTBP1 repression. variable. We found that 60% of PTBP1 repressed exons are We developed a multinomial logistic regression model and predicted to have strong binding sites within the upstream intron. trained it on three classes of regulated exons (Figure 4A) [38]. The Most of these exons also have strong binding sites within either the training set included PTBP1 repressed exons, PTBP1 enhanced exon or the downstream intron, patterns that were observed exons, and non-regulated exons. Each exon in each class was previously [13,19,27]. However, other patterns of binding site scored for the four features found to correlate with PTBP1 placement are also seen, suggesting PTBP1 dependent exons are regulation (x through x ), including the 39 splice site strength, and 1 4 following multiple rules. Some repressed exons score highly for the PTBP1 binding scores for each of three regions: the 250 PTBP1 binding only within the exon or in both the exon and the nucleotides upstream of the exon, the exon itself, and the 100 downstream intron. About half of PTBP1 enhanced exons have nucleotides downstream of the exon. These intron lengths strong PTBP1 binding sites downstream (Figure 3B). These can encompass the regions of binding site enrichment for PTBP1 co-occur with upstream intron-binding sites, but rarely with exon dependent exons (Figure 3). binding sites. Interestingly, there are exons enhanced by PTBP1 The PTBP1-enhanced exons are fewer in number and show with strong upstream binding in the absence of other sites. These more limited enrichment of PTBP1 binding sites than PTBP1- data demonstrate the heterogeneity in the position of PTBP1 repressed exons making the prediction for these exons less accurate. binding sites for its target exons. This heterogeneity needs to be We first tested models that considered just PTBP1-repressed exons considered for predicting PTBP1 dependent regulation. relative to control groups. However, we found that including the PTBP1 repressed exons exhibited significantly higher average enhanced exons as a separate training group improved the binding scores in both the upstream intron and in the exon itself, prediction of repressed exons, even though enhanced exons than either the control group of alternative exons or the PTBP1 themselves are not as easily identified (data not shown). enhanced exons (Figure 3C). The average binding scores in the The trained model yielded values for the b coefficients that downstream introns were higher for both the PTBP1-repressed weight the different features contributing to the regulation. As and PTBP1-enhanced exons than the control group (Figure 3C), expected the upstream binding score was weighted most heavily in although not at the same statistical significance. The variability of predicting PTBP1 repression (Table S1), although binding scores in binding site placement within the smaller group of PTBP1- all three regions contributed to the score for PTBP1 repression. In enhanced exons presumably contributes to the weaker statistical contrast, we found that only the downstream binding score was correlation of binding scores with positive regulation. significantly associated with PTBP1 enhancement. The upstream We also compared the three exon sets for other features that score generated a b coefficient close to zero making it essentially might contribute to their ability to be regulated by PTBP1, neutral in the prediction of enhanced exons. The exon binding score including exon length, flanking intron length, and 59 and 39 splice was subject to a negative b coefficient, indicating that exon binding site strength. Most of these features were not statistically different reduces the probability of PTBP1 enhancement. Using these b among the three-exon groups. However, both PTBP1 enhanced coefficients, the trained models for repression or enhancement each and PTBP1 repressed exons were found to carry significantly yield a value of the g-function (logit) for an exon (x) given by the log weaker 39 splice sites than the control exon set, as measured by the of the ratio of the probability of repression or enhancement over the Analyzer Splice Tool (Figure 3C) [36,37]. probability that the exon is not regulated. From this, the probability These results indicate that PTBP1-repressed exons, and perhaps that an exon is repressed by PTBP1 can be determined from the two PTBP1-enhanced exons, exhibit an ensemble of sequence features g-values as shown in Figure 4A. that define them as PTBP1 regulated and that should allow their We assessed the multinomial logistic regression model by identification by sequence alone. recursively retraining on exon sets with one exon left out and then scoring the missing exon. This leave-one-out cross validation Prediction of PTBP1 repressed exons enabled assessment of the overall performance of the model [38] Alternative exons are generally regulated by multiple factors (Figure S6). The PTBP1 dependent exon repression logit showed that act both positively and negatively on their ability to be spliced. good prediction, with an area under the curve (AUC) value of Thus, an exon controlled by a regulator in one context might not 0.72, substantially greater than random guessing (AUC = 0.5). As be affected by it under other conditions where counteracting expected, the enhanced exon logit was not as accurate as the factors are present, or required cofactors are absent. This means repression logit (AUC = 0.57), although it was better than random that the most accurate predictions of splicing regulation will need (Figure S6A). Using these data, we assessed the sensitivity and PLOS Computational Biology | www.ploscompbiol.org 7 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 4. Scheme of the PTBP1 splicing regulation model and its application to an exon in Ptbp3.A. The PTBP1 splicing regulation model was trained on known PTBP1-regulated and non-regulated exons and used to predict new PTBP1-dependent exons. Prediction results were compared to changes in exon inclusion (PSI) measured by RT-PCR and RNA-seq. An exon from Ptbp3 is presented as a prediction example. From intron and exon sequences, PTBP1 binding scores and 39 splice site strength were calculated and fed into the regulation model. B. The model predicts exon 2 of Ptbp3 as repressed by PTBP1 with high probability (0.89). Ptbp1 knockdown in mouse neuroblastoma cells (N2A) confirmed de-repression of the exon (from PSI = 45 to PSI = 70). doi:10.1371/journal.pcbi.1003442.g004 specificity across the range of scores to define a decision threshold numbers of repressed exons. We sought to choose a threshold that for exon repression scores (Figure S6B). Increasing the threshold gave a low false positive rate over one that yielded more regulated increases the specificity by eliminating many false positives, but exons. We found that above a threshold score of 0.65 the false decreases the sensitivity of the model in identifying maximum positive rate was 10% or lower (Figure S6B). PLOS Computational Biology | www.ploscompbiol.org 8 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Applying the model to 4494 alternative cassette exons from (2.3%) with probability scores above 0.65 for being repressed by UCSC genome browser database, we found 243 exons (5.4%) that PTBP1. Among other activities, these exons were enriched in yielded a PTBP1 repression probability score greater than 0.65 and genes that function in calcium ion transport, cytoskeletal which were not in the training set. The 50 top-scoring cassette exons organization, intracellular transport, and synaptic transmission, are listed in Table 1. These included two exons that were reported all functions affected by previously known PTB targets (Table S2). previously to be PTBP1 targets. An exon of Gabrg2 yields a To assess splicing of this large set of predicted PTB targets, we probability score of 0.92. Although we could not confirm its used RNA-seq to generate a large dataset of exons that change after repression in N2A cells because of low expression of the transcript, Ptbp1 knockdown. RNA from control and PTBP1-depleted N2A the orthologous exon in rat is a well-characterized PTBP1 cells was subjected to high density short read sequencing on the repression target [39]. Exon 2 of Ptbp3 (Rod1), another known Illumina HiSeq platform using a strand specific, paired end protocol PTBP1 target [40], yielded a repression probability score of 0.89 [41]. Exons whose inclusion changed between the two samples were and was confirmed by RT/PCR to show increased inclusion after identified by alignment to an exon database and quantification of Ptbp1 knockdown (Figure 4B). We performed additional RT-PCR exon inclusion using the SpliceTrap program [42]. After filtering for validation in triplicate on a series of high and low scoring exons from read coverage and removing the training set, we identified 573 transcripts expressed in N2A cells (Figure 5 & Figures S7, S8 and alternative exons whose splicing was assayable in N2A cells. These S9). Seven of ten exons scoring above 0.65 were de-repressed after exons exhibit changes in percent exon inclusion (delta PSI) ranging Ptbp1 knockdown in N2A cells, yielding a validation rate of 70%. from 229% to 62% upon PTBP1 depletion. The exons were The actual false positive rate is difficult to estimate because exons binned by their PTBP1 repression probability scores and plotted for with high repression scores that are not affected by Ptbp1 depletion their change in PSI (Figure 6). The average changes in splicing were in N2A cells might be regulated by PTBP1 in other cells. An significantly correlated with the repression probability. Exons indication that this might be occurring is that the average inclusion scoring below 0.5 distributed around zero change in PSI, but above level (or percent spliced in value, PSI) of the putative false positives is this score the average exon inclusion is altered by PTBP1 depletion. significantly higher than the confirmed true positives in N2A cells, Most notably, exons with a repression probability score above 0.65 indicating that they will be less prone to change upon Ptbp1 exhibited significantly larger changes in splicing than exons with depletion and be more difficult to validate (Figure S8B). Thus, the lower scores. Exons with intermediate scores and hence weaker true positive rate may be greater than 70%. Importantly, the high binding sites show smaller changes in splicing than high scoring validation rate for exons scoring above 0.65 indicates that the exons. Setting a threshold of a 5% change in PSI as validation, 22 of binding model and the regulation model based upon it can identify 33 exons (67%) that scored above 0.65 for PTBP1 regulation were many new PTBP1 targets that were not previously known (Table1). confirmed as PTBP1 repression targets in N2A cells. At least some High scoring exons might also fail to be validated because of of the other 11 exons are presumably PTBP1 targets in other cells. regulation by other proteins. Knockdown of Ptbp1 induces To test the model in another cell type, we examined exons expression of its close homolog Ptbp2, which targets some of the reported to change after Ptbp1 knockdown in mouse C2C12 same exons [20] (Figure S7). To test whether PTBP2 was also myoblasts, as measured on splicing sensitive microarrays [43]. Very targeting the predicted PTBP1 repressed exons, we knocked down similar to what was observed in N2A cells, we found that exons with Ptbp2 or both Ptbp1 and Ptbp2 expression in N2A cells and re- high repression probabilities showed significant de-repression upon assayed the exons in triplicate (Figures S10, S11 & S8A). Although the Ptbp1 knockdown compared to exons with low repression some exons showed greater inclusion in the double knockdown probabilities (Figure S12). Of 29 exons assayed on the arrays with a compared to depletion of Ptbp1 alone, this did not validate any repression probability above 0.65, 19 exons were confirmed as additional predicted PTBP1 repressed exons. We did identify some PTBP1 repressed on the array (q-value,0.05), yielding a validation high and low scoring exons showing more complex regulation by rate of 66%. Thus the model performed very similarly in C2C12 the two PTB proteins (Figure S10 & S11). and N2A cells. Among the 11 high scoring exons identified as We also examined a set of low scoring exons (probability score# unchanged after PTBP1 knockdown in N2A cells only 3 were 0.2) by RT-PCR after Ptbp1 and/or Ptbp2 depletion (Figure 5B assayed on the array and expressed in C2C12 cells. These again and Figure S11). All of these exons (8 of 8) failed to respond to the showed high inclusion in C2C12 prior to knockdown and so were loss of PTBP1 and are likely true negatives. Thus, PTBP1 difficult to assay for derepression. Thus, it is difficult to use the repression scores above 0.65 and below 0.2 were highly predictive C2C12 data to draw conclusions about the false positive rate. for regulation and its absence, respectively. As expected, interme- The logistical model gives us a new tool for studying the diate scores were less consistent in their predictive value (Figure regulation of alternative splicing. Using it, we can now scan genomic S9). Some exons in the intermediate scoring group were affected sequence to score exons for PTBP1 regulation. Applying the model by PTB proteins and will be interesting to assess further. genomewide, the PTBP1 repression probability scores were The prediction of PTBP1-repressed exons was improved by integrated into the UCSC genome browser. These data, displayed treating PTBP1-enhanced exons as a separate class, but the with the RNAseq data from N2A cells are available at our website probability scores for PTBP1 enhancement did not consistently (http://www.mimg.ucla.edu/faculty/black/ptbatweb/). A novel identify new PTBP1 target exons (data not shown). This is likely in PTBP1 repressed exon in the Kcnq2 gene is shown in Figure 6B. part due to the smaller number of exons in the training set and The logistic model thus allows the assessment of any exon across the their heterogeneity, with some possibly being indirect targets. transcriptome for likely PTBP1 regulation. These predictions will likely improve with training on larger numbers of PTBP1 enhanced exons as they are identified. Discussion However, it is possible that simply the presence of the PTBP1 binding site is not sufficient for predicting PTBP1 enhancement New features of PTBP1 binding sites and that binding sites for other factors will need to be considered. We have developed two computational models, one that allows We next tested the model on a genomewide scale, by applying it accurate prediction of PTBP1 binding sites and another that to a set of 168,111 mouse internal exons and ranking them by their predicts likelihood of PTBP1 repression of exons across the probability of PTBP1 repression. This analysis yielded 3824 exons transcriptome. These models uncovered several new features of PLOS Computational Biology | www.ploscompbiol.org 9 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Table 1. PTBP1 repressed exons identified by the splicing model. PTB Binding Scores Gene Name Gene Description mm9 coordinates 39 p(Repressed) Upstream Downstream Splice site Intron Intron (250 nt) Exon (100 nt) Strength Pax6 paired box gene 6 chr2:105523985–105524115(+) 8.35 21.49 20.80 0.27 0.99 Mbd5 methyl-CpG binding chr2:49134101–49135303(+) 6.46 20.93 2.27 20.32 0.98 domain protein 5 Arhgap24 Rho GTPase activating chr5:102981145–102981338(+) 6.47 0.10 20.05 0.46 0.97 protein 24 Tle1 transducin-like enhancer chr4:71819247–71819451(2) 4.71 0.05 21.26 22.56 0.94 of split 1 Acsl6 acyl-CoA synthetase chr11:54150438–54150515(+) 4.16 1.40 20.03 20.82 0.94 long-chain family Ryr1 ryanodine receptor 1, chr7:29829938–29829955(2) 4.78 20.08 20.88 21.71 0.94 skeletal muscle Ankhd1 ankyrin repeat and KH chr18:36784163–36784921(+) 4.37 0.03 1.44 20.64 0.93 domain containing 1 Slc39a14 solute carrier family 39 chr14:70713408–70713577(2) 3.51 1.16 21.27 23.20 0.92 (zinc transporter) Gabrg2 gamma-aminobutyric chr11:41727472–41727495(2) 1.95 2.77 0.56 23.79 0.92 acid (GABA) A receptor Itga7 integrin alpha 7 chr10:128378878–128378997(+) 4.14 0.29 1.13 20.25 0.92 Iqsec2 IQ motif and Sec7 chrX:148615540–148615635(+) 4.88 0.68 20.13 0.71 0.91 domain 2 Smarca2 SWI/SNF related, chr19:26825612–26825646(+) 3.94 20.09 1.23 20.86 0.91 matrix associated, actin dependent regulator of chromatin Zfand3 zinc finger, AN1-type chr17:30197755–30197795(+) 4.17 2.34 0.90 1.80 0.91 domain 3 Agap2 ArfGAP with GTPase chr10:126527198–126527257(+) 3.57 0.06 20.53 23.08 0.90 domain, ankyrin repeat and PH domain 2 Ttn Titin chr2:76723554–76723832(2) 2.93 1.06 1.19 21.83 0.90 Ptbp3 ROD1 regulator chr4:59559021–59559054(2) 3.80 0.57 1.66 0.73 0.89 of differentiation 1 (S. pombe) Mapk8 mitogen-activated chr14:34203859–34203930(2) 2.35 1.17 1.01 23.48 0.89 protein kinase 8 Snap91 synaptosomal-associated chr9:86693373–86693534(2) 2.60 1.89 20.35 22.17 0.88 protein 91 Fmnl1 formin-like 1 chr11:103059449–103059547(+) 3.93 20.60 21.36 22.70 0.88 Phldb1 pleckstrin homology-like chr9:44509029–44509169(2) 3.20 0.57 1.05 20.66 0.87 domain, family B 2310035C23Rik RIKEN cDNA 2310035C23 chr1:107637012–107637094(+) 2.03 1.76 0.85 22.46 0.87 gene Arnt aryl hydrocarbon receptor chr3:95270715–95270759(+) 3.48 20.36 2.53 0.09 0.87 nuclear translocator Smyd2 SET and MYND domain chr1:191723697–191723807(2) 3.33 0.33 20.28 21.64 0.86 containing 2 Ap2a1 adaptor protein complex chr7:52158832–52158897(2) 3.35 20.12 20.89 22.69 0.86 AP-2, alpha 1 subunit Klra killer cell lectin-like chr6:130329011–130329100(2) 2.82 3.18 0.62 0.77 0.86 receptor, subfamily A Spag9 sperm associated chr11: 93942054–93942068(+) 0.99 3.01 1.62 23.03 0.86 antigen 9 Col4a3bp collagen, type IV, chr13:97386949–97387026(+) 2.81 1.23 0.74 20.93 0.86 alpha 3 binding protein Garnl3 GTPase activating chr2:32941395–32941464(2) 4.15 0.38 0.36 0.94 0.86 RANGAP domain-like 3 PLOS Computational Biology | www.ploscompbiol.org 10 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Table 1. Cont. PTB Binding Scores Gene Name Gene Description mm9 coordinates 39 p(Repressed) Upstream Downstream Splice site Intron Intron (250 nt) Exon (100 nt) Strength Dennd1a DENN/MADD domain chr2:37982049–37982168(2) 3.37 0.80 1.35 0.59 0.86 containing 1A Ms4a7 membrane-spanning chr19:11400297–11400353(2) 2.79 2.35 0.37 20.05 0.86 4-domains, subfamily A BC030307 cDNA sequence BC030307 chr10:86169981–86170089(+) 2.75 20.16 1.95 22.40 0.85 Phactr1 phosphatase and actin chr13:43154940–43155146(+) 2.73 1.25 0.50 20.97 0.85 regulator 1 R3hdm2 R3H domain containing 2 chr10:126902187–126902240(+) 1.66 1.98 1.96 21.57 0.84 Cdc14b CDC14 cell division chr13:64306579–64306725(2) 1.42 2.75 2.44 20.58 0.84 cycle 14B Ubqln1 ubiquilin 1 chr13:58282183–58282266(2) 2.88 0.98 20.06 21.17 0.84 Ttn Titin chr2:76739898–76740179(2) 2.63 20.07 1.49 22.38 0.84 Stx3 syntaxin 3 chr19:11857290–11857400(2) 3.00 21.12 2.26 23.62 0.84 Slc8a3 solute (sodium/calcium) chr12: 82310340–82310458(2) 1.84 1.25 1.76 22.25 0.84 carrier family 8 Zfp62 zinc finger protein 62 chr11:49028057–49028156(+) 3.27 1.98 20.51 0.19 0.83 Dlg1 discs, large homolog 1 chr16:31771843–31771941(+) 1.53 1.98 1.87 21.65 0.83 (Drosophila) Nrxn2 neurexin II chr19:6463824–6463847(+) 3.35 21.37 1.33 22.26 0.83 Klra7 killer cell lectin-like chr6:130179953–130180042(2) 2.68 2.11 20.39 20.63 0.83 receptor, subfamily A Picalm phosphatidylinositol chr7:97330729–97330878(+) 1.15 2.15 2.30 22.37 0.83 binding clathrin assembly Acad8 acyl-Coenzyme A chr9:26798168–26798277(2) 2.61 0.88 20.31 21.86 0.83 dehydrogenase family Epn1 epsin 1 chr7:5033620–5033723(+) 3.92 0.65 0.07 1.06 0.82 Grip1 glutamate receptor chr10:119422530–119422685(+) 2.66 20.74 2.61 23.13 0.82 interacting protein 1 Csmd3 CUB and Sushi multiple chr15:47587514–47587627(2) 2.42 2.16 0.20 20.45 0.82 domains 3 Lrrfip1 leucine rich repeat (in FLII) chr1:92990137–92990214(+) 2.02 0.40 3.44 22.21 0.82 interacting protein 1 Srsf11 serine/arginine-rich chr3:157703405–157703586(+) 1.09 2.05 20.56 24.75 0.82 splicing factor 11 Tmem209 transmembrane chr6:30441087–30441184(2) 3.82 0.16 0.10 0.64 0.82 protein 209 The 50 highest scoring exons predicted to be repressed by PTBP1 based on sequence alone. doi:10.1371/journal.pcbi.1003442.t001 RNA recognition by PTBP1 and the properties of its target exons. The base-specific contacts that PTBP1 makes with Guanosine The PTBP1 binding model was based on triplets following the are not yet clear. Recent studies of RNA recognition by SRSF2 structures of the PTBP1 RRM domains, whose sequence specific (SC35) protein have shown that the element GGAG can be contacts are each primarily to three nucleotides. We find that the recognized by the same RRM as CCAG by flipping the initial two set of triplets that increase the probability of binding includes the G nucleotides to the syn conformation [45]. It will be very expected pyrimidine motifs, particularly those with alternating interesting to investigate whether a similar anti to syn switch cytosines and uridines. However, many triplets with guanosine occurs in RNA bound by PTBP1, when C residues are replaced residues also increase binding probability. In contrast, adenosine with G. residues have a negative effect on binding. Thus, RNA recognition Previous characterizations of PTBP1 binding sites have focused by PTBP1 is not solely dependent on pyrimidine nucleotides. The on finding enriched short motifs within populations of bound recognition of G residues by PTB was unexpected, although some RNAs or regulated exon sequences [13,44,46,47,48]. These previously characterized PTB binding sites did contain G residues methods generally identify elements whose short length will allow [13,44]. With this model, we can now predict PTBP1 binding interaction with only one RRM domain. Searching for new affinity to any site in the transcriptome. binding sites comprised of clusters of these short elements can PLOS Computational Biology | www.ploscompbiol.org 11 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets PLOS Computational Biology | www.ploscompbiol.org 12 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 5. Validation of novel PTBP1-repressed exons by RT-PCR. A. Candidate PTBP1-repressed exons with probability greater than 0.65 were validated by RT-PCR following Ptbp1 knockdown. Data shown are averages 6 standard error of PSI (Percent of Spliced In) from biological triplicates. Statistical analysis was performed using paired one-tailed Student’s t-test (p-values,0.01**, ,0.05*). B. Exons with low PTBP1 repression probabilities (#0.2) were also validated by RT-PCR following Ptbp1 knockdown in biological triplicates. doi:10.1371/journal.pcbi.1003442.g005 identify higher affinity sites but does not consider all elements or to achieve accurate predictions from the model. However, treating rank them. Crosslinking-immunoprecipitation experiments allow these as a separate exon class improves the prediction of PTBP1 large numbers of binding regions to be identified. However, not all repression. We find that for probability scores above 0.65 the the sequence within a CLIP tag will be contacting the protein and model is strongly predictive of PTBP1 repression. Applying this it is difficult to relate CLIP signals to binding affinity. The HMM criterion across the transcriptome, we identified hundreds of new allowed the individual assessment of different short elements PTBP1 target exons. within the CLIP clusters, showing that they segregated into two Alternative exons are generally regulated by multiple proteins states. The ranking of the triplets for their contributions to one of acting in combination, and a particular exon will often be subject these states yielded a model where complex clusters of short to both positive and negative regulation by antagonistic factors. elements could be assessed for binding and yielded accurate For a model based on one factor, these other proteins will predictions of binding affinity. Many RNA binding proteins are confound predictions. Exons with high PTBP1 binding scores may similar to PTBP1 in having multiple domains that may each make be counteracted by antagonistic factors in some cell types. different base specific contacts with RNA. The widespread Alternatively, synergistic factors may allow an exon with a generation of CLIP-seq datasets will allow the modeling of RNA relatively weak binding site to still recruit PTBP1. Thus, a model recognition by almost any protein based on a large number of based on one factor will be limited in its predictive power. In this known binding sites. study, our intent was to measure the effect of PTBP1 binding alone Using the same modeling approach, we also developed a before considering the contributions of other factors. The logistic binding model for PTBP2 (neuronal PTB) using a published modeling allowed the contributions of different binding site PTBP2 CLIP dataset [49]. PTBP2 is about 70% identical to placements to PTBP1 regulation to be measured. PTBP1 in sequence, and has only two amino acid changes among Several studies have used Bayesian models to dissect the the residues making direct contact with RNA [17]. We found that regulatory properties of exons [7,54]. These models can generate the binding models for two PTB proteins were also nearly identical accurate predictions by incorporating a wide variety of sequence, indicating that the two proteins are likely to differ more in their expression and conservation data. However, because so many protein/protein interactions than in their RNA binding sites (Data disparate variables are incorporated, it can be difficult to draw not shown). mechanistic conclusions from these models regarding any one protein. For example, the presence of high pyrimidine density upstream from the branch point can be predictive of exons Defining PTBP1 target exons showing neuronal specific inclusion [7,55]. This is presumably in Several PTBP1 target exons have been analyzed in detail part due to many neuronal exons being regulated by PTBP1 and [17,50]. These exons vary in the placement and action of their PTBP2. However, a subset of these exons may be regulated by PTBP1 binding sites. It is common for PTBP1-repressed exons to other factors with pyrimidine rich binding sites. In the long term, it have a binding site upstream, often encompassing the branch will be most accurate to develop predictive binding models for point of the 39 splice site [39]. Exons can also be repressed by each protein, similar to the PTBP1 model here, and then to PTBP1 binding within the exon [19,32,33]. Other exons contain incorporate each of these binding models into a larger network downstream binding sites that are needed in conjunction with an model. Such an approach will allow the analysis of the many upstream site to achieve splicing repression [51,52,53]. Although overlapping regulatory programs controlled by RNA binding acting as a repressor for most of its targets, PTBP1 also activates proteins. the splicing of a group of exons. There have been divergent reports about placement of PTBP1 binding sites needed to mediate PTBP1 enhancement of splicing. The PTBP1 binding model Materials and Methods allowed us to examine PTBP1 binding site placement across a Hidden Markov Model for PTBP1 binding affinity large set of known PTBP1 target exons. Nearly all exons had prediction predicted high affinity PTBP1 binding sites nearby. We found that A Hidden Markov Model (HMM) was designed and trained by more than half of PTBP1 repressed exons have high affinity binding sites upstream, and a fraction of PTBP1 enhanced exons an expectation–maximization (EM) method (Baum-Welch algo- have high affinity sites downstream. These exons fit with recent rithm) using published PTBP1 CLIP data [13,29,30]. In total, results on several other splicing regulators where the placement of 48,604 PTBP1-CLIP cluster sequences were used to train model the binding site determines the direction of the regulatory effect parameters. During the training step, multiple initial values were [12,14,34]. However, for PTBP1 these rules are not so clear. Some tested to avoid a local maximum problem. Trained parameters PTBP1 repressed exons have their strongest predicted binding site included emission probabilities for nucleotide triplets, initial downstream or within the exon. These results indicate that there probabilities and transition probabilities between states [29,30]. are fundamental differences between the mechanisms of PTBP1 The trained model was used to score RNA sequences. The raw mediated splicing regulation, and those governing regulation by PTBP1 binding score is defined as a log-odds ratio that compares certain other splicing factors. the score of a sequence from the HMM over the score from a To quantify the predictive value of the PTBP1 binding scores background model. Since CLIP experiments do not have an inherent corresponding negative dataset, we generated computa- for PTBP1 repression, we built a logistic model for PTBP1 regulation. For exons repressed by PTBP1, binding scores for the tional negative datasets and tested different background models upstream, downstream and exon sequences all contribute to the (Figure S4). We found that a background model that values all probability of repression. Exons enhanced by PTBP1 were too few triplets equally yielded the most accurate binding scores [29]. Raw PLOS Computational Biology | www.ploscompbiol.org 13 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets Figure 6. Large-scale validation of novel PTBP1-repressed exons by RNA-seq. A. Validation of the PTBP1 splicing model using RNA-seq. After Ptbp1 knockdown, we performed RNA-seq experiments and estimated changes in PSI (Percent of Spliced In) for 573 cassette exons. The graph shows average delta PSI values for exons, grouped by their probabilities to be repressed by PTBP1. The number of exons in the corresponding probability bin is given by n. P-values were calculated from one-tailed Student’s t-test. B. A genome browser screenshot of a novel PTBP1-regulated exon: exon 2 of the Kcnq2 gene. For whole internal mouse exons, we created custom genome browser tracks to visualize the PTBP1 splicing model and mapped RNA seq reads. doi:10.1371/journal.pcbi.1003442.g006 PLOS Computational Biology | www.ploscompbiol.org 14 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets scores were further normalized and converted to z-scores. For the with canonical splice sites (GU-AG). An exon was classified as 69 mer RNA sequences used in binding assays, scores were PTBP1 repressed or enhanced when 1) the inclusion level (PSI) of normalized by 100,000 random sequences with same length its minor isoform was greater than 5% in both the control and (Figure S3). This yielded very accurate predictions of binding knock-down samples and 2) the inclusion level of its minor isoform affinity (Figure 2). was changed by 30% or more in the Ptbp1 knock down condition When considering binding scores in genomic sequence, exons compared to the control sample. Next, we collected sequence features for each exon and its flanking exons. The features and upstream or downstream intron regions have different base compositions and will yield different average binding scores. Thus, included PTBP1 binding scores, 59 and 39 splice site strengths, exon/intron lengths, and word frequencies. The PTBP1 binding to score binding sites adjacent to possible regulated exons, it is scores were calculated from the PTBP1 binding model described more informative to score sites relative to equivalent sequence above. The strength of splice sites was calculated by the splice-site regions. From the annotated mouse genome, we retrieved 168,111 analyzer tool [37]. Using a mouse whole internal exon set, we internal exons and their flanking introns as separate sequence sets normalized features and fed them into the model. The PTBP1 using a python library, Pygr. We scored log odds of these sequences splicing model is based on a multinomial logistic regression with the trained model. Since the lengths and base compositions of framework using the following steps: 1) selection of initial variables intronic and exonic sequences are different, and binding scores with a moderate level of association (p-value from t-test,0.25), 2) automatically increase with length (Figure S13) [29], we grouped removal of outlier exons, 3) stepwise variable selection [38]. We sequences by their location and sequences in each group were scored mouse internal exons with the trained PTBP1 splicing sorted according to length into bins of 1000 sequences each. The model and validated candidate exons with RT-PCR and RNA-seq average score and standard deviation were determined for each experiments. Exons from the training set were excluded from the bin. These values were used to transform the raw scores into validation. z-scores for each upstream intron, downstream intron, and exon sequence. We localized the PTBP1 binding sites along each RNA sequence using the Viterbi algorithm [29,30]. Validation of exon candidates by RT-PCR and RNA-seq To test alternative splicing events for candidate exons, we assayed exon inclusion levels in cells following Ptbp1, Ptbp2, and Validation of PTBP1 binding model scores by binding double Ptbp1 & Ptbp2 knock down. The knockdown experiment assay was performed as described previously with minor modification To test predicted PTBP1 binding scores, we selected thirteen [20]. Mouse neuroblastoma (N2A) cells were cultured in DMEM mouse exon/intron RNA sequences (69 nucleotides) exhibiting a with 10% FBS and 2 mM L-glutamine. At 70 to 80% confluency, range scores. In the selection, other sequence features such as cells were trypsinized and suspended in the growth medium. secondary structure were not considered. Target RNAs were DNA–Lipofectamine 2k (Invitrogen) complexes were prepared transcribed in vitro from dsDNA using T7 RNA polymerase and and mixed with cells in a tube according to manufacturer’s subjected to an electrophoretic mobility shift assay (EMSA). instructions. Tubes were incubated for 5 h with mixing every half During the transcription, radioactive a-32P UTP was incorporat- hour. Then cells were centrifuged and cultured in plates for 3 d. ed into RNA to visualize the probes. The RNA probes were then Proteins and RNA was extracted from collected cells. Protein denatured for 2 min at 85uC and cooled down on ice immediately samples were subjected to fluorescence immunoblotting to to reduce secondary structure formation. Binding assays were monitor knockdown efficiency of Ptbp1 and Ptbp2. Total RNA carried out as previously described with some modifications [27]. was collected using Trizol (Invitrogen) according to the manufac- Specifically, each gel mobility shift reaction (10 mL) contained the turer’s instructions. The RNA was further treated with DNase I to indicated amounts of recombinant human PTBP1 in 6 mLDG avoid DNA contamination. For RT-PCR (Reverse Transcription- buffer (20 mM Hepes-KOH ph 7.9, 20% glycerol, 80 mM PCR) assays, the RNA was reverse transcribed to cDNA with potassium glutamate, 0.2 mM EDTA, 0.2 mM PMSF), 1 mL random hexamers using SuperScript enzyme (Invitrogen) follow- 22 mM MgCl2, 1 mL 0.5 mg/ml tRNA, 0.5 mL RNase inhibitor ing the manufacturer’s instructions. PCR reactions were per- (20 unit, RNaseOut from invitrogen), 0.5 mL DEPC treated H2O, formed to assay alternative splicing of particular target exons. and 1 mL 100 nM RNA probe. At first, all reaction components First, forward and reverse PCR primers were designed for the excluding RNase inhibitor, tRNA, and RNA probes were mixed flanking exons using PRIMER3 program [56]. To label PCR and incubated for 8 min at 30uC. Then RNase inhibitor and products, a 59 fluorescent-labeled universal primer (59-FAM- tRNA were added and mixed. RNA probe was then added and CGTCGCCGTCCAGCTCGACCAG-39) was added to the PCR the reaction was incubated for an additional 15 min. The reaction and a universal priming site was introduced to the 59 end reactions were put on ice for 5 min and mixed with 1.2 mL of the forward primer (59-CGTCGCCGTCCAGCTCGACCAG- glycerol loading dye (30% glycerol). They were separated on 8% Forward Primer-39). Each PCR reaction (15 mL) was carried out native polyacrylamide gels with 25 mM Tris-Gly running buffer in with 1.5 picomole of the forward primer and 6.75 picomole of the a cold room. Gels were dried and exposed to a phosphor screen. reverse and universal primers [57]. PCR amplification proceeded Then images were scanned using Typhoon 9410 and quantified with an initial denaturation at 94uC for 4 m followed by 24 cycles using ImageQuant TL program (GE Lifesciences). The apparent of 94uC for 30 s, at a melting temperature of the reverse primer for Kd values were estimated by fitting the data to non-linear curves 45 s, and 72uC for 45 s, with a final extension step at 72uC for using Prism software. 10 m. The samples were mixed with 26 formamide buffer (Formamide with 1 mM EDTA pH 8.0) and denatured at 95uC Logistic regression model for PTBP1 dependent exon for 5 min. Then samples were chilled on ice and run on 8% prediction denaturing polyacrylamide gels. Gels were directly scanned by An exon training set was compiled from previous microarray Typhoon and quantified by ImageQuant program. and RT-PCR experiments [20,35]. The training set was composed RNA-seq libraries were constructed following standard proto- with 68 PTBP1 repressed, 37 PTBP1 enhanced, and 69 non- cols (Illumina TruSeq RNA Sample Prep Kit). To make strand- PTBP1 regulated simple cassette exons. We only considered exons specific libraries, we added two extra steps to the protocol [41]. PLOS Computational Biology | www.ploscompbiol.org 15 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets After first strand cDNA synthesis, remaining dNTPs were slightly improves the linear fit (20.91 vs. 20.95) for some strong removed by a size selection on beads (AMPure XP). Second- binders. However, it wrongly predicted some binders as non- strand cDNA was synthesized with a dNTP mix containing dUTP binders, which reduced the rank correlation (20.95 to 20.90). instead of dTTP. The reaction contained samples eluted in 50 ml The two shuffled models did not perform well. resuspension buffer, 2 ml56 FS buffer, 1 ml 50 mM MgCl2, 1 ml (TIF) 100 mM DTT, 2 ml 10 mM dUTP nucleotides mix, 15 ml Second Figure S5 Correlations of particular sequence features Strand Buffer (Invitrogen), 0.5 ml E.coli DNA Ligase (10 U/ with PTBP1 repression. For PTBP1-repressed and PTBP1 ml;NEB), 0.5 ml RNase H (2 U/ml;Invitrogen), 2 ml DNA E.coli non-regulated exons, we calculated scores for sequence features Polymerase I (10 U/ml;NEB). The reaction was incubated for 2 h and determined the fraction of PTBP1-repressed exons in each at 16uC. After sequencing adaptors were ligated, 1 ml USER score bin. Shown are the graphs for the 39 splice site score, and the (Uracil-Specific Excision Reagent enzyme; NEB) was added to PTBP1 binding scores in the upstream intron, the exon, and the reactions to degrade the second strand cDNA. The samples were downstream intron plotted against the percent of exons within the incubated for 15 min at 37uC and the reaction were inactivated at score bin that are PTBP1-repressed. 94uC for 5 min. The samples were put in ice and then subjected to (TIF) PCR amplification. Average size of inserts was about 225 bp and the libraries were subjected to 100 bp paired-end sequencing Figure S6 Performance of PTBP1-dependent splicing (Illumina HiSeq2000 platform). Using SpliceTrap [42], 60–65% models. A. Receiver Operating Characteristic Curves in a leave- of reads were mapped to exon duos or trios. In total, 180M one-out cross validation for each logit: exon repression (left) and (179,511,116) and 145M (145,334,711) paired end reads were exon enhancement (right). B. Sensitivity and specificity plotted used to infer exon inclusion ratios in the control and Ptbp1 across the whole threshold range. Sensitivity is defined as the knockdown conditions, respectively. The data have been deposited percent of true repressed exons that are correctly predicted as in NCBI’s Gene Expression Omnibus [58] and are accessible repressed at the corresponding threshold. Specificity is defined as through GEO Series accession number GSE45119. the percent of actual non-repressed exons that are correctly predicted as non-repressed at the corresponding threshold. Supporting Information (TIF) Figure S1 Five-fold cross-validation of the PTBP1 Figure S7 ShRNA mediated depletion of PTBP1 and binding model using Hela CLIP clusters. To test the two PTBP2. Duplicate immunoblots after shRNA knockdown of state model of binding and non-binding triplets, we divided the PTBP1, PTBP2 or both proteins. Note that depletion of PTBP1 CLIP-data to the five subsets. In each plot, four sub sets were used induces expression of PTBP2 as observed previously. Numbers in training the model and one subset was subjected to scoring. We above each lane indicate the fluorescence intensity for PTBP1 or then compared scores from the CLIP-subset sequences to random PTBP2 relative to the control lane. sequences picked from same genic regions that contained the (TIF) CLIP clusters. As shown, sequences from CLIP-subset generated Figure S8 Characteristics of false positive exons. A. significantly higher scores than random. The results indicate that Exon inclusion was measured for three false positives exons after the triplets identified by the HMM as predictive of state 1 are Ptbp1 knockdown (left), or Ptbp2 knockdown, and Ptbp1 & Ptbp2 predictive of PTBP1 CLIP sites and thus of protein binding. double knockdown (right). P-values were calculated from biolog- (TIF) ical triplicates using paired one-tailed t-tests. B. False positive Figure S2 The density of PTBP1 binding triplets exons exhibit higher PSI values prior to PTBP1 depletion. Box predicted by the HMM peaks at the aligned crosslink plot of exon inclusion for twenty-nine false positive exons showing sites from Human ESC iCLIP clusters. little change in splicing by RNA-seq after PTBP1 depletion (delta (TIF) PSI,5%) that score with high probability to be repressed (.0.55). Exon inclusion levels prior to PTBP1 depletion are compared to Figure S3 PTBP1 binding model scores and validation. thirty-six true positive exons. A. Summary statistics and the distribution of raw and normalized (TIF) PTBP1 binding scores for 100,000 random sequences. B. Electrophoretic mobility shift assay of RNAs with various PTBP1 Figure S9 Two exons with intermediate scores for binding scores. RNAs were transcribed in vitro, incubated with PTBP1 repression show complex responses to PTBP1 increasing concentrations of purified PTBP1 (0 to 200 nM), and and PTBP2 depletion. Exon inclusion was measured after the bound and unbound RNA separated on native gels. Arrows Ptbp1 depletion (left), or after Ptbp2 and Ptbp1/Ptbp2 double indicate RNA-protein complexes. The fraction of PTBP1-bound depletion (right). P-values were calculated from biological RNA is plotted below for each RNA. triplicates using paired one-tailed t-tests. (TIF) (TIF) Figure S4 Evaluation of different background models Figure S10 PTBP2 dependence of predicted PTBP1 for scoring PTBP1 binding. Four background models were target exons. RT-PCR of high probability PTBP1 exon targets evaluated. The uniform distribution model assumes equal following Ptbp2 knockdown or Ptbp2/Ptbp1 double knockdown. frequencies of triplets. The PTBP1 target gene set model used Relative band intensities of the gels in triplicate on the right are random sequences from genes containing PTBP1 CLIP clusters. plotted on the left to show the average delta PSI 6 SE (Percent The two shuffled models used shuffled CLIP cluster sequences Spliced In). P-values were calculated from paired one-tailed t-tests maintaining mono or di nucleotide ratios. We calculated PTBP1 with PSI values in control samples. binding scores based on each background model and compared (TIF) the scores to the measured dissociation constants. Based on the rank correlation, the uniform distribution model worked best. The Figure S11 PTBP2 dependence of predicted non-PTBP1 PTBP1 target gene set model showed comparable performance. It target exons. RT-PCR of exon with low probabilities for PTBP1 PLOS Computational Biology | www.ploscompbiol.org 16 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets repression (#0.2). Ptbp2 knockdown and Ptbp1/Ptbp2 double knock Table S1 Trained PTBP1 splicing regulation model. The table presents a summary of the multinomial logistic regression down with data analysis as in Figure S10. (TIF) model for PTBP1 splicing regulation, including estimated coefficients and their statistics. Figure S12 Boxplot of PSI (Percent of Spliced In) values (TIF) estimated from splicing sensitive microarray data for Table S2 Enriched gene ontology categories for novel exons expressed in mouse C2C12 myoblasts. Exons with PTBP1-repressed exons. The table lists ontology entries higher (.0.65) and lower (,0.2) repression probabilities are enriched in genes with predicted PTBP1-repressed exons (prob- compared. Exons used in the original training set were excluded ability score of exon repression .0.65). Whole mouse internal from the plot. The number of exons in the corresponding exons were used as the control set, and p-value cut off was 0.05. probability bin is given by n. The p-value was calculated from a Gene ontology analysis was performed using the GOTM web one-tailed Student’s t-test. server. (TIF) (TIF) Figure S13 Distribution of PTBP1 binding scores of exons and introns before and after normalization. Acknowledgments 168,111 exons and their flanking introns from the set of annotated We thank Manny Ares, Matteo Pellegrini, Yi Xing, and Christopher Lee mouse internal exons were subjected to scoring and normalization. for thoughtful suggestions on the experiments, the analysis and manuscript, Raw PTBP1 Binding scores are affected by sequence length and and all members of the Black laboratory for helpful discussion. We thank base composition. To account differences in these features Shalini Sharma and Batoul Amir-Ahmady for advice on the EMSA between introns and exons in normalized scores, we grouped the experiments, and Paul Boutz for advice on the PTBP1 knock down exons and their upstream and downstream introns separately. The experiments. We especially thank Shalini Sharma for the recombinant sequences in each group were sorted according to length into bins PTBP1 proteins. of 1,000 sequences each. The average scores and standard deviations were determined for each bin. These values were used Author Contributions to transform the raw scores into z-scores for each sequence per Conceived and designed the experiments: AH DLB. Performed the bin. experiments: AH PS AJL. Analyzed the data: AH DLB. Contributed (TIF) reagents/materials/analysis tools: YZ XDF. Wrote the paper: AH DLB. References 1. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, et al. (2008) Alternative 19. Llorian M, Schwartz S, Clark TA, Hollander D, Tan LY, et al. (2010) Position- isoform regulation in human tissue transcriptomes. Nature 456: 470–476. dependent alternative splicing activity revealed by global profiling of alternative 2. Black DL (2003) Mechanisms of alternative pre-messenger RNA splicing. splicing events regulated by PTB. Nature structural & molecular biology 17: Annual review of biochemistry 72: 291–336. 1114–1123. 20. Boutz PL, Stoilov P, Li Q, Lin CH, Chawla G, et al. (2007) A post- 3. Gabut M, Chaudhry S, Blencowe BJ (2008) SnapShot: The splicing regulatory transcriptional regulatory switch in polypyrimidine tract-binding proteins machinery. Cell 133: 192 e191. reprograms alternative splicing in developing neurons. Genes & development 4. Matlin AJ, Clark F, Smith CWJ (2005) Understanding alternative splicing: 21: 1636–1652. Towards a cellular code. Nature Reviews Molecular Cell Biology 6: 386–398. 5. Blencowe BJ (2006) Alternative splicing: New insights from global analyses. Cell 21. Boutz PL, Chawla G, Stoilov P, Black DL (2007) MicroRNAs regulate the 126: 37–47. expression of the alternative splicing factor nPTB during muscle development. Genes & development 21: 71–84. 6. Wang Z, Burge CB (2008) Splicing regulation: from a parts list of regulatory 22. Makeyev EV, Zhang J, Carrasco MA, Maniatis T (2007) The MicroRNA miR- elements to an integrated splicing code. RNA 14: 802–813. 124 promotes neuronal differentiation by triggering brain-specific alternative 7. Barash Y, Calarco JA, Gao W, Pan Q, Wang X, et al. (2010) Deciphering the pre-mRNA splicing. Mol Cell 27: 435–448. splicing code. Nature 465: 53–59. 8. Ule J, Jensen K, Mele A, Darnell RB (2005) CLIP: a method for identifying 23. Zheng S, Gray EE, Chawla G, Porse BT, O’Dell TJ, et al. (2012) PSD-95 is protein-RNA interaction sites in living cells. Methods 37: 376–386. post-transcriptionally repressed during early neural development by PTBP1 and PTBP2. Nature Neuroscience 15: 381–U202. 9. Darnell RB (2010) HITS-CLIP: panoramic views of protein-RNA regulation in 24. Tang ZZ, Sharma S, Zheng SK, Chawla G, Nikolic J, et al. (2011) Regulation of living cells. Wiley interdisciplinary reviews RNA 1: 266–286. the Mutually Exclusive Exons 8a and 8 in the CaV1.2 Calcium Channel 10. Hafner M, Landthaler M, Burger L, Khorshid M, Hausser J, et al. (2010) PAR- Transcript by Polypyrimidine Tract-binding Protein. Journal of Biological CliP–a method to identify transcriptome-wide the binding sites of RNA binding Chemistry 286: 10007–10016. proteins. Journal of visualized experiments : JoVE. 25. Oberstrass FC, Auweter SD, Erat M, Hargous Y, Henning A, et al. (2005) 11. Konig J, Zarnack K, Rot G, Curk T, Kayikci M, et al. (2010) iCLIP reveals the function of hnRNP particles in splicing at individual nucleotide resolution. Structure of PTB bound to RNA: specific binding and implications for splicing Nature structural & molecular biology 17: 909–915. regulation. Science 309: 2054–2057. 26. Lamichhane R, Daubner GM, Thomas-Crusells J, Auweter SD, Manatschal C, 12. Wang ET, Cody NA, Jog S, Biancolella M, Wang TT, et al. (2012) et al. (2010) RNA looping by PTB: Evidence using FRET and NMR Transcriptome-wide regulation of pre-mRNA splicing and mRNA localization spectroscopy for a role in splicing repression. Proc Natl Acad Sci U S A 107: by muscleblind proteins. Cell 150: 710–724. 4105–4110. 13. Xue Y, Zhou Y, Wu T, Zhu T, Ji X, et al. (2009) Genome-wide analysis of PTB- 27. Amir-Ahmady B, Boutz PL, Markovtsov V, Phillips ML, Black DL (2005) Exon RNA interactions reveals a strategy used by the general splicing repressor to repression by polypyrimidine tract binding protein. RNA 11: 699–716. modulate exon inclusion or skipping. Molecular cell 36: 996–1006. 14. Yeo GW, Coufal NG, Liang TY, Peng GE, Fu XD, et al. (2009) An RNA code 28. Kafasla P, Mickleburgh I, Llorian M, Coelho M, Gooding C, et al. (2012) for the FOX2 splicing regulator revealed by mapping RNA-protein interactions Defining the roles and interactions of PTB. Biochemical Society transactions 40: in stem cells. Nature structural & molecular biology 16: 130–137. 815–820. 29. Durbin R (1998) Biological sequence analysis : probabalistic models of proteins and 15. Licatalosi DD, Mele A, Fak JJ, Ule J, Kayikci M, et al. (2008) HITS-CLIP yields nucleic acids. Cambridge, UK New York: Cambridge University Press. xi, 356 p. p. genome-wide insights into brain alternative RNA processing. Nature 456: 464– U422. 30. Rabiner LR (1989) A Tutorial on Hidden Markov-Models and Selected 16. Auweter SD, Oberstrass FC, Allain FHT (2006) Sequence-specific binding of Applications in Speech Recognition. Proceedings of the Ieee 77: 257–286. single-stranded RNA: is there a code for recognition? Nucleic Acids Research 31. Hofacker IL (2003) Vienna RNA secondary structure server. Nucleic acids 34: 4943–4959. research 31: 3429–3431. 32. Izquierdo JM, Majos N, Bonnal S, Martinez C, Castelo R, et al. (2005) 17. Keppetipola N, Sharma S, Li Q, Black DL (2012) Neuronal regulation of pre- Regulation of Fas alternative splicing by antagonistic effects of TIA-1 and PTB mRNA splicing by polypyrimidine tract binding proteins, PTBP1 and PTBP2. on exon definition. Molecular cell 19: 475–484. Critical reviews in biochemistry and molecular biology 47: 360–378. 18. Spellman R, Smith CWJ (2006) Novel modes of splicing repression by PTB. 33. Shen HH, Kan JLC, Ghigna C, Biamonti G, Green MR (2004) A single Trends in Biochemical Sciences 31: 73–76. polypyrimidine tract binding protein (PTB) binding site mediates splicing PLOS Computational Biology | www.ploscompbiol.org 17 January 2014 | Volume 10 | Issue 1 | e1003442 Prediction of PTBP1 Binding and Splicing Targets inhibition at mouse IgM exons M1 and M2. Rna-a Publication of the Rna 46. Ray D, Kazan H, Chan ET, Pena Castillo L, Chaudhry S, et al. (2009) Rapid Society 10: 787–794. and systematic analysis of the RNA recognition specificities of RNA-binding proteins. Nature biotechnology 27: 667–670. 34. Ule J, Stefani G, Mele A, Ruggiu M, Wang X, et al. (2006) An RNA map 47. Singh R, Valcarcel J, Green MR (1995) Distinct binding specificities and predicting Nova-dependent splicing regulation. Nature 444: 580–586. functions of higher eukaryotic polypyrimidine tract-binding proteins. Science 35. Xing Y, Stoilov P, Kapur K, Han A, Jiang H, et al. (2008) MADS: a new and 268: 1173–1176. improved method for analysis of differential alternative splicing by exon-tiling 48. Perez I, Lin CH, McAfee JG, Patton JG (1997) Mutation of PTB binding sites microarrays. RNA 14: 1470–1479. causes misregulation of alternative 39 splice site selection in vivo. RNA 3: 764– 36. Carmel I, Tal S, Vig I, Ast G (2004) Comparative analysis detects dependencies among the 59 splice-site positions. RNA 10: 828–840. 49. Licatalosi DD, Yano M, Fak JJ, Mele A, Grabinski SE, et al. (2012) Ptbp2 37. Shapiro MB, Senapathy P (1987) RNA splice junctions of different classes of represses adult-specific splicing to regulate the generation of neuronal precursors eukaryotes: sequence statistics and functional implications in gene expression. in the embryonic brain. Genes & development 26: 1626–1642. Nucleic acids research 15: 7155–7174. 50. Wagner EJ, Garcia-Blanco MA (2001) Polypyrimidine tract binding protein 38. Hosmer DW, Lemeshow S (2000) Applied logistic regression. New York: Wiley. antagonizes exon definition. Molecular and cellular biology 21: 3281–3288. xii, 373 p. p. 51. Carstens RP, Wagner EJ, Garcia-Blanco MA (2000) An intronic splicing silencer 39. Ashiya M, Grabowski PJ (1997) A neuron-specific splicing switch mediated by an causes skipping of the IIIb exon of fibroblast growth factor receptor 2 through array of pre-mRNA repressor sites: evidence of a regulatory role for the involvement of polypyrimidine tract binding protein. Molecular and cellular polypyrimidine tract binding protein and a brain-specific PTB counterpart. biology 20: 7388–7400. RNA 3: 996–1015. 52. Chan RC, Black DL (1995) Conserved intron elements repress splicing of a 40. Spellman R, Llorian M, Smith CWJ (2007) Crossregulation and functional neuron-specific c-src exon in vitro. Molecular and cellular biology 15: 6377– redundancy between the splicing regulator PTB and its paralogs nPTB and ROD1. Molecular cell 27: 420–434. 53. Gooding C, Roberts GC, Smith CW (1998) Role of an inhibitory pyrimidine 41. Parkhomchuk D, Borodina T, Amstislavskiy V, Banaru M, Hallen L, et al. element and polypyrimidine tract binding protein in repression of a regulated (2009) Transcriptome analysis by strand-specific sequencing of complementary alpha-tropomyosin exon. RNA 4: 85–100. DNA. Nucleic acids research 37: e123. 54. Zhang C, Frias MA, Mele A, Ruggiu M, Eom T, et al. (2010) Integrative 42. Wu J, Akerman M, Sun S, McCombie WR, Krainer AR, et al. (2011) modeling defines the Nova splicing-regulatory network and its combinatorial SpliceTrap: a method to quantify alternative splicing under single cellular controls. Science 329: 439–443. conditions. Bioinformatics 27: 3010–3016. 55. Castle JC, Zhang C, Shah JK, Kulkarni AV, Kalsotra A, et al. (2008) Expression 43. Hall MP, Nagel RJ, Fagg WS, Shiue L, Cline MS, et al. (2013) Quaking and of 24,426 human alternative splicing events and predicted cis regulation in 48 PTB control overlapping splicing regulatory networks during muscle cell tissues and cell lines. Nature genetics 40: 1416–1425. differentiation. RNA 19: 627–638. 56. Rozen S, Skaletsky H (2000) Primer3 on the WWW for general users and for 44. Reid DC, Chang BL, Gunderson SI, Alpert L, Thompson WA, et al. (2009) biologist programmers. Methods in molecular biology 132: 365–386. Next-generation SELEX identifies sequence and structural determinants of 57. Schuelke M (2000) An economic method for the fluorescent labeling of PCR splicing factor binding in human pre-mRNA sequence. RNA 15: 2385–2397. fragments. Nature biotechnology 18: 233–234. 45. Daubner GM, Clery A, Jayne S, Stevenin J, Allain FH (2012) A syn-anti 58. Edgar R, Domrachev M, Lash AE (2002) Gene Expression Omnibus: NCBI conformational difference allows SRSF2 to recognize guanines and cytosines gene expression and hybridization array data repository. Nucleic acids research equally well. The EMBO journal 31: 162–174. 30: 207–210. PLOS Computational Biology | www.ploscompbiol.org 18 January 2014 | Volume 10 | Issue 1 | e1003442

Journal

PLoS Computational Biology – Public Library of Science (PLoS) Journal

Published: Jan 30, 2014

There are no references for this article.