Computational Identification of Tissue-Specific Splicing Regulatory Elements in Human Genes from RNA-Seq Data

Alternative splicing is a vital process for regulating gene expression and promoting proteomic diversity. It plays a key role in tissue-specific expressed genes. This specificity is mainly regulated by splicing factors that bind to specific sequences called splicing regulatory elements (SREs). Here, we report a genome-wide analysis to study alternative splicing on multiple tissues, including brain, heart, liver, and muscle. We propose a pipeline to identify differential exons across tissues and hence tissue-specific SREs. In our pipeline, we utilize the DEXSeq package along with our previously reported algorithms. Utilizing the publicly available RNA-Seq data set from the Human BodyMap project, we identified 28,100 differentially used exons across the four tissues. We identified tissue-specific exonic splicing enhancers that overlap with various previously published experimental and computational databases. A complicated exonic enhancer regulatory network was revealed, where multiple exonic enhancers were found across multiple tissues while some were found only in specific tissues. Putative combinatorial exonic enhancers and silencers were discovered as well, which may be responsible for exon inclusion or exclusion across tissues. Some of the exonic enhancers are found to be co-occurring with multiple exonic silencers and vice versa, which demonstrates a complicated relationship between tissue-specific exonic enhancers and silencers.


Introduction
Alternative splicing (AS) is an essential cellular process in eukaryotes that pre-mRNA usually undergoes to produce multiple mRNA isoforms of the same gene with likely different functions [1]. In a typical splicing process, introns within the pre-mRNA are removed, and the exons are joined together to form the mature mRNA [2,3]. AS supports the joining of different combinations of exons; said another way, what sequence constitutes an exon is readily redefined. As a consequence, different proteins are produced from the same gene. Recent studies show that AS occurs in more than 95% of human genes [3,4]. AS is regulated by specific proteins, called splicing factors. They bind to certain short sequences on the pre-mRNA, called splicing regulatory elements (SREs). Identifying these SREs and their combinatorial effects are crucial to understanding AS. SREs are classified as enhancers if they promote exon inclusion and silencers if they inhibit exon inclusion [3,5,6]. Accurate splicing is crucial, as it is believed that mutations either in the core splicing signals or in the SREs contribute to approximately 50% of human genetic diseases [4,7,8].
Alternative splicing plays a key role in generating tissue-specific proteins [3,9]. Tissue-specific alternative splicing is regulated by a combination of tissue-specific and ubiquitously expressed RNA-binding factors [10]. They interact with the splicing regulatory elements to affect the spliceosome assembly (splicing machinery) and consequently the transcribed isoforms. There are several splicing factors that activate or repress splicing in different contexts [10]. Skipping exons is one of the notable alternative splicing events between different tissues that results in different proteins. Wang, et al. [10] suggested that having this switch-like regulation between tissues requires additional splicing regulatory elements to be present. Most work in mammalian systems revealed that AS decisions are often made by a combinatorial action of general and tissue-specific regulators [11]. Even simple tissue-specific decisions can involve additional layers of complexity, where regulatory elements cooperate or compete with each other [12].
Several studies have identified tissue-specific regulatory elements. Mouse RNA-Seq data for three tissues (brain, liver and skeletal muscle) were utilized to calculate the expression level of each isoform of genes for a set of predefined cassette exons [3]. The authors identified 456 putative enhancers and silencers. Among these, 45 were common to all tissues.
Kim et al. [13] utilized a distribution-based quantitative association rule mining to find exonic and intronic sequence motifs in 10 mouse tissues. Combinatorial cis-regulatory motifs were also discovered. The author identified statistically significant associations between sequence motifs and tissue-specific exon skipping rates. Ninety-seven interesting association rules were identified, of which three contain multiple 7-mers.
Focusing on human tissues, a varying effect regression model on splicing elements (VERSE) was developed to predict genome-wide intronic SREs [14], where RNA-Seq data for 16 human tissues was used. Approximately half of the SREs (55.68%) were found to be significant only in one tissue.
Castle, et al. [15] designed microarrays monitoring exons and exon-exon junctions in 17,939 human genes. For the regulatory elements, they extracted sequences in eight neighborhoods around regulated exons where 135 motifs were identified.
Brudno, et al. [16] identified intronic regulatory elements that are brain-specific, while, in [17], the identified regulatory elements were muscle-specific. In [8], a probabilistic approach was utilized and several intronic regulatory elements in different human tissues were identified.
Wang, et al. [18] developed a linear regression model to estimate the effect of various splicing factors on exon inclusion between two tissues, and, hence, the binding sites of these splicing factors are predicted. They applied their model on data from liver and heart tissues and predicted 15 motifs that contribute to exon skipping events. The work was extended to 11 human tissues in [9].
Wang, et al. [10], analyzed 10 human tissues. A high frequency of tissue-specific regulation was observed for each of various alternative splicing event types, including over 60% of the analyzed skipped exons.
Ke and Chasin [19] utilized a hypergeometric test to discover sequence pairs that are overrepresented in intronic regions flanking human exons. They identified more than 60,000 5-mer sequence pairs with p 10 −4 . They showed that some pairs are associated with tissuespecific genes.
Taking into consideration the combinatorial effect of SREs, a biophysical principals based regression model for the regulation of AS was developed in [20]. It captures both the main effects of individual SREs and the combinatorial effects of SRE pairs. Overall, 619 different SREs and 196 SRE pairs were detected from different tissues. Their model was limited to the interaction of at most two SREs.
As discussed above, most of the human tissue-specific studies focus on identifying individual SREs in intronic regions [8,14,16,17]. Wang, et al. [9,18] identified SREs in both intronic and exonic regions. However, the focus was only on individual regulatory elements as well.
There are other studies that identify individual SREs in both exonic and intronic regions [21][22][23][24][25][26][27][28], which we use for comparing our results. However, the focus in these studies was not on tissue-specific SREs. For more detailed review about these approaches, please refer to our previous paper [21].
The tissue-specific studies that focus on the combinatorial SRE effect were [19,20]. In [20], out of 196 identified SRE pairs, only two pairs have both SREs in the exonic region. In [19], the identified pairs reside in the intronic regions.
In this work, we performed genome-wide analysis to study alternative splicing on multiple tissues (brain, heart, liver, and muscle). The RNA-Seq data set from the Human BodyMap project [29] was utilized. We used DEXSeq [30] to identify differential exons across tissues. Then, we applied our algorithms, GenSRE [21] and CoSREM [31], to identify both individual and combinatorial exonic regulatory elements responsible for exons that exist in one tissue but not in other tissues. Putative tissue-specific exonic enhancers were discovered and a complicated exonic enhancer regulatory network was revealed. Multiple exonic enhancers were found across multiple tissues, while some were found only in specific tissues. Putative combinatorial exonic enhancers and silencers were discovered as well that may be responsible for exon inclusion or exclusion across tissues. This is, to our knowledge, the first analysis to focus only on discovering exonic regulatory elements (individual and combinatorial) across tissues.

