Large-scale identification and characterization of alternative splicing variants of human gene transcripts using 56,419 completely sequenced and manually annotated full-length cDNAs.

We report the first genome-wide identification and characterization of alternative splicing in human gene transcripts based on analysis of the full-length cDNAs. Applying both manual and computational analyses for 56,419 completely sequenced and precisely annotated full-length cDNAs selected for the H-Invitational human transcriptome annotation meetings, we identified 6877 alternative splicing genes with 18 297 different alternative splicing variants. A total of 37,670 exons were involved in these alternative splicing events. The encoded protein sequences were affected in 6005 of the 6877 genes. Notably, alternative splicing affected protein motifs in 3015 genes, subcellular localizations in 2982 genes and transmembrane domains in 1348 genes. We also identified interesting patterns of alternative splicing, in which two distinct genes seemed to be bridged, nested or having overlapping protein coding sequences (CDSs) of different reading frames (multiple CDS). In these cases, completely unrelated proteins are encoded by a single locus. Genome-wide annotations of alternative splicing, relying on full-length cDNAs, should lay firm groundwork for exploring in detail the diversification of protein function, which is mediated by the fast expanding universe of alternative splicing variants.

INTRODUCTION

Alternative splicing is a phenomenon in which different combinations of exons are spliced to produce distinct transcripts (1). Especially in higher eukaryotes, alternative splicing is frequently used as versatile means of producing diverse transcripts from a single gene locus. The alterations of the exons may cause changes in the encoded amino acid sequences and, at least in some cases, produce functionally divergent proteins by modifying, for example, the binding site of a growth factor receptor or an activation site of transcription factor (2,3). The most striking example of this is the Drosophila DSCAM gene, which is an axon guidance receptor gene. This gene consists of 17 exons, with 12, 48, 44 and 2 mutually exclusive alternative splicing exons for exons 4, 6, 9 and 17, respectively. Thus, even if not all combinations of these exons are allowed, this gene can encode thousands of protein products, which should enable functional diversification of the protein that could in turn assure precise axonal trajectory (4,5). Since the initial draft sequence of the human genome revealed that there seem to be an unexpectedly small number of genes embedded in the human genome (6), it has been hypothesized that alternative splicing is one of the most significant processes giving rise to the functional complexity of the human genome and that it might be indispensable for generating highly complex organisms such as humans (7).

However, despite the growing interest on the impact of alternative splicing in various aspects of the biological processes, our understanding of alternative splicing is still very primitive and its mechanisms of control are mostly unknown (8). In order to advance our understanding of the biological significance of alternative splicing in humans, it is essential to identify and characterize the genes that are subject to alternative splicing and which splicing patterns are used in what context in a genome-wide manner.

For this reason, large-scale attempts to identify alternative splicing have been initiated by several groups [i.e. (9)], mainly using bioinformatics analysis of partially sequenced cDNAs (ESTs). The EST sequences are clustered and compared to evaluate if the differences in their sequences might have been resulted from alternative splicing. So far, millions of human ESTs have been analyzed, and the newly identified alternative splicing variants are presented in several databases, such as AceView (130 576 alternative splicing variants from 19 557 genes; ; Y. Kohara et al., in preparation), ASAP [30 793 alternative splicing variants from 7991 genes; ; (10)] and ASD [73 340 alternative splicing variants from 16 236 genes; ; (11)]. However, these recent high-throughput approaches have limitations. The first is the skewed coverage of the ESTs over the mRNAs. Since the ESTs are generally scarce around the 5′ ends of an mRNA, previous data could have been biased towards the 3′ ends. Secondly, some combinations of the alternative splicing exons may not be allowed. In those cases, analysis of only partially sequenced cDNAs would have no chance to identify mutual dependence of the combinations in alternative splicing.

Full-length cDNAs provide an ideal solution to all of these limitations. Moreover, not only full-length cDNAs can be utilized to extract complete cDNA sequence information, but are also useful as physical reagents indispensable for experimental analysis to examine the functional consequences of the identified alternative splicing. In this study, we used 56 419 cDNA sequences of human genes selected for the H-Invitational human transcriptome annotation meetings (12). These cDNAs were enriched for full-length cDNAs using various methods (13–16). In addition, they were fully sequenced with a sequence reliability higher than 99% [Phred values greater than 30; (17)]. Also, potentially problematic sequences, such as vectors and poly(A) tails were precisely trimmed. Thus, this cDNA collection constitutes an outstanding resource for comprehensive studies of alternative splicing. Out of the 56 419 cDNAs, 55 036 were successfully mapped onto the human genomic sequence (UCSC hg16; ) and clustered into 24 425 loci. Of these, 10 127 loci contained two or more cDNAs and both manual and computational inspection allowed us to identify 18 297 alternative splicing variants encoded in 6877 loci [Table 1; (12)]. General statistics and a part of the related information have been published as a part of H-Invitational paper (12). In this paper, we describe further detailed features of alternative splicing. Here, we report the large-scale genome-wide identification and analysis of human alternative splicing based on fully sequenced and precisely annotated full-length cDNAs.

