A correlation with exon expression approach to identify cis-regulatory elements for tissue-specific alternative splicing

Identification of muscle-enriched alternative exon and control exon datasets

Total RNA from three biological replicates (three separate individuals) of 16 normal adult human tissues was purchased from BioChain (Hayward, CA, USA). Labeled target was generated from ∼200 ng of total RNA and hybridized to a prototype version of the Affymetrix Human Exon Array as described (26). The set of microarrays contain ∼1.4 million probesets designed to interrogate, as comprehensively as possible, more than 1 million exon clusters derived from a variety of input sources including annotated genes, cDNA sequences and exon prediction algorithms. Design information and microarray data is available at the GEO database (http://www.ncbi.nlm.nih.gov/geo/; accession number: GSE5791).

Candidate muscle-enriched probesets were identified using the splicing index approach (26,36,37). Exon-level expression was normalized to the expression level of the parent gene by dividing probeset intensities by the median intensity of probesets from exons supported by RefSeq or Ensembl annotations. Exons that exhibited statistically significant differences in inclusion rate were identified using a student's t-test on the gene-level normalized probeset intensities (NI). NI values from the three biological replicates of heart and skeletal muscle tissues were compared as a group to the replicates of 14 other non-muscle tissues as a second group. The magnitude of inclusion rate change (splicing index) was estimated by calculating a log ratio (base 2) of the median muscle NI and the median non-muscle NI (26,37). After filtering out non-expressed probesets and genes with low expression, probesets with t-test P-values <0.001 and splicing index magnitudes of >0.5 were considered candidates for muscle-enriched exons.

Manual filtering of the initial list was performed to select further for high confidence internal cassette exons, by mapping candidate muscle-enriched probeset to their genomic context using the BLAT tool (38) at the UCSC genome browser (http://genome.ucsc.edu). Probesets that overlapped annotated alternative transcriptional starts, alternative polyadenylation sites, or regions with alternative 5′ or 3′ splice sites, were removed from consideration in this study. Exon-level probeset intensities were additionally observed using BLIS (Biotique Systems, Inc. Reno, NV), an integrated genome browser that enables exon expression data from the microarray to be viewed in genomic context. Only probesets that showed clear patterns of muscle enrichment were kept for further analysis. Probesets had to demonstrate higher intensity levels in the muscle tissues and have exon-level data for surrounding probesets consistent with exon skipping in a majority of non-muscle tissues. Probesets were subsequently mapped to the May 2004 human genome (NCBI Build 35) using BLAT (38). Exact exon boundaries were determined by comparison to EST and mRNA sequences requiring consensus splice sites.

For phylogenetic analysis, the orthologous exons were identified in another mammalian genome (mouse; Mus musculus), in an avian genome (chicken; Gallus gallus) and in an amphibian genome (frog; Xenopus tropicalis) using VISTA alignment tools. Automatic alignment was successful at finding most of the longer alternative exons directly, but in a few cases the alignments were adjusted manually. The upstream 200 nt (U200) intronic region was selected as the base 1 to base 200 adjoining the exon in the upstream direction, while downstream 200 nt (D200) intronic region was selected as the base 1 to base 200 downstream of the exon. Alignments of orthologus introns and exons sequences were generated by LAGAN using default parameters (39).

The ‘tissue-non-specific alternative’ exon dataset was derived as described previously (8) from the European Bioinformatics Institute database of human alternative exons (http://www.ebi.ac.uk/asd/altextron/index.html). ‘Control exon datasets’ were generated from randomly selected chromosomal regions by extraction from RefSeq annotation databases to get exon coordinates. Control groups for the mammalian and chicken genomes were described previously (8). The muscle-enriched datasets and the control datasets is available at: http://vision.lbl.gov/People/ddas/NAR_SPLICE1/

Validation of muscle-enriched expression

A random subset of candidate muscle-enriched exons was selected for validation by RT–PCR, focusing (for ease of amplification) on those ⩽155 nt in length. RNAs from different human tissues, including heart, skeletal muscle and six non-muscle sources, were purchased from Clontech. One microgram of each RNA source was transcribed into cDNA using random hexamer primers in a total volume of 10 μl. Then, 2 μl cDNA was amplified in a volume of 25 μl, using primers located in the flanking constitutive exons (Supplementary Table 2), for 35 cycles under the following conditions: 30 s at 94°C; 30s at 55°C; 45 s at 72°C. The identity of PCR products was confirmed by DNA sequence analysis.

Correlation with expression

Linear correlation

Counts of hexamers were obtained in a specific pre-mRNA sequence region (upstream or downstream proximal intron). For each region, a linear model was fitted between the logarithm of ratios of gene-normalized exon expression levels and count of each 6-mer word w across a set of exons, :

NI_e is the gene-normalized expression level of exon e in muscle, and C refers to a reference sample. The reference data was taken as the average NI across all tissues. The coefficients a_w and b_w were obtained by a least squares fit. P-values were calculated using an F-test, as described previously (40).

The best fit was obtained for a set of sequences that included the muscle-specific exons (foreground set) and a background set of m sequences (m = 300), drawn randomly from a set of manually curated 957 cassette exons across the human genome (11). Since we started with a prioritized set of tissue-enriched sequences, a background set was necessary to model the correct dependence of log ratios on word count. n such random draws were performed (n = 25), and a linear fit was obtained for each such draw. A geometric mean of the P-values from all iterations reflects the overall significance of the word.

Linear splines

Linear splines differ from lines by introducing a threshold, called knot, below which the function is constant and linear above it (19,22). A significant difficulty in modeling binding sites via PWMs is that they give rise to a continuous distribution of scores across all possible binding sites. Consequently, a cutoff score needs to be determined to discriminate the true sites from false sites. Such cutoff scores are often based on predetermined background sequences and thresholds, and as a result, are complicated by subjective choices (41). In a linear spline model, the cutoff score corresponds to the knot, and thus, is learnt directly from the input data (22). For each PWM μ of width L, each L-mer in the input sequence was assigned a probability score M:

where p_i(b_i) is the probability of observing the base b_i at the position i. Thus, the score M always assumes a value between 0 and 1. It is related to binding affinity (42). PWM scores across exons for a given motif μ were fitted to the splicing ratios {log (NI_e / NI_eC)} using the following model:

where θ(x,0) is a linear spline: it is x, when x⩾0, and zero, otherwise. ξ_μ, termed knot, corresponds to the cutoff score. The coefficients a_μ and b_μ and the location of the knot ξ_μ were determined by a least squares fit. This leads to an unbiased and adaptive determination of the knot ξ_μ for any given PWM. Importantly, in contrast to previous approaches (22), where contribution from only the maximum scoring site was considered, we systematically accounted for contribution from active sites with weaker scores as well. The number of such active sites is adaptively learnt, as displayed in the equation above. Thus, both binding affinity and occurrences of active motifs are accounted for in our approach. The significance of the fit was assessed using an F-test (40). The overall significance of each PWM was enumerated using the same iterative procedure as for the linear regression discussed above.

Over-representation analysis

Identification and characterization of muscle-enriched alternative exons

The human muscle-enriched exon dataset analyzed in this study (Supplementary Table 1) was derived from exon microarray hybridization data using a platform designed to provide a comprehensive genome-wide analysis of annotated and predicted exons (see Methods section). In order to identify motifs that regulate tissue-specific alternative splicing, it is critical to identify a set of alternative exons having similar expression patterns indicative of regulation by a shared splicing program. Therefore, the group of exons studied here was carefully selected by analysis of exon microarray data from a panel of 16 normal adult human tissues. Probesets that exhibited gene-level NI that were significantly higher in heart and skeletal muscle, relative to 14 other tissues, were first identified. For this part of the analysis, we grouped the heart and skeletal muscle exon expression together to enhance the power of the statistical tests. Then a manual filtering process was performed so as to retain only probesets representing cassette exons, and to eliminate probesets corresponding to alternative first and last exons or to exon regions generated from alternative 5′ and 3′ splice sites. The final dataset consisted of 56 muscle-enriched, internal cassette exons. Most of these exons (∼80%) are integral multiples of 3 nt in length, with a median length of 84 nt, consistent with the notion that alternative exons are smaller than average constitutive exon length [∼145 nt (47,48)]. However, the genes with such alternative exons have a median size of 123 kb, much longer than the average gene length. To explore evolutionary conservation of candidate splicing regulatory elements, we also identified highly conserved orthologs for most of these human muscle-enriched exons in mouse, chicken and frog (Supplementary Table 1). It is important to note that while many of these exons show evidence of alternative splicing in Genbank, most were not previously known to exhibit muscle-enriched splicing and were not identified in the pilot study of muscle-enriched exons by Sugnet et al. (11). Therefore, analysis of this dataset should yield novel insights into the vertebrate muscle alternative splicing program, and should provide an opportunity to explore computationally the regulatory motifs that carry out this program.

Muscle-enriched splicing patterns for a random subset of these exons were validated experimentally in the human dataset by RT-PCR (Figure 1). Although splicing patterns were not absolutely muscle specific, in almost every case the efficiency of exon inclusion was highest in heart and skeletal muscle, confirming the predictions of the exon microarray. Importantly, mRNA and/or EST evidence from the genetic databases (data not shown) demonstrates that the majority of these exons are alternatively spliced in at least one of the other species examined (mouse, chicken or frog), suggesting that the incidence of conserved alternative exons in this specialized dataset is higher than the reported rate for general alternative exons (49). Taken together, these results indicate that the muscle-enriched exons constitute a special class of highly conserved alternative exons.

Intron sequences flanking orthologous alternative exons in the mouse and human genomes tend to be evolutionarily conserved (24), consistent with the observation that cis-regulatory elements for tissue-specific alternative splicing are often located in those proximal intron regions. We used VISTA genome alignment tools to compare the proximal intron sequences in this muscle-enriched dataset and extended the evolutionary comparison to include chicken and frog. In the proximal 200-nt upstream (U200) and downstream (D200) introns, mouse sequences were highly similar to their human orthologs (median identity of 61 and 58%, respectively), while chicken and frog introns were much less homologous. The full quantitative data are shown in Supplementary Table 3 and representative alignments of exons with relatively high conservation (FXR1), or lower conservation restricted mainly to the exon (LRRFIP1), are displayed in Figure 2. The reduced overall homology of chicken and frog introns suggests that conserved motifs in these regions are likely conserved specifically for their function as cis-regulatory elements for muscle-specific alternative splicing, rather than being passively conserved as part of a larger conserved element.

Frequent occurrence of muscle-enriched exons in genes encoding proteins with functions in cytoskeletal organization

Previous studies have demonstrated that the brain-specific alternative splicing factor, NOVA1, modulates the splicing of many components of the neuronal synapse (50). We hypothesized that the muscle alternative splicing program might similarly coordinate the expression of a particular class of genes that share a common pathway or cellular process. Using the method described previously (22,51), to examine the gene ontology (GO) terms associated with each parent gene for the muscle-enriched exons, we found a strong association with cytoskeleton organization and biogenesis, microtubule stabilization and muscle development (Supplementary Table 4). These associations were statistically significant (P < 0.001), suggesting that the muscle alternative splicing program is critical for proper expression of the unique cytoskeleton characteristic of vertebrate muscle.

Correlation with exon expression identifies splicing regulatory elements

Alternative splicing regulatory elements responsible for tissue-specific splicing are often located in proximal intron sequences (25). To search for candidate intronic regulatory motifs for the muscle-specific splicing program, we correlated the frequencies of hexamers in specific intronic regions with the logarithm of ratios of gene-normalized exon expression levels in skeletal muscle, across the 56 muscle-enriched exons in the human dataset. The ratio for any exon was enumerated against its average gene-normalized expression level across all the tissues. Thus, it is similar to the splicing index used above. The cis-elements exhibiting significant correlation with expression were considered potentially functional in regulating muscle-specific splicing. These were further examined for relative over-representation in introns of muscle-enriched exons, compared to a background set of introns flanking constitutive exons, using a hypergeometric distribution based on word counting in the oligonucleotide sequences. Finally, we examined their spatial conservation through vertebrate evolution (8) by testing whether the motif is over-represented in the other species using exactly the same statistical measures as used for humans.

Here we consider ugcaug in the downstream 200 nt (D200) of intron sequence as an example. This hexamer represents the binding site for mammalian Fox-1 and Fox-2 splicing factors (31), which have identical RRM domains. ugcaug has been reported as a common motif in proximal introns adjacent to tissue-specific exons. In a few cases, functional splicing assays have confirmed the importance of this motif in regulation of splicing (27–33). In the large group of muscle-enriched exons studied here, we found a highly significant correlation of ugcaug frequency with muscle expression (P = 6.8E−05). The distribution and the linear fit for a single iteration of correlation analysis (see Methods section) are shown in Figure 3A. Similar analysis shows that muscle expression does not correlate with ugcaug occurrences in the upstream intron, whereas in the downstream intron the magnitude of correlation decreases with distance from the exon as demonstrated by the increasing P-values (Figure 3B). These dependencies are further corroborated by strong over-representation of this motif in proximal downstream introns in human and other vertebrates (Figure 3C). Indeed, almost half of the muscle-enriched exons in all four datasets (23/56 in human, 21/54 in mouse, 20/43 in chicken and 19/36 in frog) possessed at least one ugcaug motif in the first 200 nt of the downstream proximal intron. Together, these results strongly support the hypothesis that ugcaug is potentially an important regulatory element for muscle-specific alternative splicing, as predicted by the correlation with expression analysis.

Analysis of upstream and downstream intron sequences

We extended the above analysis to identify additional muscle-specific cis-elements in upstream and downstream intron sequences. In the downstream 200 nt sequence, we searched all possible hexamers and identified a total of 35 hexamers that were significantly correlated with expression (P⩽0.05) and also over-represented in the human dataset (P ⩽ 0.05, q⩽0.2); nine of these were also over-represented in at least one other species (Table 1A and Supplementary Table 5A). Several of these elements have been previously characterized experimentally as regulators of muscle-specific splicing. These motifs fell into three distinct classes: (i) the Fox1/2-binding motif ugcaug (P = 6.8E−05) and two closely related hexamers (gcaugg, uuugca); of note, in all four species the majority (58–76%) of GCAUG motifs in the D200 region occurred in the context of the full UGCAUG hexamer. (ii) ug-rich elements gugugu and uguguc (correlation P-values = 0.032 and 0.005, respectively), that resemble binding sites for the CELF family of splicing factors; and (iii) the novel motif acuaac (P = 0.0006) and related hexamers cuaacc (P = 0.004) and cacuaa (P = 0.04). The latter class is similar to the uacuaac element noted in a recent study of a small group of muscle-specific exons in mouse (11). The distribution of these elements in flanking introns of exons in the human, mouse, chicken and frog datasets is shown in Figure 4A. Importantly, this analysis revealed that ugcaug was the most over-represented hexamer in all four datasets, and both gugugu and acuaac were also consistently in the top ∼1% of the most over-represented hexamers in these species.

For upstream intron sequences (200 nt), we found a total of 27 hexamers that were significantly correlated with expression and also over-represented in human, of which three were over-represented in at least one other species (Table 1B and Supplementary Table 5B). Many such elements are strongly pyrimidine-rich, characteristic of binding sites for PTB protein, an inhibitor of splicing for many alternative exons (35). In all four species, the muscle-enriched datasets showed strong over-representation of the reported PTB-binding sites, cucucu and ucuu, in the proximal upstream intron (Figure 4). cucucu was concentrated mainly in the U200 region. ucuu was focused even more tightly in the U100 region (Figure 4B), where it was consistently among the top five over-represented tetramers in all four species. Lesser over-representation of ucuu motifs over a broad area of downstream intron sequences was also noted, perhaps consistent with previous findings that optimal splicing repression by PTB requires binding sites both upstream and downstream of the regulated exon (52,53).

Many of the remaining significant hexamers, both for upstream and downstream introns, have low similarity to the previously discovered elements. Although these may represent novel elements, given that splicing elements are often degenerate, they can also be specific examples of known degenerate motifs. Our analysis using degenerate motifs presented below suggests that the latter possibility is more likely. Finally, for some of the major splicing regulatory elements described above, we observed that the profiles of positional over-representation have been conserved through vertebrate evolution: mouse, chicken and frog. This is displayed in Figure 4 for Fox, CELF and PTB-binding sites. Such strong positional conservation of motifs lends additional support to our findings using correlation with expression.

The Fox-binding site, ugcaug, was previously shown to be over-represented downstream of brain-enriched alternative exons (7,8), raising the possibility that the brain- and muscle-specific alternative splicing programs might exhibit functional similarities by sharing related components of the splicing machinery. To determine which of the candidate muscle cis-regulatory elements might be shared with brain-specific alternative splicing, and which are unique to the muscle program, we compared the frequency of several key cis-regulatory motifs in muscle (this study) and brain (8) datasets. Two elements, acuaac and ugugug, were clearly muscle specific since their frequencies were consistently higher in the D200–D300 region adjacent to muscle-enriched exons compared with the intronic region downstream of brain-enriched exons (Supplementary Figure 1; positive contrast scores). These motifs were also not over-represented in brain relative to control exons (data not shown). In contrast, the motifs ugcaug and cucucu occurred at even higher frequencies in the proximal introns of the brain-enriched dataset than they did in the muscle-enriched dataset (Supplementary Figure 1; negative contrast scores for ugcaug in the D200–D400 region, and for cucucu in the U100 region). Essentially equivalent distribution patterns were observed in the mouse, chicken and frog datasets (data not shown). These results strongly suggest that tissue-specific alternative splicing programs may utilize a combination of unique and shared cis-regulatory motifs that will require much additional analysis in the future.

Motifs identified via PWM analysis are consistent with word analyses

Because many splicing factors bind degenerate oligonucleotide sequences in RNA, we extended our analyses to include degenerate motifs through the use of PWMs (3,54). PWMs are probabilistic representations of degenerate binding sites. Over-represented PWMs in introns of 56 muscle-specific exons in the human dataset were obtained using the DME algorithm (44,45). We scanned multiple parameter settings of DME in order to obtain a large number of PWMs and reduce bias from DME. To identify the functional PWMs, we assessed their correlation with muscle expression using linear splines (19,22). Linear splines are among the simplest non-linear variants of linear models. In contrast to many other approaches, they facilitate adaptively learning the cutoff scores of PWMs that discriminate true targets from false targets of trans-factors. Previous regression approaches have used either maximum score of the PWM (22) or a global average of PWM scores for all potential binding sites (20) on an input sequence as the predictor variable. However, realistically, a small number of sites, sometimes >1, are bound by the corresponding trans-factor. Here we overcame this limitation by including both strength of PWM and the number of putative binding sites in our linear splines approach (see Methods section).

Degenerate 6-nt and 4-nt sequences that were over-represented in the proximal downstream intron sequence are shown in Table 2A and B and Supplementary Table 6A and B. Notably, all of the top 10 over-expressed PWM hexamer motifs in the D200 region are consistent with the major over-expressed unique motifs identified above. Among these, the six most statistically significant motifs represent close matches to the Fox-binding site, ugcaug; two (nhcuaa and hcuaan) are very similar to the novel acuaac element; and the remaining motifs (sukugs and cugysr) resemble ug-rich-binding site for CELF proteins. Analysis of over-expressed 4-mers in the D50 region revealed that the top-scoring motif is ugcm. While this motif is included in the Fox recognition sequence, other considerations suggest that it would not be sufficient for Fox binding. Instead, ugc likely represents the cug-rich sequences, characteristic of some CELF-binding sites.

Table 2.

Significant position weight matrices (PWMs) identified by correlation analysis using linear splines.

Open in a new tab

Panel A: Top 10 PWMs of width 6 nt in proximal downstream intron sequences (length = 200 nt).

Panel B: All significant PWMs of width 4 nt in proximal downstream intron sequences (length = 50 nt).

Panel C: Top 10 PWMs of width 6 nt in upstream intron sequences (length = 200 nt). Complete list of significant PWMs in downstream and upstream 200 nt introns is shown in Supplementary Table 6.

In the U200 region, all of the statistically over-represented motifs were quite pyrimidine-rich relative to the control group (Table 1B and 2C). Further investigation will be required to determine whether these elements are primarily bound by the PTB protein or by additional splicing factor(s). All remaining PWMs that have high significance in upstream and downstream introns exhibit at least partial similarity to the above three elements. For PWMs with only partial similarity, the similarity is observed either at the 5′ or at the 3′ end of the motif, indicating that the remainder of the motif most probably represents the flanking region. For example, for the PWM rrwgca, the last four bases match the 5′ end of the ugcaug, and hence, the first two bases are presumably the flanking region of this putative Fox-binding motif.

Furthermore, in contrast to previous work (22), the new formulation of linear splines used here allowed us to obtain not only the potential target exons, but also the binding sites of the above splicing factors (see Methods section). The results for a representative set of motifs are summarized in the Supplementary Table 7. For the putative Fox-binding motif wgcauk, we find ugcaug as the most frequently occurring oligonucleotide sequence, as expected of Fox-binding sites. We have observed similar accuracy in binding site prediction in the context of transcriptional regulation (Das,D., unpublished data). Interestingly, we notice that not all possible combinations of nucleotides of a degenerate PWM are realized in the set of 56 muscle-specific exons. For example, for the candidate CELF-binding motif, cugysr, only cuguga is predicted as the binding site. These are consistent with the previous observations made in the context of transcriptional regulation (55,56).