Materials and Methods Data
RNA-Seq data from the Human BodyMap 2.0 project is utilized [29]. This data originates from 16 different human tissues. It contains 50 bp paired-end reads, 75 bp single-end reads, and 100 bp single-end reads. We focus on four tissues, namely, brain, heart, liver, and muscle. All the samples provided for tissues of interest have been utilized.
For comparing our results with previously published results, several databases are utilized. SpliceAid-F [22] is a recent comprehensive database that includes all the experimentally verified splicing factors and their binding sites. It contains 71 splicing factors and 655 binding sites for human. We also used AEdb [23], which is a database for alternative exons and their properties from various species; it is the manually curated component of the Alternative Splicing Database (ASD). The exon data in AEdb have been experimentally verified.
In addition, we compared our putative exonic splicing enhancers (ESEs) with five other computational data sets. The RESCUE-ESE [24] data set contains 238 6-mers for human exons. Another data set is PESE [25], where 2096 8-mers were identified. The third data set is from [26] and contains 4-and 5-mers as potential ESEs. In the data set from [27], the authors concentrated on 5-mer putative ESEs.
For exonic splicing silencers (ESSs), we compared our results with FAS-ESS [28], and PESS [25]. The FAS-ESS data set contains 130 10-mer sequences that were identified utilizing the mini-gene approach. PESS is another data set where the authors compared the frequencies of 8-mers (allowing one mismatch) in constitutively spliced non-coding exons with those in pseudo-exons and the 5 0 untranslated regions (UTRs) of intronless genes.

Overview of our proposed pipeline
In our pipeline, different tools are utilized to identify tissue-specific exonic regulatory elements. The first stage (Fig 1) is to identify differential exons. To do that, we utilized DEXSeq [30]. DEXSeq identifies exons that are differentially used between two tissues. We compare each tissue of interest with the other three tissues. The output of this stage is three sets of differentially used exons in the tissue of interest but not in the remaining tissues. We then apply GenSRE [21] separately on the exon sets to identify exonic enhancers. We identify tissue-specific exonic enhancers by determining the common ESEs across the sets. We also apply CoS-REM [31] to identify co-occuring exonic enhancers and silencers that may be responsible for exon inclusion or exclusion across tissues.