Table 1

Statistics of the data processing and of the alternative splicing variants and exons identified

	#Locus	#cDNA	#Total exon	#Alternative exon	#Constitutive exon
H-Invitational total	25 585	56 419	389 895a	44 727	345 168
Successfully mapped	24 425	55 036	389 895	44 727	345 168
≥2 cDNAs per locus	10 127	35 030	331 924	44 727	287 197
Identified alternative splicing	6877	18 297	176 505	37 670	138 835
5′ end alternative splicing	4568	7494	18 297	7494	10 803
Internal alternative splicing	5565	11 156	139 911	25 236	114 675
3′ end alternative splicing	2933	4940	18 297	4940	13 357
5′-UTR alternative splicing	3216	4750	18 262	6398	11 864
CDS alternative splicing	6005	13 409	148 242	28 728	119 514
3′-UTR alternative splicing	797	1034	5877	1401	4476

aUnmapped cDNAs' exons could not be counted.

MATERIALS AND METHODS

Dataset description

In the present study, the set of the 56 419 cDNAs, selected for the H-Invitational human transcriptome annotation meetings was used. Contributors and attributes of each of the cDNA data subsets are described in Supplementary Table 1, the reference (12) and further references therein.

Computational procedures to identify and characterize the alternative splicing variants

We mapped the full-length cDNAs to the human genome (UCSC hg16; ). Alignments were generated using EST2GENOME (); alignments having at least 95% identity and 90% coverage were selected. The cDNAs mapped on the same genomic region (at least one base overlap) were clustered and regarded as a putative ‘locus’. For further details of the mapping and clustering, see Ref. (12).

Since this study focuses on complete transcript variants, we conducted a sequence inspection of the mapped cDNAs to identify and discard 5′ end-truncated cDNAs from the dataset. To this end, we excluded cDNAs whose 5′ ends had been located inside the second or later exons of any other cDNAs with compatible exon structure in the same locus. We accepted the cDNAs whose 5′ ends were located inside of the first exons and considered as variations in the exact transcriptional start sites. We also assumed that those cDNAs whose 5′ ends were located outside of the exonic regions of any other clones could not be truncated forms of any known types of transcripts, at least. This assumption is based on the fact that the combination of multiple errors, for example, truncation followed by erroneous oligo-capping that occurred on an immature form, would be required to erroneously generate such cDNAs [for further detailed discussion of this subject, see Ref. (18)]. Using similar concepts, we examined the completeness of the 3′ end using the same procedure and similarly removed all possible 3′ end-truncated cDNAs. The 3′ ends located inside of the last exons were allowed and considered as alternative polyadenylation sites.

Computational identification of the alternative splicing variants was then performed using the resulting filtered set of full-length cDNAs as follows: (i) The genomic position of each exon–intron boundary was compared with those of the other transcripts belonging to the same locus. For the comparison, a 10 bp allowance was made; (ii) If a cDNA had a part of the exonic sequence in the first/last exon inside confirmed intronic regions of the other cDNAs, it was regarded as being a ‘5′/3′ end’ alternative splicing variants. (iii) If a cDNA had a part of an internal exonic sequence inside a confirmed intronic region of other cDNAs, it was recognized as being an ‘internal’ alternative splicing variants (for a schematic representation, see Figure 1).

Figure 1

Schematic representation of the identification of the alternative splicing. Essentially, the illustrated patterns of the exon pairs were searched for and selected as alternative splicing exons in both computational and manual annotations.

In order to evaluate and characterize the outcomes of the identified alternative splicing events on the encoded protein sequences, we used the information of the ORFs (i.e. positions and reading frames) annotated during the H-Invitational meetings. Differences in the length of the ORFs were evaluated in a pair-wise manner between all alternative splicing variants within the locus and the average ORF length difference was calculated for each locus. Possible targets for nonsense-mediated decay were selected as variants in which the stop codon mapped more than 50 bp upstream of the last exon junctions. For the detection of Alu-like elements, RepeatMasker was run with default settings and for the detection of exonic splice enhancers (ESEs) (19), the RESCUE-ESE program was run as described previously (20).

Based on the deduced amino acid sequences for alternative splicing variants, protein motifs and Gene Ontology (GO) terms were predicted using InterProScan () with default parameters. GO terms were automatically added to each of the variants when a protein motif(s) was recognized and could be associated with functional annotation. Enrichments of the motifs and GO terms were statistically evaluated using a hypergeometric distribution by using the following equation: Here, N = 24 425 (number of loci containing successfully mapped cDNAs), n = 6877 (number of loci containing alternative splicing variants), M = 12 764 (number of motif-containing loci), k = 5523 (number of motif-containing loci with alternative splicing variants) in the case of Table 2. Similar calculations were done to evaluate the statistical enrichments for particular motifs and GO terms in Tables 3 and 4.

Table 2

Relation between alternative splicing genes and motifs

	#Motif-related locus	#NOT motif-related locus	Total
Alternative splicing locus	5523a (3015; motif-changed)	1354	6877
NOT alternative splicing locus	7241	10 307	17 548
Total	12 764	11 661	24 425

aP-value < 10−16.

Table 3

Most frequently observed motifs which were affected by alternative splicing variants

InterPro ID	Motifs in alternative splicing locus	Motifs in all locus	Ratio	Significance of enrichment (P-value)	Definition
003598	417	495	0.84	<10−16	Immunoglobulin C-2 type
000005	237	245	0.97	<10−16	Helix–turn–helix, AraC type
000867	73	79	0.92	<10−16	Insulin-like growth factor-binding protein (IGFBP)
000345	114	211	0.54	10−15	Cytochrome c heme-binding site
003962	55	78	0.71	10−15	Fibronectin, type III subdomain
002017	56	88	0.64	10−12	Spectrin repeat
000379	62	103	0.6	10−11	Esterase/lipase/thioesterase
002035	42	60	0.7	10−11	von Willebrand factor, type A
000595	22	25	0.89	10−10	Cyclic nucleotide-binding domain
003034	31	42	0.74	10−9	DNA-binding SAP

Table 4

Most frequently observed GO terms which were affected by alternative splicing variants

GO ID	GOs in alternative splicing locus	GOs in all locus	Ratio	Significance of enrichment (P-value)	GO term
0003676	451	1112	0.41	<10−16	Nucleic acid binding
0003700	327	518	0.63	<10−16	Transcription factor activity
0003677	276	603	0.46	<10−16	DNA-binding
0004713	164	318	0.52	<10−16	Protein tyrosine kinase activity
0005215	164	299	0.55	<10−16	Transporter activity
0008270	148	276	0.54	<10−16	Zinc ion binding
0005520	73	79	0.92	<10−16	Insulin-like growth factor-binding
0005524	379	967	0.39	10−14	ATP binding
0003824	190	429	0.44	10−13	Catalytic activity
0016491	116	237	0.49	10−11	Oxidoreductase activity

For the subcellular localization signal predictions, PSORT II () was run as indicated previously (21) and for the predictions of the transmembrane domains, TMHMM () was run with the cut-off value of 0.8. In cases where a protein motif(s) and/or a predicted localization signal(s), including transmembrane domains, were altered between two alternative splicing variants belonging to the same locus, the corresponding locus was defined as a ‘Motif/GO/Subcellular localization/Transmembrane domain-changed’ locus.

Manual procedures to inspect the identified alternative splicing variants

The results of the computational identification and annotation of the alternative splicing were visually inspected by the members of the alternative splicing annotation team of H-Invitational meetings by using the G-integra human genome browser [(12); Supplementary Figure 1] to verify the accuracy of the detected alternative splicing. Several items were manually evaluated, including correct discrimination of truncated cDNAs, proper identification of the alternative splicing types, and the corresponding functional annotations related to them. Further details of the annotation items and the records of the manual inspection are available at our web site (see below). The results of the manual and computational annotations were compared with each other, and in cases where the results were consistent between the two approaches, the alternative splicing and the related annotations were defined as ‘validated’. The results obtained for the annotations for each of the loci were made public and freely available at our web site ().

Identification of uncommon patterns of alternative splicing

Three ‘uncommon’ patterns of alternative variations were defined as follows: (i) ‘bridged’: a locus in which two alternative splicing variants are arrayed tandemly without sharing any exons and another transcript ‘bridged’ these two variants, sharing at least some of its exons with both of them; (ii) ‘nested’: a locus in which protein coding sequence (CDS) of one alternative splicing variant is not shared with another variant and (iii) ‘multiple CDS’: a locus in which different ORFs >200 amino acids in length are annotated independently for different alternative splicing variants having overlapping CDSs of different reading frames.

RT–PCR of the ‘bridged’ transcript