Identifying differential exons across tissues
Our main goal is to identify splicing regulatory elements that are tissue-specific. To do that, we utilized DEXSeq [30] to identify exons that are differentially used between two tissues. Differential exon usage analysis aims at identifying the changes in relative usage of exons caused by a certain condition [32].
Let R g = {t 1 , t 2 , . . ., t n } be the set of all the transcripts from a specific gene g. Let X g = {e 1 , e 2 , . . ., e m } be the set of all exons that constitute the transcripts of gene g. Let T(e, g) be the set of all the transcripts from R g that contain an exon e. Exon usage U(e, g) of exon e is defined to be Uðe; gÞ ¼ Tðe; gÞ n In DEXSeq, generalized linear models are utilized to model read counts, and the χ 2 likelihood ratio test is then used to get an analysis of deviance p-value. The null hypothesis in this test is that none of the conditions influence exon usage. Rejecting the null hypothesis indicates that the count of sequencing reads that map to the exon under the test differs significantly between the different conditions. One of the advantages of this model is accounting for biological variability when the data has replicates for different conditions in contrast to other methods [33].
Before conducting differential exon usage analysis, flattening gene models and counting the reads were performed. Flattening gene models means aligning the sequencing reads to a reference genome and accumulating all the reads for each exon in each tissue from all the transcripts that contain this exon. As some of the transcripts may contain only a part of an exon, the exon counting bins term is used to refer to an exon or a part of an exon [32,34]. It is notable that having exons with differential usage does not mean that their corresponding transcripts are differentially expressed. The output of DEXSeq is a table that contains the differential counting bins, their genes, their read counts in both tissues, and the p-values to determine significance.
For each pair of tissues, we applied DEXSeq. We identify differentially used exons by choosing the exons with p-value 0.05. We further filter these exons by choosing the exons with log 2 fold change ! 2 or −2. To focus on exons that exist in one tissue but not in the other, we only used exons that have reads in one tissue and no reads at all in the other tissue. The output of this stage is a set of exons in each tissue that are differentially used and present in one tissue but not the other.
These exons were then retrieved from the ENCODE project [35]. The February, 2009, human genome assembly (GRCh37/hg19) was used. The complete exonic sequences were retrieved. The 200 intronic nucleotides upstream and the 200 intronic nucleotides downstream of each exon were also retrieved.

Identifying tissue-specific enhancers
We utilized a formalism based on de Bruijn graphs that we previously reported [21]. It combines genomic structure, word count enrichment analysis, and experimental evidence to identify enhancers found in exons [21]. In this model, a six-dimensional de Bruijn graph is constructed G = (V,E) over the DNA alphabet σ = {A,C,G,T}, and each vertex is associated with its rank based on LEIsc (log of the enrichment index, scaled) scores from Ke et al. [36]. Then, a subset U & V is chosen to be the 400 vertices with the highest LEIsc values. An SRE graph G U is then constructed and its weakly connected components are calculated. The algorithm Gen-SRE is then applied to each weakly connected component to determine a set of potential enhancers. Finally, these sequences are submitted to word count enrichment analysis.

CoSREM
CoSREM (Combinatorial SRE Miner) is an algorithm for discovering combinatorial SREs that we reported in [31]. CoSREM is a two-level graph mining algorithm that is applied to the SRE graphs from [21] to identify co-occurring sets of SREs. The focus is on identifying sets of exonic splicing regulatory elements whether they are enhancers or silencers. Experimental evidence is incorporated through the SRE graphs to increase the accuracy of the results. The identified SREs do not have a predefined length, and the algorithm is not limited to identifying only SRE pairs.

GO enrichment analysis
We utilized the command-line version of Ontologizer [37], with the goal of determining the enriched GO annotations of the genes that contain the identified enhancers that appears only in one tissue.
For each tissue, the genes, whose exons were identified as differential exons, are utilized as a background data set (population set). For each exonic splicing enhancer in a specific tissue, the exon data set is searched to allocate each splicing enhancer, and the corresponding gene set is identified to form the study set. GO annotation files gene_ontology_edit.obo and gene_association.goa_human were downloaded. GO enrichment analysis is performed using the Topology-Elim algorithm. The Westfall-Young Single Step multiple testing correction procedure is then applied.
We are interested in the biological process annotations. Therefore, we choose the biological process category with the minimum adjusted p-value, where we consider only terms with p 0.05 to be significant.