RT–PCR was performed using the primers: primer A 5′-CGTGAGCTCGCCCGCCAGAAG-3′; primer B 5′-TCCAACTCCAGCTCCACATC-3′; primer C 5′-CGAGATGACGGGCTTTCTGC-3′; primer D 5′-GGAATGCCATCGGTGCTGG-3′; primer E 5′-CCGACTATGCAGAGGAGAAG-3′; primer F 5′-GCGTTCTGCTGCTGCTCGAG-3′; primer (GAPDH fw) 5′-TCGGAGTCAACGGATTTGGT-3′; primer (GAPDH rv) 5′-TGACGGTGCCATGGAATTTG-3′, using ABI Prism 7900 Real Time PCR (ABI) with standard reaction conditions. The template RNAs (50 ng for each PCR) used were an RNA panel (BD Biosciences). For a negative control, 50 ng of human genomic DNA (Promega) was used as a template.

RESULTS

Identification of alternative splicing variants by manual and computational methods

A total of 56 419 human full-length cDNAs (Supplementary Table 1) were mapped onto the human genome and 55 036 cDNAs were unambiguously mapped. The mapped cDNAs were clustered into 24 425 loci, resulting in 2.3 cDNAs per locus on average. Single exon transcripts and the sole transcript in the locus were then removed. As a result, 10 127 loci contained at least two cDNAs (for further details on the mapping and clustering procedures, see Materials and Methods). These 10 127 loci, containing 35 030 cDNAs in total, were subjected to computational and manual inspection schemes to find if alternative splicing variants (defined as ‘complete forms of the transcripts’, with fully sequenced cDNAs) were included. We used this strategy since we considered that both the manual and computational methods should have advantages and disadvantages. Concerning the manual annotation, human errors are inevitable and for computational methods, detection of spurious alternative splicing due to various errors/ambiguities inherent to automated analyses is problematic. In order to maximize the accuracy of the analysis, the results of the manual and computational analyses should be compared with each other. When the results from computational analyses were ‘approved’ by manual inspections, the results were regarded as ‘validated’ (for the criteria of the annotation, see Figure 1).

Out of 35 030 cDNAs, the cDNAs annotated to be derived from ‘5′/3′ end-truncated’ or ‘immature’ transcripts either by manual or computational annotation were also excluded. In total, 5308 (15%), 706 (2%) and 787 (2%) cDNAs were defined as ‘5′ end’, ‘3′ end’ and ‘5′/3′-both side’ truncated, respectively (for further confirmation of the 5′ end completeness of each of the transcripts, see Supplementary Figure 2). A total of 1297 (4%) cDNAs were defined as ‘immature’. Also note that as for 3913 (11%) cDNA sequences, some of the annotators reported concerns that they might contain sequence problems due to cloning errors. The largest population of them was suspected spontaneous deletion of a part of cDNA insert in bacteria. The cDNAs supported by no ESTs were also identified. However, we did not remove these indecisive cDNAs, because it was not always straightforward to distinguish them from non-canonical rarely-occurring, but, biologically interesting alternative splicing (also see below). We provide independent statistics using the dataset which did not include them in Supplementary Table 2. Also, caveats for possible errors were precisely annotated for each of them in our web site ().

As a result, a representative ‘validated’ dataset of 6877 loci (68%), in which 18 297 ‘unique’ alternative splicing variants (2.7 variants per locus) were made of 37 670 alternative splicing exons (2.1 exons per variant; also see Supplementary Figure 3), was obtained and used for the subsequent analyses (Table 1). Also, with respect to each of the analyses described below, we chose the same strategy; to perform computational calculations first and then manually check the results. Therefore, each of the subsequent statistical analyses was based on the data which had passed both computational and manual inspections. A schematic representation of the computational calculations and the web interfaces used for the manual checks are shown in Supplementary Figure 1. The final results of each of the annotations were made public and freely available from our web site.

Interestingly, a surprisingly large population of the alternative splicing variants identified in the present study (as full-length forms) did not match in Ensembl (). If H-Invitational transcripts were defined as identical to Ensembl transcripts, when all of the exon–intron boundaries were corresponded to those of Ensembl transcripts with 10 bp allowance, we found 11 704 out of 18 297 (64%) of the H-Invitational transcripts were represented in Ensembl. When this was counted at the locus level, 6284 out of 6877 (91%) of the alternative splicing locus presented here contained novel alternative splicing relationship. This low level of overlap might be reflecting the fact that Ensembl is mainly based on the conservative analyses of the EST sequences and do not put much stress on the full-length cDNAs.

Patterns of the identified alternative splicing variants

Using the obtained dataset of the 18 297 alternative splicing variants from 6877 loci, we first examined the genome-wide features of alternative splicing in terms of their complete variants. Although similar analyses have been carried out by several groups (22–24), most of them were based either on partial EST sequences or analyses using smaller datasets or smaller numbers of cDNAs subjected to limited analyses (25).