Differentially used exons between tissue pairs
In this study, we analyze four tissues from the RNA-Seq data of the Human BodyMap project [29]. DEXSeq [30] is utilized to identify differentially used exons between pairs of tissues. Table 1 illustrates a part of the DEXSeq output for each pairwise comparison.
For each tissue pair comparison, DEXSeq produces a list of exons that are differentially used in one tissue but not in the other tissue. These exons may be whole exons or parts of exons. Being conservative, we define tissue-unique exons as differentially used exons with p 0.05, log 2 fold change ! 2 or −2, and exons that have reads in one tissue and zero reads in the other tissue. We preferred to have a strict set of exons that exist in one tissue against the other. Table 2 illustrates the number of exons resulting from each pairwise comparison. The complete list is listed in S1 File.For example, applying DEXSeq on brain and heart tissues resulted in 1423 exons that are brain-specific. In the same manner, 857 exons were found to be heart-specific while having zero reads in brain tissue. It should be noticed that brain tissue has the largest number of tissue-unique exons, which is consistent with results reported in [13] that the brain has a large number of tissue-specific alternative spliced exons. We performed 6 pairwise tissue comparisons. Therefore, we have 3 different sets of exons for each tissue with a total of 12 tissue-unique exon sets for all 4 tissues.

Tissue-specific exonic enhancers
We used a de Bruijn graph based algorithm called GenSRE to identify ESEs in each tissue [21].
Having three exon sets for each tissue, we applied GenSRE on each set separately. Table 3 indicates the number of the identified ESEs in each set. The complete list is listed in S2 File. We focused on identifying putative exonic enhancers that may play a role in the inclusion of these exons within these specific tissues. In GenSRE, the exonic flank size is determined to be 50 nucleotides. Therefore, all exons that don't meet this criterion are excluded. Table 4 contains the number of exons that were utilized. GenSRE identifies variable length SREs in the exonic flanking regions. The identified ESEs lengths range from 6 to 15 nucleotides. To identify tissue-specific ESEs, we extracted the common set of exonic enhancers across the three exon sets for each tissue. This behavior suggests that these ESEs are tissue-specific as they repeatedly appeared in one tissue against all the other tissues. Table 3 illustrates the number of tissue-specific ESEs for brain, heart, liver, and muscle tissues. S1 Table includes the list of tissue-specific ESEs for all tissues. We also wanted to assess the accuracy of this pipeline in identifying tissue-specific ESEs. Therefore, for each tissue, we identified the set of exons that is present in only that tissue. This is done by finding the intersection between the different exon sets of each tissue. For example, in the case of the brain tissue, each exon set represents exons that are present in the brain and not in the heart, liver, and muscle tissues, respectively. The common set of exons between these three sets represents brain-unique exons (Fig 2).
Then, we searched for tissue-specific ESEs that we previously identified in these unique sets of exons, to see if they appear in the exonic flanking regions or not. Table 5 illustrates the number of ESEs that were found for each tissue. For example, out of 205 exonic enhancers that we Each circle represents the number of brain-specific exons that resulted from brain pairwise comparisons with other tissues (heart, liver, and muscle). The intersection represents brainunique exons against all other tissues.  identify as brain-specific ESEs, we found 173 ESEs in the brain-unique exons, which represents about 85% of all the identified ESEs. It is notable that, for instance, as indicated in Table 5 for the heart tissue, although there are 85 ESEs identified as heart-unique, the number of heartspecific exons is only 23 exons. In these 23 exons, we found 19 ESEs, which is about 83% coverage for the exons. This result suggests that our approach can be utilized to identify tissue-specific ESEs.

Comparison with SRE databases
We compared our results with databases from [22][23][24][25][26][27]. We extracted only exonic enhancers and silencers to compare our results with. We applied the same approach we used in [21]. We first compared our results with exonic binding sites from SpliceAid-F [22]. SpliceAid-F contains 330 different binding sited for human, 59 of them are exonic enhancers. Since our predicted ESEs are of variable length, as are SpliceAid-F binding sites, we calculated the overlap between the two sets by finding whether each sequence in the first list is entirely contained in at least one sequence in the second list or vice versa.
Another database is AEdb [23], where we considered only the 64 human enhancers.
In addition, we compared our ESE list with three other computational data sets, such as the RESCUE-ESE [24] data set, the PESE [25] data set, and the data set from [26]. The data set from [26] contains only 4-and 5-mers as potential enhancers. As a result, we could only test if sequences in our list include any of its sequences. This also applies to the data set from [27]. Table 6 summarizes the overlapping results. Overall, about 61.8% of the predicted ESEs can be matched to one of the previously published databases.