The alternative splicing variants were first classified with regard to the splicing types, as often employed in previous studies (8). As shown in Figure 2, ‘cassette’ type alternative splicing variants were most frequently observed in our dataset, which is consistent with a previous result for the alternative splicing genes located on chromosome 22 (26). This type of alternative splicing variants may be preferred, because diversification of the transcripts can be achieved more flexibly than with the other categories of splice variations. ‘Retained intron’ type alternative splicing variants were observed in 1970 loci. However, it became a concern to us that cDNAs derived from unspliced, immature forms of transcripts might be classified as ‘retained intron’ types. To address this concern, we further checked how many of the alternative splicing variants of this category might be subject to nonsense-mediated decay [NMD; see Materials and Methods for the procedure to identify possible NMD-target variants (27)]. In 402 cases, one of the alternative splicing variants were annotated as a possible NMD-target, thus, might not have biological relevance. However, in the remaining cases, the alternative splicing variants with a ‘retained intron’ seemed to encode proteins with no explicit defects. So, it may indicate that most variants with ‘retained introns’ were beneficial for the diversification of human proteome.

Figure 2

Patterns of the identified alternative splicing.

We then examined the sequences of the alternative splicing exons and found that remnants of Alu-like elements were detected in a significant population (12%), which is considerably higher than the frequency in constitutive exons (2%; Table 5). This is consistent with previous results demonstrating that integration and subsequent alteration of the Alu-like elements play a significant role in the birth of alternative splicing exons (28). We also examined exon–intron junctions and detected canonical GT-AG (22) splice junctions in 156 181 (98.7%) out of 158 208 sites in total. Non-canonical AT-AC sites, the processing of which is using a splicing machinery with different and specific components (29), were also found in 89 cases. Another non-canonical splice site, GC-AG, were found in 57 cases. Taken together, we found that non-canonical junctions were enriched in alternative splicing exons (3.4%) compared to constitutive exons (0.8%). It is also intriguing those ESEs (20), which are suggested to play a role in efficient splicing, were also less frequent in alternative splicing exons. These findings suggest that the alternatively spliced junctions may be less deterministic than constitutive junctions, thus allowing versatile patterns of splicing events.

Table 5

Characteristics of the identified alternative splicing exons

	Exon-intron junction type		Containing Alu-like element	Containing ESE	Total
	Canonical	Non-canonical
Alternative splicing exon	26 888 (96.6%)	954 (3.4%)	12%	8%	27 842
Constitutive exon	129 293 (99.2%)	1073 (0.8%)	2%	10%	130 366
Total	156 181 (98.7%)	2027 (1.3%)	4%	9%	158 208

Positions of the identified alternative splicing exons

In our dataset of the 6877 alternative splicing genes, alternative splicing events located at the 5′ end, internal and 3′ end exons were observed in 4568, 5565 and 2933 genes, and consisted of 7494, 11 156 and 4940 ‘5′ end’, ‘internal’ and ‘3′ end’ alternative splicing variants, respectively (Table 1). Considering that the total numbers of the exons examined were 18 297, 139 911 and 18 297, the average frequency of the alternative splicing exons were 0.41, 0.08 and 0.27 for each of the positions, respectively. It was intriguing that the ‘5′ end’ alternative splicing variants formed the most frequent category. We further examined these 5′ end alternative splicing exons and found that, among the 7494 ‘5′ end alternative splicing genes’, 3495 genes (47%) contained 5′ end alternative splicing variants, which were separated by more than 500 bp from each other. It is likely that these alternative splicing exons were produced as a consequence of the use of alternative promoters (30). Therefore, the biological significance of those 5′ end alternative splicing variants should be comprehensively analyzed together with the results of recent studies demonstrating that alternative use of promoters is prevalent in human genes (18,31,32). It is possible that the 5′ end alternative splicing exons were most abundant because of the requirements for the diversification of the transcriptional modulation and vice versa.

In order to characterize the biological consequences of the identified alternative splicing, we examined the relative position of the alternative splicing variants compared to their CDSs in 6555 genes in which two or more alternative splicing variants are annotated as protein coding. We found that alternative splicing variants were located in the 5′-untranslated regions (5′-UTRs), CDSs, and 3′-UTRs in 4750, 13 409 and 1034 alternative splicing variants in 3216, 6005 and 797 genes, respectively. In the majority of the alternative splicing genes, alternations in their CDSs were observed. Notice that 3′ end alternative splicing was counter-selected in our analyses, by performing in silico filtration of NMD. Although in 80% of the alternative splicing variants, the alteration of the polypeptide length was less than 200 amino acid, 3% of the alternative splicing variants resulted in alteration of the polypeptide length of more than 500 amino acid (Figure 3). For this population, alternative splicing should have the most significant impact on the functions of the encoded proteins. It is even possible that these loci encode two functionally different proteins simultaneously, in which the fraction of the transcripts that is shared and their biological functional similarities may widely vary. Extreme examples of these cases will be discussed later.

Figure 3

Distribution of the length difference between the alternative splicing variants. The percentages show the populations belonging to the corresponding groups.

Possible biological relevance of the identified alternative splicing variants to the diversification of protein functions

The influence of the identified alternative splicing variants on the encoded protein sequences and their possible biological functions was further evaluated from the following viewpoints. The complete records of each of the analyses are available for each of the entries from our web site ().

Motif

Protein motifs were frequently affected by the alternative splicing located within the CDSs. Out of 6005 alternative splicing variants residing in the CDSs, 3015 were located in protein motifs (Table 6). As shown in Figure 4, the alternative splicing variants influenced a wide variety of protein motifs. For example, we identified a novel alternative splicing variant, AK093798 in the IκB kinase epsilon gene (IKKɛ; NM_014002). This novel variant specifically lacks the kinase domain, while the rest of the CDS remained intact (Figure 4A). The IKK complex plays a pivotal role in immune and inflammatory responses by transmitting a variety of signals to a transcription factor, NF-κB (33). It has been reported that a kinase-deficient dominant-negative mutant of IKKɛ blocks the induction of NF-κB invoked by T cell receptor but has no effect on its activation invoked by TNFα (34). The natural kinase-deficient variant of IKKɛ represented by AK093798 could serve as a modulator between these two signaling pathways, and thus provide cells with the opportunity to regulate the relative amounts of the signals that the two pathways receive at a single point.

Table 6

Numbers of the genes in which alternative splicing variants should influence the possible protein functions

	#Locus	#cDNA
Alternative splicing affecting functional annotation (total)	4481	12 542
Motif-changed	3015	8727
Subcellular localization-changed	2982	8624
GO-changed	1779	5179
Transmembrane domain-changed	1348	3933
Uncommon alternative splicing pattern (total)	316	1033
Bridged	129	604
Nested	172	390
Multiple CDS	27	56

Figure 4

Examples of the alternative splicing variants detected as ‘motif-changed’ (A), ‘subcellular localization-changed’ (B) and ‘transmembrane domain-changed’ (C). Exons and introns are represented by green boxes and lines, respectively. The violet boxes are protein coding regions and yellow boxes are alternative splicing exons. The positions of the detected motifs and transmembrane domains are shown beneath the transcripts. In the uppermost panel, GO terms attached to the transcript indicated by the lower line are shown.

In total, among 6877 genes with alternative splicing, alterations of the motifs due to the alternative splicing were observed in 3015 genes (44%; Table 2). In 2508 genes (36%), alternative splicing did not change the annotated motifs. The remaining 1354 genes (20%) contained no annotated motif. The alternative splicing genes had a higher frequency of motifs overall (80%), as among 17 548 non-alternative splicing genes, only 7241 genes (41%) contained annotated motifs (P < 10−16; Table 2). We also found that alternative splicing exons were enriched in motifs. On average, one motif was contained per 1.6 alternative splicing exons, while the average frequency of motifs was only one per 3.0 constitutive exons. A similar tendency was also suggested from a recent EST-based study (35). The motifs should be actively associated with alternative splicing, and in many cases direct switching of the motifs is enabled by alternative splicing.

Subcellular localization

We also examined whether the predicted protein motifs determining the subcellular localization signals of the proteins, such as secretion signal peptides, mitochondria targeting signals and transmembrane domains, were affected by the alternative splicing (Table 6). For the subcellular localization signal predictions, PSORT II was used. The subcellular localization signal was predicted for each of the alternative splicing variants, and in cases where the alternative splicing variants were predicted to localize in different subcellular compartments, the alternative splicing were categorized as ‘subcellular localization-changed’ (shown in Figure 4B). Similarly, ‘transmembrane domain-changed’ alternative splicing were identified using TMHMM. In total, 2982 subcellular localizations and 1348 transmembrane domains were affected (Table 6). Figure 4C shows a case in which the transmembrane domain was altered. The most frequently observed switching of the predicted subcellular localization signal was between ‘nuclear’ and ‘cytoplasm’ (2455 cases). Switching between ‘secretory’ and ‘plasma membrane’ was detected in 1145 cases. These findings suggested that the proteins produced from the same loci as a result of the alternative splicing are actively utilized in a multi-faceted manner at different locations or compartments in the cells. A similar tendency of the use of several localization signals in the protein isoforms was also indicated by a recent bioinformatics study (36).

GO