Tissue-specific ESE regulatory network
To understand the complex relationship among the identified exonic enhancers found in multiple tissues, we constructed a regulatory network (Fig 3). It is a bipartite graph with two types of nodes: the circular nodes represent ESEs and the rectangular ones represent tissues. The size of any node type is proportional to its degree (number of incident edges). It is clear that there are ESEs that are specific to one tissue and others that may regulate more than one. That is in accordance with the suggestion in [11] that AS decisions are often made by a combinatorial action of general and tissue-specific regulators.
Focusing on the ESEs that are involved in multiple tissues (Fig 4), we noticed a hierarchical relationship, where some ESEs regulate two tissues (10% of all ESEs). A smaller number regulate 3 tissues (1%), and only one ESE is found in all four tissues. The other 89% of the SREs are tissue-specific where brain-specific exonic enhancers represent 64% of all the identified ESEs. These results are consistent with the conclusions from [13], that brain tissue exhibits a very large number of tissue-specific SREs and a limited number of general ones. However, in [13], the focus was on intronic SREs.

Tissue-specific combinatorial SREs
Identifying individual cis-regulatory elements does not suffice to explain tissue-specific or condition-specific AS. For instance, in the case of exon skipping events, if an exon has both enhancer and silencer elements in proximity and in case of having an SR splicing factor with great affinity (SR factors are proteins that bind to enhancers and play various roles in spliceosome assembly [38]), the SR protein will bind to the enhancer and stimulate exon inclusion. This is through recruiting other spliceosome proteins, such as U1 and U2, to the core splicing signals. Consequently, the spliceosome machinery is assembled, and the exon is included.
On the other hand, if an inhibitory splicing factor such as hnRNP A1, which acts as a splicing repressor, is also present, it may inhibit the exon inclusion by binding to the silencer sequence and recruiting the binding of other inhibitory factors. These factors extend to the exon boundary and prohibit the binding of the SR protein. As a result, the exon will be skipped [11,38]. We utilized our graph mining algorithm CoSREM [31] to identify co-occurring exonic enhancers (ESEs) and silencers (ESSs) that may cause exon inclusion in one tissue and its exclusion in another tissue. To do that, the identified exon set for each tissue was utilized. Table 7 illustrates the number of exons used for each tissue and the number of identified combinatorial exonic enhancers and silencers.
The actual combinatorial SRE sets are given in S3 File. We notice that these SRE sets appear in most of the specified exons. For example, in the brain tissue, the 366 combinatorial SRE sets appeared in 8753 out of 8858 exons, which may explain the inclusion of these exons in the brain tissue and their exclusion in the other tissues. Although the number of combinatorial enhancers and silencers in the brain tissue was quite large, the number of unique enhancers and silencers was surprisingly small. There are 63 unique ESSs and 30 unique ESEs, of which 26 are identified by GenSRE in the previous section as brain-specific ESEs. Table 8 illustrates the number of ESEs that were identified and verified previously. Table 9 illustrates the number of ESEs and ESSs identified as a part of combinatorial SREs for all the tissues.  To understand the relationship between these ESEs and ESSs, we constructed a regulatory network for exonic enhancers and silencers in the brain tissue (Fig 5). The red nodes represent ESEs, and the blue ones represent ESSs. The node size is proportional to node degree. The network illustrates the many-to-many relationship between the exonic enhancers and silencers. In other words, one ESE can co-occur with multiple ESSs and vice versa. The numbers in parentheses are the numbers of ESEs discovered as a part of a combinatorial set.