In a large population of the genes (1779 genes; 27%; Table 6), the GO terms attached to the transcripts were altered between the alternative splicing variants even within the same locus because the protein motifs and subcellular localizations were affected by the alternative splicing variants, as mentioned above. For example, in the case shown in Figure 4A, the GO term, ‘kinase activity’ (and related terms), was assigned to the kinase domain-containing variants, but not to the kinase negative-variants. Some of the motifs and GO terms were significantly enriched in alternative splicing genes (Table 4). It is especially noteworthy that motifs and GO terms associated with signal transduction and transcriptional regulation were most frequently influenced by alternative splicing (Tables 3 and 4). Fine adjustments of the resultant protein functions might be one of the most important functions of alternative splicing in higher organisms, such as humans.

Uncommon patterns of alternative splicing

During manual annotation, we noticed that a number of transcripts undergo novel patterns of alternative splicing. Definitions of them are described in the legend forFigure 5.

Figure 5

Examples of the ‘uncommon’ patterns of alternative splicing; ‘bridged’ (A), ‘nested’ (B) and ‘multiple CDS’ (C). These ‘uncommon’ patterns of alternative variations were defined as following: i) ‘bridged’: a locus in which two alternative splicing variants were arrayed tandemly without sharing any exons and another transcript ‘bridged’ these two variants, sharing at least some of its exons with both of them; ii) ‘nested’: a locus in which CDS region of one alternative splicing variant was not shared with another variant and iii) ‘multiple CDS’: a locus in which different ORFs >200 amino acid in length were annotated independently for different alternative splicing variants having overlapping CDSs of different reading frames. In the lower panel of (A), the results of RT–PCR are shown. Each photograph shows the amplicons of the indicated RT–PCR using the indicated primers. Tissue origins of the template RNAs are shown in the margin. The asterisk indicates a non-specific band. The coloring of the figures is the same as in Figure 4.

Bridged

As shown in the upper panel of Figure 5A, a cDNA, AK000438, was identified as an alternative splicing variant of AK000479 and AF161485. However, it was more likely that AK000438 was overlaid on two adjacent loci rather than ‘a variant of either of the loci’, representing a ‘bridging’ transcript from two genes, the SERF2 gene (NM_005770) and the HYPK gene (NM_016400). It is noteworthy that both genes are associated with neuromuscular diseases. These genes were originally identified respectively as a candidate modifying gene for spinal muscular atrophy (37) and a Huntingtin interacting protein (38). We excluded the possibility that AK000438 was derived from artifacts, such as chimeric transcripts produced during the cDNA cloning process because of the facts that: (i) genes immediately adjacent to each other are bridged; (ii) the transcript was connected exactly at exon–intron junctions satisfying the GT-AG rule; (iii) there are supporting dbEST sequences which correspond to the ‘bridging’ transcript and (iv) RT–PCR using primers sets at the SERF2 and HYPK transcripts produced direct evidence that the ‘bridging’ transcript exits in vivo (Figure 5A lower panel). A previous study demonstrated that this kind of ‘bridging’ transcript is produced between the GALT gene and the IL-11Ra gene, and results in an mRNA encoding a protein retaining the functions of both proteins (39). This kind of ‘bridging’ transcript has been reported to be especially frequent in Caenorhabditis elegans, (179 out of 3829 alternative splicing genes) as are complex genes generating two completely distinct proteins in addition to a fusion of both. Recent reports suggested that such transcripts seem abundant in humans as well [(40); also see ]. Although, further experimental characterization is required to clarify their biological consequences, it is possible that such bridged transcripts play some roles by combining the respective gene functions in general.

Nested

Figure 5B shows another uncommon type of alternative splicing (‘nested’). In this case, all exons of AB508739 except for a shared 3′ end exon, are embedded in the last intron of AK130874 (Figure 5B, upper panel). Neither a mapping error of the cDNAs nor an assembly error of the genome sequence is likely to account for this observation, since the sequence identity/coverage between the genome and cDNA sequences are almost 100% for each of the exons and genomic sequences are ‘finished’ in this region. In this case, these transcripts may share some modulatory elements, such as binding sites of regulatory proteins or non-coding RNAs, but their biological functions should be completely different. A similar specific case of alternative splicing was also described with regard to the 5′ end exons, on which promoters seemed to be shared. (Figure 5B, lower panel). In a very recent paper, two otherwise independent genes were shown to be co-expressed from one promoter, giving yield to tissue-specific expression of an alternative form of an otherwise B-cell specific gene (41).

Multiple CDS

Alternative splicing even seemed to allow a particular genomic sequence to encode two distinct amino acids in different reading frames. The last exon of AK097244 overlapped with the second exon of AK000272 (Figure 5C, upper panel). Both of the exons were coding exons, however, their reading frames were different. Similarly, the reading frames of the second exon of AK096258 and fifth exons of BC029781 were different (Figure 5C, lower panel). Interestingly, in this gene, another alternative splicing variant, BC043484 used the latter-type reading frame for the N-terminal half of the protein and the former-type reading frame for the C-terminal half. In this variant, sixth exon seemed to serve as a ‘frame switching’ exon. Although the gene function of this locus still remains unknown, it is possible that the novel proteins discovered in the present study should also play roles which are distinct from but, at the same time, correlate with each other. Consistent with this possibility, emerging evidence also suggested that the cases in which mRNAs transcribed from a single locus are used as templates for two independent proteins seem not extremely rare (42). Further experimental validation should clarify which gets translated to produce functional proteins and which have regulatory roles.

The numbers of cases identified for each of the above-mentioned ‘uncommon’ patterns of the alternative splicing are shown in Table 6. In any case, it would be more appropriate to regard these cases as two ‘genes’ merged into a single ‘locus’ rather than as mutually ‘alternative’ variants occurring from the same locus, since, obviously, it is not probable that these variants share the same protein function. If we consider the evolutionary origins of these loci, they could originally have been two neighboring genes that then evolved to share some of their exons as a result of mutational changes. In each case, the extents to which multiple loci are merged vary. In this study, we applied strict criteria to select them (see Materials and Methods). Therefore, the numbers shown in Table 6 should be the minimum values, suggesting that it is not extremely rare to observe such ‘uncommon’ alternative splicing events in cells. Such mechanisms may allow for further diversification of the transcriptome of human genes. Although careful evaluations should be necessary (i.e. whether the cDNAs are originated from cancerous cells) before concluding whether the particular examples described here should have biological relevance, if any, such mechanism would enable the multi-faceted use of the biological information encoded in a given genomic region.

It should also be noted that without a combination of manual and computational analyses, this kind of alternative splicing variants would have been overlooked. By performing detailed manual and computational analyses, we could precisely identify these novel patterns of alternative splicing. Also, if no more than partial information of the alternative splicing exons had been available, these findings would have been utterly impossible. Judging from these results, alternative splicing is likely to be used in a more versatile manner to enable the diversification of the gene functions of human genes than was previously thought.

DISCUSSION

In this paper we described the genome-wide identification and characterization of alternative splicing variants of human gene transcripts based on completely sequenced human full-length cDNAs in details. Starting with 56 419 cDNAs, we identified 18 297 complete alternative splicing variants at 6877 loci of human genes and precisely annotated each of them.

Lander et al. (6) proposed that there should be five alternative splicing variants per locus in human genes. Although, our dataset used here is the largest one from a full-length cDNA collection, it is much smaller than that from dbEST, which includes more than 7 million ESTs. Therefore, some of the alternative splicing variants that actually exist in human genes might not be represented in our dataset (we identified alternative splicing variants from 68% of the loci examined, with 2.7 variants per locus). However, our dataset should have two major advantages that largely compensate for the shortage of coverage. First, our dataset had been validated by heuristic annotations, and therefore various computational errors which could be easily discriminated by the human eye were excluded. Secondly, each of the identified alternative splicing variants was supported by full-length cDNAs whose sequences had been completely determined. This enabled us to assure that each of the alternative splicing variants identified in this study corresponds to a complete form of a particular transcription unit. This feature was especially important when the relevance of the alternative splicing variants to the protein motifs or various subcellular localizations was evaluated. If comprehensive information of the cDNA sequences had not been available, those analyses would have been less reliable. Sometimes protein motifs are embedded over a wide region of the protein sequences, and all of the combinations of the alternative splicing exons may not be allowed. Besides, for certain types of subcellular targeting signals, such as signal peptides, the position within the protein sequence is critical.

Recently, it was reported that many of the alternative splicing variants seemed not to be evolutionarily conserved, and thus the biological significance of their existence might be questionable (43,44). The same argument may be used in the other direction as well. Alternative splicing could be the easiest road to diversification of the species. Future analyses of our full-length alternative splicing data should be useful to clarify which alternative splicing should be evolutionarily wobbly and which should provide a molecular basis for various species-specific biological features. Indeed, it will be important in future studies to discriminate which of the identified transcript variants have a general raison d'etre and which play species-specific roles, because lack of such knowledge would severely restrict the potential power of the comparative genomic approaches utilizing the genomic sequences of many organisms which will be determined in the next few years. The information described here together with the availability of the accompanying physical full-length cDNA clone resources should lay firm groundwork for exploring how alternative splicing generates the functional diversification of the human transcriptome and proteome.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR online.