Discussion
We propose a graph-based computational approach to identify exonic SREs that contribute to exon skipping events across tissues. We utilized DEXSeq [30] to identify differential exons across tissues, and then we applied GenSRE [21] and CoSREM [31] algorithms to identify both individual and combinatorial SREs that are tissue-specific. We utilized Ontologizer [37] to assess the significance of our predicted tissue-specific ESEs and whether they are involved in tissue-specific biological processes. Therefore, for each tissue, we determined the enriched GO annotations of the genes that contain the identified ESEs. We focused on identifying significant biological processes with adjusted p-value 0.05. We identified several brain-related processes, as illustrated in Table 11. We also identified some heart-related processes such as "regulation of cardiac muscle contraction by regulation of the release of sequestered calcium ion" with p-value 0.00295. Other biological processes were identified as liver and musclerelated such as "digestive tract morphogenesis", and "regulation of muscle system process" with p-values 0.00202, and 0.00253, respectively. S4 File includes the complete list of the enriched GO annotations for the four tissues of interest. We applied the same analysis to identify the enriched GO annotations of the genes that contain combinatorial SREs. Various functionally related biological processes to the tissue of interest are recognized. The complete list of the enriched GO annotations of the four tissues is included in S5 File. They resulted of GO enrichment analysis of gene sets that contain putative brain-specific ESEs identified by our approach. Population total is the total number of genes in the population set. In this case, we have 5096 genes in the population set. Pop. count is The number of genes in the population set that are annotated to the GO term in question. Study total is the number of genes in the study set. Study count is the number of genes in the study set that are annotated to the GO term in question.
We also identified tissue-specific combinatorial SREs. These are sets of co-occuring exonic enhancer and silencer elements in each tissue. As our focus here is exons that are differentially used in one tissue with respect to other tissues, identified combinatorial SREs may play a role in exon inclusion in one tissue and its exclusion in another tissue. Mayeda et al. [39] showed in vitro that having different ratios of SF2/ASF to hnRNP A1 splicing factors promotes exon skipping or inclusion by binding to different enhancers or silencers. In other words, it was shown that the ratio of SF2/ASF to hnRNP Al can affect whether the internal exon is included or excluded. An excess of SF2/ASF promotes exon inclusion, while hnRNP Al excess promotes exon exclusion.
Therefore, we investigated this hypothesis by incorporating the splicing factor information. We identified splicing factor proteins (enhancer factors and inhibitor factors) from SpliceAid-F [22] that bind to SREs in our combinatorial SRE sets. We then calculated the expression levels for these splicing factors from the RNA-Seq data (Fig 6). The actual expression levels in FPKM is in S2 Table. Then, for each combinatorial SRE set, we calculated the ratio between the expression levels of its splicing factors, if any. Table 12 provides an example from the combinatorial SREs identified in the brain tissue. We focused on results where the ratio ! 1 in one tissue and 1 in at least another tissue, which suggests that these splicing factors and their binding sites may play a role in regulating exon inclusion or exclusion between the tissues.
An interesting result involved the splicing factors FMRP and hnRnpLL where the FMRP to hnRNPLL expression level ratio was (% 1.83) in brain tissue while it was (% 0.52) in muscle tissue as illustrated in Table 12. More interestingly, These splicing factors were associated with multiple SRE sets. We identified multiple co-occuring exonic enhancer and silencer elements that all are potential binding sites to these splicing factors in the brain tissue. For example, one ESS identified as a potential binding site to hnRNPLL protein is CAGCCA. It co-occurs with different exonic enhancer elements (AAGAGA,AAGGAA,CAAGGA,GAGAGC,GTGGAG, AGAGGA). All ESEs were identified as binding sites to FMRP splicing factor. Other ESSs elements were identified as well as binding sites to hnRNPLL such as CCACCA, CCAGCA. All the identified binding sites include CA dinucleotide repeats, which is known to be preferentially recognized by hnRNPLL [40]. This suggests the hypothesis that these two splicing factors may have an antagonistic behavior that results in some exons being included in the brain tissue and excluded in the muscle tissue. Although this hypothesis needs further experimental validation, our approach can highlight interesting results for more experimental testing.

Conclusion
In summery, we report a genome-wide analysis to study alternative splicing on multiple tissues. That includes brain, heart, liver, and muscle tissues. Our proposed pipeline identifies differentially used exons across tissues and hence tissue-specific exonic SREs. Both individual and combinatorial SREs are identified which may be the reason for exon inclusion or exclusion across different tissues. Our approach helps understand the relationship between different types of exonic SREs. It yields interesting results that can open new directions to study alternative splicing and how it may influence tissue specificity.
Supporting Information S1 File. It contains 6 tables. Each table contains a list of genomic coordinates of significant differential exons for each pairwise comparison between tissue pairs. (XLSX) The ratio columns contain the ration between the expression level of the enhancer factor and the inhibitor factor in the specified tissue. The cells with 'non' indicates that we could not identify an associated splicing factor and hence no ratios are provided.