Uso de técnicas computacionais no estudo da transcrição e regulação gênica em Homo sapiens e Mus musculus

Galante, P.A.F., Use of computational methods to study the transcription and gene regulation in Homo sapiens and Mus musculus. 2007. (134 pag) PhD Thesis – Graduate Program in Biochemistry. Instituto de Química, Universidade de São Paulo, São Paulo. Genes, nucleotide sequences necessary for the synthesis of functional molecules, are transcribed and regulated by extremely complex cellular and molecular processes. To understand when and in which tissues the genes are expressed, their functional isoforms, control regions and the factors involved in gene regulation is one of major challenges of modern molecular biology. Today, the availability of complete genome sequences and transcriptomes, together with the development of new computational methods allows the study of phenomena related to the transcription and gene regulation in a large scale. For example, it is possible to quantify, concomitantly, gene expression of thousands of genes in different tissues and analyze different aspects of their regulation. In this work we developed and applied computational methods to the study of four key aspects of gene transcription and regulation. In the first study, we addressed tissue specific gene expression through the study of genes that are preferentially expressed in the brain and ten different mouse brain regions. In the second study, we identified sequences that are potentially involved in the control of gene transcription through the study of motifs that are over represented in the promoter region of olfactory receptor genes. In the third study, we browsed the human for the presence of intron retention, a type of alternative splicing. In the fourth study, we addressed the transcriptoma complexity and gene expression regulation through the study of pair of sense-antisense genes in human and mouse. In all studies, our results allowed us to make specific conclusions about each phenomenon analyzed which showed us the importance of a large scale approach. In addition, we verified that our computational methods can be efficiently applied to the study of transcription and gene regulation in Homo sapiens and Mus musculus.


Dados de hibridização in-situ
Os dados de hibridização in-situ foram obtidos do Allen Brain Institute. A expressão de cada gene foi analisada manualmente com base nas imagens de cortes coronais e sagitais do cérebro completo de camundongos C57BL/6J adultos.
Um experimento de SAGE e de MPSS gera uma biblioteca que contém um conjunto de tags e a freqüência de cada uma. Então, o primeiro passo para se analisar estes dados é identificar (mapear) o gene que cada tag representa.

Gene Name
Gria2
Transcrição bidirecional múltipla corresponde aos loci contendo dois ou mais pares S-AS. Os conjuntos for subdivididos de acordo com a presença de mRNAs ou de ESTs no clusters.

Estudo dos eventos de retenção de introns
Primeiro, os nossos resultados mostraram que a retenção de intron é o tipo menos freqüente de splicing alternativo, ocorrendo em ~15% dos genes humanos. Em seguida, mostramos que os introns retidos apresentam características similares aos exons e não aos introns não retidos. Por fim, mostramos que a seleção está atuando no introns retidos, deixando-os similares às seqüências codificadoras. Agrupando e analisando todos os resultados, concluímos que estes introns retidos não são simplesmente produtos de transcrição espúria ou artefato experimental, mas sim eventos de splicing alternativo com uma função biológica no contexto celular.  Mouse odorant receptors (ORs) are encoded by >1000 genes dispersed throughout the genome. Each olfactory neuron expresses one single OR gene, while the rest of the genes remain silent. The mechanisms underlying OR gene expression are poorly understood. Here, we investigated if OR genes share common cis-regulatory sequences in their promoter regions. We carried out a comprehensive analysis in which the upstream regions of a large number of OR genes were compared. First, using RLM-RACE, we generated cDNAs containing the complete 5Ј-untranslated regions (5Ј-UTRs) for a total number of 198 mouse OR genes. Then, we aligned these cDNA sequences to the mouse genome so that the 5Ј structure and transcription start sites (TSSs) of the OR genes could be precisely determined. Sequences upstream of the TSSs were retrieved and browsed for common elements. We found DNA sequence motifs that are overrepresented in the promoter regions of the OR genes. Most motifs resemble O/E-like sites and are preferentially localized within 200 bp upstream of the TSSs. Finally, we show that these motifs specifically interact with proteins extracted from nuclei prepared from the olfactory epithelium, but not from brain or liver. Our results show that the OR genes share common promoter elements. The present strategy should provide information on the role played by cis-regulatory sequences in OR gene regulation.

Estudo da transcrição gênica através de técnicas computacionais
[Supplemental material is available online at www.genome.org. The sequence data from this study have been submitted to GenBank under accession nos. DR065530-DR065963.] Mammalian olfactory sensory neurons select, from >1000 possible choices, one single olfactory receptor (OR) allele to express (Ressler et al. 1993;Vassar et al. 1993;Malnic et al. 1999; but see . The receptor type that is chosen will not only determine the range of odorants to which this neuron will respond, but also axonal targeting to specific glomeruli in the olfactory bulb (Mombaerts et al. 1996;. OR gene choice is therefore fundamental for the functional organization of the olfactory system. How this choice is accomplished is, however, still unclear. Little is known about the role of cis-regulatory sequences in the regulation of OR gene expression. In studies using transgenic mice, different sizes of genomic DNA segments containing OR genes were tested for their ability to drive an OR expression similar to that of the endogenous gene. It was demonstrated that short pieces of DNA located upstream of the coding region, ranging from 460 bp to 6.7 kb, are sufficient for expression of the ORs M4, M71, and MOR23 (Qasba and Reed 1998;. However, large segments of ∼200 kb are required to obtain expression of MOR28 (Serizawa et al. 2000). Sequence comparison of the mouse and human genome revealed a 2-kb conserved sequence located ∼75 kb upstream of the MOR28 gene cluster. This region, denominated the H region, was proposed to work as a cis-acting locus control region (LCR) that would activate the expression of one single OR gene member from within the MOR28 gene cluster . Altogether, these re-sults indicate that cis-regulatory sequences may play important roles in OR gene choice.
To date, different combinations of transcription factor binding sites (TFBSs) have been identified in promoters of OR genes (Hoppe et al. 2000(Hoppe et al. , 2003Lane et al. 2001;; however, there is no evidence yet that these sites are directly involved in OR gene choice. A strong consensus sequence, the Olf-1 site (O/E-like site), was identified in the promoters of several olfactory specific genes, such as GnaI (formerly known as Golf), adenylyl cyclase III (AcIII), olfactory cyclic nucleotide gated channel (Cnga2), and olfactory marker protein (Omp) (Kudrycki et al. 1993;, and was also found in the promoter regions from some OR genes (Glusman et al. 2000b;Hoppe et al. 2003) but not from other OR genes (Hoppe et al. 2000;Lane et al. 2001).
There are >1000 OR genes dispersed throughout the genome (Young et al. 2002;. A genomic approach to identify potential regulatory cisacting sequences is to search for DNA sequence elements that are conserved in a large number of OR gene promoters. Promoter sequences, which are usually located proximal to and upstream of the transcription start site (TSS), can be retrieved from the available mouse genome sequence. This can be done by aligning full-length OR mRNA sequences with their counterpart genomic sequences. It is important that full-length mRNA sequences are used, because the transcriptional start sites can be located far away from the translational start sites. It has been demonstrated, for example, that ORs in the mouse P2 OR cluster have 5Јuntranslated regions (5Ј-UTRs) that range from 1.7 to 9 kb (Lane et al. 2001). Besides, it is known that many OR genes have 5Ј non-coding exons, and some undergo 5Ј alternative splicing Lane et al. 2001).
To date, full-length 5Ј cDNA sequences are available only for ∼30 different mouse OR genes, distributed over five different chromosomes (Bulger et al. 2000;Hoppe et al. 2000Hoppe et al. , 2003Lane et al. 2001;. Another study produced cDNA sequences representing >400 OR genes, many of them containing 5Ј-UTR sequences. However, they may not be full-length sequences, since they originated from a regular cDNA library, and not from a full-length cDNA library ).
Here we have used RLM-RACE with degenerate primers to produce cDNAs containing the complete 5Ј-UTRs for a total number of 198 mouse OR genes. Then, we aligned these cDNA sequences to the mouse genome so that sequences corresponding to one same OR gene are organized in one same cluster, making a total of 198 clusters. For each OR gene cluster, we determined the 5Ј-gene structure (exon and intron distribution) and TSS. The 198 promoter sequences were retrieved and browsed for common elements. Motifs that are common to a large percentage of the OR genes were found. Most motifs resemble O/E-and homeodomain-like sites. The O/E-like sites are localized within 200 bp upstream of the TSS and specifically interact with proteins extracted from nuclei prepared from the olfactory epithelium, but not with proteins extracted from brain or liver.

Amplification of OR cDNA 5Ј-ends
To obtain 5Ј-end sequences from OR cDNAs, we performed RNA ligase-mediated rapid amplification of 5Ј cDNA ends (RLM-RACE) using total RNA purified from mouse olfactory epithelium. This method has the advantage that only full-length transcripts (authentic 5Ј-end capped mRNAs) are amplified (Fig. 1A). In order to obtain 5Ј cDNA ends for a large number of OR genes, we used degenerate primers matching to conserved regions in ORs (Fig. 1B) because these primers can amplify the majority of the members of the OR family Malnic et al. 1999). The initial PCR reaction was performed using the 5Ј GeneRacer primer together with the reverse P27 degenerate primer. In order to eliminate PCR artifacts, a secondary (nested) PCR was performed using the first PCR product as a template and the 5Ј GeneRacer nested primer together with the P26R degenerate primer (Fig. 1A). The primary PCR product contained a heterogeneous mixture of cDNA fragments ranging from 0.8 to 2 kb in length (Fig. 1C, lane 2). This result was expected, since different OR cDNAs must have different 5Ј-UTR sizes. The secondary PCR reaction produced a similar range of cDNA fragments, except that their sizes were ∼350 bp smaller (Fig. 1C, lane  3). This was also expected, since the P26R primer matches to a region in the OR coding sequence that is located ∼350 bp upstream to the region matched by P27 (Fig. 1B). Since the region between the first AUG codon and TM-III in OR genes is ∼380 bp long, only the nested PCR products Ն380 bp long were gelpurified, cloned, and sequenced.
The same procedure was also performed using primers 5Ј GeneRacer and P8 degenerate primer for the primary PCR reaction (Fig. 1C, lane 5), and primers 5Ј GeneRacer nested and P27 for the secondary PCR reaction (Fig. 1C, lane 6). In this case, since the region between the first AUG codon and TM-VI in ORs is ∼750 bp long, only the nested PCR products Ն750 bp were analyzed.

OR gene clusters
We sequenced 1012 clones from their 5Ј-ends, and 80% of them correspond to ORs, indicating that our strategy preferentially amplifies OR sequences. In addition, 96% of the OR cDNAs contain 5Ј sequences upstream of the predicted initial AUG codon, indicating that full-length mRNAs were amplified. Sequence analysis showed that 5Ј-RACE products were obtained for a total number of 198 different OR genes, corresponding to ∼17% of the complete mouse OR gene repertoire (Table 1). Only nine of the OR genes are pseudogenes.
OR sequences are classified into two phylogenetic classes, referred to as Class I (fish-like) and Class II (terrestrial-specific) ORs Freitag et al. 1995;Glusman et al. 2000a). The Class I ORs constitute ∼12% of the mouse OR repertoire . Nine of the 198 ORs (4.5%) are Class I ORs, indicating that the method we used amplifies members of the two OR Classes, although it may favor amplification of Class II ORs.
The OR genes for which cDNA sequences were obtained are distributed among all of the mouse chromosomes previously shown to contain OR genes , except for chromosome 3 ( Table 1). The majority of the sequences corre- Figure 1. The strategy used for the generation of OR complete 5Ј-end cDNA sequences. (A) Truncated mRNAs were dephosphorylated using calf intestine phosphatase (CIP) so that they cannot participate in subsequent ligation reactions. The RNA was then treated with tobacco acid pyrophosphatase (TAP) to remove the 5Ј-cap structure (represented by black balls) from intact full-length mRNA, and the GeneRacer RNA oligonucleotide was ligated to the decapped mRNA. Reverse transcription was performed using oligo(dT) primers. To obtain 5Ј-ends, PCR was done using the GeneRacer 5Ј-primer and a degenerate primer directed toward conserved OR regions (P8 or P27). Nested PCR was then done using the GeneRacer 5Ј nested primer and another OR degenerate reverse primer (P27 or P26R). spond to genes in chromosomes 2 and 7, which had been previously shown to contain the higher numbers of OR genes (Young et al. 2002;. In addition, the 198 OR genes can be subdivided into 102 (41%) out of the total 248 mouse OR subfamilies (where all members of a subfamily are Ն60% identical to all other members in amino acid sequence, as described by ; see Supplemental Table 1). We therefore believe that our sequences are representative of a random sample of the mouse OR genes. However, it is important to note that the present method may favor the amplification of OR cDNAs that have short 3Ј-UTRs or that are highly expressed in the olfactory epithelium.
We next used BLAST (Altschul et al. 1990) and Sim4  to align all of the cDNA sequences with the mouse genome sequence. The previously annotated mouse OR genes (Young et al. 2002; were also included in the alignment, to help with the localization of the cDNAs 5Ј-UTR regions. Each cDNA sequence aligned to one single genomic region, and the cDNAs that aligned to the same genomic region as one of the annotated ORs were considered to correspond to that particular OR gene. In this way, we obtained a total of 198 OR clusters, where each cluster corresponds to one different OR gene. Each one of the clusters contains at least one cDNA sequence, the largest cluster contains 71 sequences, and 54% of the clusters contain more than two sequences (Supplemental Table 1 Tables 1 and 2), and 39 showed alternative splicing (in this case, only clusters containing more than two cDNA sequences were analyzed).

5Ј-structure of the OR genes
On average, the 5Ј-UTR is 189 bp long, ranging from 32 to 659 bp. The 5Ј-UTR exons range from 9 to 403 bp (Supplemental Table 2). Most of the OR genes (72%) have only one 5Ј-UTR exon, 23% have two, and 3% have three (Supplemental Tables 1 and 2). The 5Ј-UTR introns show a wide variation in size, ranging from 91 bp to 22.5 kb. A large fraction of the OR genes (46%) have introns with sizes between 2 kb and 4 kb, and 34% of the genes have introns >4 kb.
It is predicted that in most eukaryotic mRNAs translation initiates at the first AUG starting from the 5Ј-cap (Kozak 1999). Therefore, one would expect not to find AUG codons or upstream ORFs (uORFs) in the 5Ј-UTRs. However, we found that 108 of the 198 OR genes have uORFs at least 10 codons long. It has been suggested that uORFs could be involved in down-regulation of protein translational efficiency . Further analysis should clarify whether these uORFs play a role in OR gene regulation.
We also found that a significant percentage of the OR genes (18%) have an in-frame upstream AUG, indicating that these genes code for OR proteins with a longer N-terminal region than the one originally predicted from their genomic sequences.

OR promoter regions
The generation of 5Ј full-length cDNAs allowed us to precisely determine the TSSs for the 198 OR genes. On average, the TSSs are located 4.3 kb upstream of the initial AUG codon, the furthest TSSs being located 22.5 kb away and the closest only 18 bp away. Sequences (600 bp) upstream of each TSS were excised from the mouse genomic sequence and analyzed. We first screened the  118 (37) 20 (1)  10 47 (13)  9  11 35 (20) 21 (2)  12 --13 12 (2)  4  14 32 (4)  6  15 6 (3) 3 (1)  16 29 (8)  18  17 36 (17)  5  18 - The number of annotated mouse OR genes from  and of mouse OR genes for which cDNA sequences were obtained in the present study. Numbers of pseudogenes are indicated in parentheses. and found that only 7% of the total sequences contain repeats or low-complexity regions. Typical TATA-boxes were found in only a small number of the OR gene promoters (35%), consistent with previous reports (Hoppe et al. 2000;Lane et al. 2001).
Because most TFBS are usually short, they can occur very frequently in the sequences, making it difficult to identify significant sites. In order to reduce the false-positive predictions, we decided to search for motifs that are common to a large fraction of the promoter sequences and thereby identify elements that are more likely to be functionally important. To do this, we used the Gibbs recursive sampler , Consensus , and Weeder (Pavesi et al. 2004) tools, which were designed to locate common elements in collections of unaligned DNA sequences. We found several motifs that are shared by the OR gene promoter sequences (Fig. 2). A closer inspection of the motifs revealed that although they are diverse, the majority of them resemble Olf-1 (O/E) like sites (Figs. 2 and 3;. The O/E-like motifs can be divided into four groups, denominated M1-M4. Motifs in each one of these four groups show different conserved nucleotide sequences (Fig. 3). We also found motifs that resemble the homeodomain sites, previously shown to be located in proximity to O/E-like binding sites in OR promoter genes (Fig. 2;. The spatial distribution of O/E-like motifs M1-M4 in the promoter sequences is shown in Figure 4A. All of the motifs are concentrated near the TSSs (between +1 and ‫002מ‬ bp). Differently, the homeodomain-like sites show a broader distribution over the entire extent of the 600-bp sequence (Fig. 4B).
The O/E-like motifs and the homeodomain-like motifs were found, respectively in 87% and 95% of the OR promoter sequences (Table 2).

DNA motifs form complexes with olfactory nuclear proteins
To evaluate the biological significance of the motifs, we performed gel shift assays. As shown in Figure 5, motifs M1-M4 formed DNA-protein complexes in the presence of the olfactory nuclear extracts, but not in the presence of brain and liver nuclear extracts. The different complexes show similar electrophoretic mobilities, but the intensities of the shifted bands vary. Formation of these complexes was inhibited by pre-incubation of the binding reaction with a 100-fold excess of the corresponding specific unlabeled oligonucleotides (Fig. 5A). We next evaluated the specificity of the DNA-protein complexes. Mutated motifs M1-M4, where the conserved nucleotides were changed into different ones, were unable to form stable DNA-protein complexes with olfactory nuclei proteins (Supplemental material 3). Altogether, these results show that proteins present in nuclei of olfactory epithelium cells, but not in liver and brain extracts, specifically bind to motifs M1-M4.

Discussion
The role of cis-regulatory sequences in the regulation of OR gene expression is not understood. In the present study, we carried out a comprehensive analysis in which the promoter regions of 198 OR genes were compared.

Common OR gene promoter motifs
By comparing a variable set of OR genes, we aimed to identify common promoter elements that may be involved in the general mechanisms of OR gene regulation. Strikingly, 87% of the OR gene promoter regions contain O/E-like sites, and 95% contain homeodomain-like sites. No other types of motifs were found to be overrepresented in these sequences, although we cannot ex- Figure 3. Nucleotide sequence alignment of the conserved sequences in O/E-like motifs M1-M4 and comparison with Olf-1 binding sites. Ten motif sequences were randomly selected for each one of the motifs (M1-M4) and manually aligned to the Olf-1 sites (from Omp, AcIII, Cnga2, and GnaI) . Identical nucleotide positions in each group of sequences are indicated by asterisks. The names of the OR genes from which the motifs were retrieved are indicated (e.g., M1_720: motif 1 from olfr720). The distribution of TAATTG homeodomain sites, previously shown to be present in many OR promoter regions  and to be involved in the regulation of M71 expression , is also shown.
clude the possibility that the algorithms we used may have precluded us from finding additional sites.
The roles of the O/E-like and homeodomain-like proteins in OR gene expression are still unclear. It has been demonstrated that expression of the M71 OR gene can be driven by a short region upstream of the TSS (a minimal promoter region) containing an O/E-like site and a homeodomain-like site . Mutational studies using transgenic mice indicate that both sites are required for normal OR gene expression . Yet, disruption of olf-1-like genes does not alter OR gene expression Wang et al. 2003), possibly because of the functional redundancy of the multiple O/E family members expressed in the olfactory epithelium (O/E1, O/E2, O/E3, and O/E4) . Nevertheless, it was demonstrated that O/E2 and O/E3 mutant mice show defects in the projection of olfactory neurons to the olfactory bulb, indicating that the O/E genes' function may not be completely redundant .
The LIM-homeodomain protein Lhx2 was shown to bind to the homeodomain site in the M71 OR promoter region (Hirota and Mombaerts 2004). Lhx2 knockout mice do not express ORs, but since they also lack mature olfactory neurons, there is no evidence so far that this homeodomain protein is directly involved in OR gene expression (Hirota and Mombaerts 2004;Kolterud et al. 2004). It is also possible that other homeodomain proteins than Lhx2 bind to the homeodomain sites in the OR gene promoters.
Although the vast majority of OR gene promoters have O/Elike sites, the structure of these sites is variable (see Figs. 2 and 3). It has been previously shown that the different O/E proteins possess similar DNA-binding properties . Our results indicate that the M1-M4 motifs interact differently with proteins from olfactory epithelium nuclei (Fig. 5). Several possibilities could explain the different DNA-protein complex affinities. The motifs could bind to different O/E proteins, or to alternatively spliced versions of these proteins . It is also possible that the same O/E protein types could bind to the M1-M4 motifs, but with different affinities.
Our findings suggest that different OR gene promoters are bound by different combinations or amounts of O/E-like proteins. The consequences of these differential interactions for OR gene regulation are unknown. It is known that different ORs are expressed in different levels . One interesting possibility is that the types of O/E-like sites in an OR gene promoter region may determine its probability of being transcribed. The identification of the proteins that interact with each one of the motifs and the analysis of the expression patterns of OR genes that have different motifs should clarify the role of these O/E-like sites in OR gene regulation.

Promoter DNA elements and OR gene regulation
Different models for OR gene regulation have been considered to date (for review, see Serizawa et al. 2005;Shykind 2005). It has been recently demonstrated that the monoallelic expression of an OR gene is regulated by a negative feedback mechanism that requires a functional OR protein . In addition, it was shown that immature olfactory neurons that express a given odorant receptor can switch receptor expression at a low frequency, while neurons expressing a mutant (nonfunctional) OR can switch expression with a greater probability . Based on these results, a new model has been proposed (Serizawa et al. 2004;Shykind 2005). In this model, after an OR gene is stochastically selected for expression by a limiting factor, its corresponding OR protein product mediates a feedback signal that results in the maintenance of the receptor choice.
Here we show that a collection of random OR genes will have the same types of cis-regulatory elements, suggesting that these common promoter elements are likely to play an important role in OR gene expression. It is possible that enhancers or LCRs interact with elements in one OR gene promoter to select that specific OR for expression. Interestingly, it was shown that the H region, which works as an LCR and is located 75 kb upstream of the MOR28 gene cluster , also contains at least one set of homeodomain-and O/E-like sites (Hirota and Mombaerts 2004). Alternatively, cis-elements and protein factors that bind to these elements could bring one given OR gene promoter to a single expression site body in the nucleus (Borst 2002;Voss et al. 2006).  However, it is important to note that other olfactory genes that are expressed in all mature olfactory neurons, such as Omp, GnaI, and AcIII, also have O/E-like sites (Fig. 3). Therefore, the mere presence of O/E-like sites in the promoter regions does not explain the mosaic pattern of OR expression in the olfactory epithelium.
In conclusion, our results indicate that intraspecies comparisons of promoter sequences are likely to be a useful strategy for identifying common regulatory motifs that may be involved in regulation of OR gene expression. A similar strategy can also be applied to other multigene families whose members are coordinately regulated, such as the pheromone receptor families (Dulac and Torello 2003).
Twenty-five-microliter nested PCR reactions were done using 1 µL of a 200-fold dilution of the primary PCR product and 25 cycles as above.

Degenerate oligonucleotide primers
The following degenerate OR primers were used in this study:

Cloning and sequencing
The RACE nested PCR products were gel-purified and cloned into the pCRII vector (Invitrogen). Colonies containing OR cDNAs were selected by colony PCR using the OR degenerate primers P26 and P27, and their orientation was determined also by colony PCR using the pair of primers T7/P27 or SP6/P27. Plasmid DNA was prepared from positive colonies using Filter Plate for high-throughput separations (Multiscreen Millipore). DNA was sequenced with the ABI PRISM Big Dye Terminator V3.1 Cycle sequencing kit using T7 or SP6 primers on an ABI PRISM 3100 Genetic Analyzer (Hitachi). Four percent of the OR sequences were truncated OR RNAs, probably because the CIP reaction during RLM-RACE was not 100% efficient.

Genomic alignment of cDNA sequences
We aligned the cDNA sequences against the mouse genome using BLAST. Only the alignments with percent identities >93% were considered. The position of each alignment was calculated, and the flanking 50-kb genomic sequences were extracted from the corresponding genomic contigs. Each sequence was realigned with its corresponding extracted genomic sequence using the Sim4 program . Only the Sim4 alignments showing average percent identity >93%, entire sequence alignment >50%, and with the best score (based on the nucleotide identity over the entire alignment) were selected. A MySQL database was loaded with the alignment information.

Clustering of cDNA sequences
The cDNA sequences were clustered based on their genomic coordinates. Sequences that share at least one same exon/intron boundary were included in the same cluster. When no exon/ intron boundaries were defined, sequences with at least 30-bp overlap in one same genomic location were included in the same cluster. The Olfactory Receptor cDNA Clusters Viewer site was generated using the Generic Genome Browser (Stein et al. 2002; http://www.gmod.org/ggb/).

Preparation of nuclear extracts
Nuclear extracts were prepared using the method described by Kudrycki et al. (1993) from olfactory epithelium dissected from 30 4-7-wk-old C57BL/6J mice. The extract was first concentrated using a Microcon centrifugal filter device (Millipore), and then the buffer was exchanged with binding buffer (10 mM Tris-HCl at pH 7.9, 1 mM EDTA, 5 mM MgCl 2 , 50 mM KCl, 10% glycerol, 3 mM DTT, 0.3 mM PMSF) using a Micro Bio-Spin P-6 chromatography column (Bio-Rad). Aliquots were stored at ‫°08מ‬C. Protein concentration was determined using the Bradford assay (Bio-Rad).

Gel shift assay
The digoxigenin (DIG) gel shift kit (Roche Applied Science) was used for gel shift assays. Binding reactions contained 2 µg of poly[d(I-C)], 0.1 µg of poly-L-lysine, 1.2 ng of labeled oligonucleotide, and 15 µg of nuclear protein extract. After a 10-min incubation on ice and a 15-min incubation at room temperature, the mixture was added with 5 µL of loading buffer (60% ‫ן52.0‬ TBE buffer, 40% glycerol, 0.2% bromophenol blue) and electrophoresed in ‫ן5.0‬ TBE on a nondenaturing 4% polyacrylamide gel in ‫ן5.0‬ TBE containing a 2-cm 15% acrylamide layer at the bottom to retain the unbound probe in the gel, as described by Bell et al. (1999). The gel was pre-electrophoresed for 1 h at 80 V before the samples were applied. Competition experiments were performed by incubating the binding reaction mixtures with a ‫ן001‬ excess of unlabeled competitor oligonucleotide for 5 min before the addition of the labeled oligonucleotide. Blotting was performed using a Bio-Rad electroblotting system, and chemiluminescence detection of the DIG-labeled DNA-protein complexes was performed using anti-digoxigenin antibody conju-gated to alkaline phosphatase and the CSPD substrate (Roche Applied Science).

INTRODUCTION
Serial analysis of gene expression (SAGE) , Massively parallel signature sequencing (MPSS)  and sequencing-by-synthesis (SBS) (www.illumina.com/pages.ilmn?ID¼201) determine the expression level of genes by measuring the frequency of sequence tags derived from polyadenylated transcripts. SAGE allows the construction of a comprehensive expression profile in which each mRNA is defined by specific 14 or 21 bp sequences adjacent to the 3 0 most NlaIII site. MPSS generates 17-20 bp sequences adjacent to the 3 0 most DpnII site from millions of mRNA molecules in a sample, providing a quantitative assessment of the transcript abundance. SBS is a new technology that is similar to MPSS (it also generates millions of 20-bp-long tags corresponding to sequences adjacent to the 3 0 most DpnII site). SAGE and MPSS have been widely used for the study of gene expression in pathological tissues (Hermeking, 2003) as well as in the study of general aspects of gene expression in different organisms (Hermeking, 2003;Jongeneel et al., 2005;Meyers et al., 2004). An essential step in the analyses of SAGE, MPSS and SBS data is the correct assignment/mapping of a tag to a gene. Basically, there are two strategies to perform this process: an annotation based on data from websites, such as SAGE Genie  and SAGEmap  or an annotation based on databases constructed using in-house computer approaches (Blackshaw et al., 2004;. SAGEmap and SAGE Genie are web-based tools that map 'tags to cDNAs' and 'tags to gene/cDNA references', respectively (these web tools are exclusive to SAGE data). These strategies for tag-togene assignments either do not produce a complete annotation or cannot batch process a large number of tags. To address these limitations, we have created a web-based application, called automatic correspondence of tags and genes (ACTG). ACTG is designed to map a large number of tags into several reference datasets. ACTG is user-friendly and generates an output file that is simple to analyze and can be integrated with other applications.

OVERVIEW OF ACTG
ACTG is a collection of several Perl scripts, Perl þ CGI scripts, shell scripts and HTML codes that performs three main tasks: (i) The assembling of virtual tag databases used in the ACTG, (ii) The uploading and processing of information submitted by the user and (iii) The mapping of the submitted tag list and generation of the output files.
ACTG datasets are composed of virtual tags (a computer prediction of tags produced by a SAGE, a MPSS or a SBS experiment) extracted from three major databases, SAGEGenie, SAGEmap and almost all public cDNAs from GenBank. The following are the main steps for the assembly of each dataset. For SAGEGenie: (i) Download of the files containing the best virtual tags matches to UniGene clusters, (ii) Parsing of the raw data, (iii) Assembly of the data and generation of the final dataset. For SAGEmap, the process is similar: (i) Download of the files containing the reliable virtual tags matches to UniGene clusters, (ii) Parsing of the raw data, *To whom correspondence should be addressed. (iii) Selection of the best tag to cDNA assignment (based on the SAGE map tag-to-gene score) and (iv) Assembly of the final dataset. For all public cDNAs sequences the process is more complex: (i) Download of all cDNAs from UniGene , RefSeq (Pruitt et al., 2005), MGC (Strausberg et al., 1999) and dbEST , (ii) Split of dbEST ESTs, RefSeq, MGC and UniGene mRNAs in subsets presenting either a poly(A) tail (at least five adenosines at the cDNA 3 0 end), or a canonical poly(A) signal (AAUAAA or AUUAAA at the 3 0 most 50 bp segment), or both, poly(A) tail and signal and (iii) Extraction of the virtual tags and assembly of the final datasets in the ACTG format (for dbEST data, only ESTs containing poly(A) tail and poly(A) tail and signal have been included). These subdivisions of the cDNA datasets, based on the 3 0 end information, are important to evaluate the tag to cDNA mapping, from which a well-defined 3 0 end produces a more reliable tag-to-gene match .
Since ACTG presents virtual tags from commonly obtained sources of cDNA data, some users may have difficulty to fully utilize the datasets and interpret the redundancy between each. We therefore created an additional virtual tag dataset, which is a non-redundant tag list that merges and removes the redundancy of every ACTG datasets. In addition, we classified the virtual tags to three categories of tag-to-gene assignment in terms of their reliability (high, medium and low), thus producing a non-redundant and ranked virtual tag dataset (details of the non redundant tag list are described on the ACTG website). These new datasets would help simplify the tag mapping process and the user's interpretation of results.
We have also created another additional data, a set of putative artifactual tags. This dataset contains tags similar (allowing 1 mismatch) to the sequence of the linker (utilized in the construction of SAGE library) and ambiguous tags, tags mapped to two or more genes. Putative tags are identified by special characters in the output files and should be considered as non-reliable in the tag-to-gene assignments. These specific features are also available in SAGE Genie . Figure 1 illustrates the main steps necessary to map a list of tags using ACTG. The first step is the submission of a file with

Select organism and tag type
Statistics of mapping Database # 1 Database # 2 All databases merged SBS tag list a list of tags (a plain text file containing the tag sequences and optional columns, for example, the tag frequencies). Second, selecting an organism (human or mouse), tag type (SAGE, MPSS or SBS), and at least one database to map the tags. The final steps include the execution of the program ('run ACTG') and downloading of the output files. ACTG produces three output files; (1) A file(s) containing the tag mapping for each selected database, (2) A file that merges all mapping results, and (3) A file containing statistics of the mapping process.
Below is a simplified guide to help best utilize the functions of ACTG and the interpretation of its mapping results. If a user is interested in: (1) Identifying known genes and produce the most 'reliable' set of tag-to-gene mapping, use databases containing virtual tags from mRNAs sequences with poly(A) tails and poly(A) tails and poly(A) signals, (2) Identifying new genes and/or new transcripts, use databases containing virtual tags from ESTs with poly(A) tails and poly(A) tails and poly(A) signals, and (3) Performing the most complete mapping, use databases containing virtual tags from SAGEGenie, SAGEmap or virtual tags from the non-redundant tag lists. In reference to the interpretation of the results, a critical aspect to be carefully inspected is the ambiguity of 'tag-to-gene' assignments. If a tag is mapped to virtual tags from two or more different cDNAs, but all sequences are in the same UniGene cluster, this 'redundancy' is acceptable and will not influence the tagto-gene match. However, the assignment is ambiguous if a tag is mapped on distinct UniGene clusters.
For further details on the functions of ACTG, the construction of the ACTG datasets, the mapping of a submitted tag list, and the interpretation of results; please visit the 'HELP & FAQ' section of the ACTG website.

APPLICATION
In order to evaluate our tool, we have selected and submitted to ACTG sets of tags from two recent publications (Blackshaw et al., 2004;Jongeneel et al., 2005). In Blackshaw et al. (2004), they were not able to map 37 813 distinct SAGE tags from 12 retina libraries. By, using ACTG, we were able to map 16 927 of these tags (databases: SAGEGenie, UniGene mRNAs and ESTs with poly(A) tail). In Jongeneel et al. (2005), they constructed an atlas of human gene expression based on tags from 32 MPSS library, where, for each library, the tags were assigned to transcribed regions (genes). For example, tags from cerebellum were mapped to 8183 genes and tags from bone marrow were mapped to 7182 genes. By using ACTG, tags from cerebellum and bone marrow were mapped to 13 773 and 12 689 genes, respectively (database: mRNAs from UniGene, RefSeq, MGC and ESTs with poly(A) tail and poly(A) signal and EST with poly(A) tail). The annotations for these datasets are available in the 'Publication' section of ACTG website as supplementary data.

SUMMARY
ACTG, a web-based tool, allows the user to quickly annotate a complete SAGE (14 or 21 bp long), MPSS (tags 17 or 20 bp long) or SBS (20 bp long) library using complete datasets for both human and mouse organisms. In addition, ACTG can filter redundant and artifact tags and generates a report of the mapping process. ACTG is a simple publicly available tag-to-gene mapping tool packaged in a user-friendly environment that addresses an essential step in SAGE, MPSS, and SBS data analyses.

INTRODUCTION
The past few years have witnessed the emergence of a picture showing a high frequency of alternative splicing in the human transcriptome (Hanke et al. 1999;. It seems that at least half of all human genes undergo alternative splicing, but the right number probably will exceed that because the detection of splicing variants of lowly expressed genes is difficult (Kan et al. 2002;Zavolan et al. 2002). The biological significance of alternative splicing is well documented for several models. A classical example is sex determination in Drosophila, in which alternative splicing acts like a genetic trigger (Baker 1989). Another classical example is the gene slo that corresponds to a calcium-activated potassium channel (Xie and Black 2001). Splice variants of this protein form a gradient in the hair cells of the inner ear of the chicken that is involved in sound detection.
Although extremely important due to their biological significance, these examples are isolated cases reflecting particular interests of individual laboratories. Due to the genomics and bioinformatics revolution, the cDNA databases remain growing at an exponential rate. This, coupled to the availability of the human genome sequence, has allowed the use of an integrated approach for the detection of splicing variants . It remains doubtful however, whether a significant fraction of the events have biological significance or whether they are spurious products from the splicing machinery (Kan et al. 2002).
There are four major types of alternative splicing: exon skipping, alternative 3Ј splice site, alternative 5Ј splice site, and intron retention. Exon skipping seems to be the most frequent type as shown by . Intron retention is certainly the least studied of all types of alternative splicing, especially because these variants are believed to be largely derived from unspliced or partially spliced pre-mRNAs. There are at least a couple of cases of intron retention events with known biological consequences. The P-element transcript in Drosophila is a transposon active in the Drosophila germ line only (Rio 1991). The blocking of the transposition in somatic cells is due to a truncation of the encoded transposase caused by the retention of the third intron in the P-element transcript. The Ret tyrosine kinase also presents a couple of intron retention events, one of them generates a truncated protein that is enriched in familial and sporadic pheochromocytomas (Le .
Here we present a large-scale analysis of intron retention for a set of 21,106 known human genes. We evaluated the distribution of intron retention events in relation to several features such as GC content, codon usage, coding potential, and conservation in mouse among others. The observed data suggest that a significant fraction of the intron retention events is not spurious and probably reflects biological significance.

Identification of intron retention events
Because of the large-scale nature of their generation, EST data are enriched with artifactual sequences like genomic contamination and unprocessed or partially processed mRNAs. These issues are critical for the detection of true splicing variants. The availability of an almost complete draft of the genome sequence has allowed the integration of both genomic and cDNA data to identify splicing variants . We have used a similar strategy by mapping all human cDNAs onto the draft sequence of the human genome . We have also developed a clustering strategy that groups all cDNAs mapping to the same genomic region into a single cluster (see Materials and Methods for more details). A nonredundant set of 21,106 different cDNA clusters containing at least one known human full-insert cDNA (by that we mean a fully sequenced cDNA clone) was used in the analyses described below. To validate our database, we performed an analysis of exon skipping and compared our results to data previously published by different groups. We found, for instance, that 52% of all known human genes undergo exon skipping and that 77% of the events are located within the coding region. Similar numbers were found by different groups (Hanke et al. 1999;Kan et al. 2002).
To be identified as a cDNA containing an intron retention event, a sequence was required to span the entire length of the respective retained intron. Furthermore, the same sequence had to contain at least one exon/exon boundary in order to minimize the chance of the retention being an artifact. This critical step eliminates unprocessed messages (although there is still the possibility of finding partially processed messages). The identification of an intron retention event was based on a pairwise comparison of a transcript to other sequences belonging to the same cDNA cluster. A retention event could be identified in different possible situations (Fig. 1). The prototype used to characterize the retention, spanning the respective exon/exon boundary, could be either an EST or a full-insert cDNA. The sequence containing the retained intron could also be either an EST or a full-insert cDNA. We found 3127 known genes with at least one intron retention event, which corresponds to 14.8% of the total number of known genes analyzed here. Table 1 shows the numbers of genes found for every possible case except the cases where both the prototype and the sequence containing the retention are ESTs. We decided to exclude those cases from our analysis because they involve only ESTs, which are known to be low-quality sequences. The majority of cases correspond to those pairs where the prototype is an EST and the sequence showing the retention is a full-insert cDNA. This is probably due to at least two factors: the average length of EST sequences and a bias for longer clones in the generation of full-insert cDNA sequences. A complete list of all GenBank accessions reporting an event of intron retention is available at http://www. compbio.ludwig.org.br/∼pgalante/IR. Four genes reporting an intron retention event were selected for experimental validation. RT-PCR reactions were performed using primers flanking the intron being retained and primers specific to the retained sequence. Figure 2 shows the products of amplification for one of the four genes in all tissues sampled. We observed the variant containing the retained intron for all genes as well as the corresponding prototypes.

Distribution of events along transcripts
Because the presence of an intronic sequence is expected to cause a dramatic effect on the corresponding protein when it occurs within the coding sequence (CDS), we evaluated the CDS/untranslated region (UTR) distribution of all events within the data set of 640 genes where both the prototype and the sequence defining the retention event were fullinsert cDNAs. For the sake of simplicity, from now on we will call these pairs of sequences the "elite group." Because the distribution of events is primarily dependent on the frequency of introns in each part of the sequence, we first obtained the number of introns in the 5Ј UTR, CDS, and 3Ј UTR. This analysis was done using a set of 16,951 full-insert cDNA sequences that contained annotated 5Ј and 3Ј UTRs. Based on that, we calculated the number of expected intron retention events in each of the three sections of the transcript where IR is the total number of retention events, Ii is the number of introns in the respective transcript section, and I is the total number of introns in the data set, and compared that with the observed numbers of events ( Table 2). We noticed a significant excess of intron retention events within both UTRs, with a corresponding decrease of observed events within the CDS ( 2 = 229, 1 degree of freedom, p < 10 −49 ).
It is well known that the presence of premature stop codons in mRNA triggers its degradation through a process of nonsense-mediated decay (NMD; . It is, therefore, possible that the observed biased distribution for retention events in the UTRs is simply a product of NMD. This would actually suggest that the rate of intron retention is higher than observed and a fraction of the events is filtered off by NMD. Another possibility is underreporting of sequences containing a premature stop codon. To test whether this was a factor affecting the above described distribution, we selected all cases of intron retention reported by the Mammalian Gene Collection (MGC) initiative, which reports sequences to GenBank even when they contain premature stop codons . Intron retention events reported in MGC sequences are still biased for UTRs, showing that underreporting does not seem to affect the above distribution (Table 2). Those cases confirmed by at least two cDNA sequences are shown in parentheses. "Total" corresponds to the nonredundant sum of the values in the corresponding row or column. The redundance is due to the fact that the same gene may have intron retention events reported by both ESTs and full-insert cDNAs (in those cases, the gene was counted only once).

FIGURE 2. (A)
RT-PCR validation of a representative cDNA (BC004239) showing intron retention. When primers flanking the retained intron are used, two bands are visible corresponding to both variants (even numbers). If a primer specific to the intronic sequence is used, a single band corresponding to the variant with the retention is observed (odd numbers). cDNAs from colon, brain, prostate, and lung were used in the reactions. Template was excluded in the control reactions. (B) Experimental validation for all remaining genes in all tissues. (+) The retained intron was observed in the corresponding tissue. A gel containing all validation data is presented as Supplemental Figure 2, which can be found at http://www.compbio.ludwig. org.br/∼pgalante/IR/. We considered only those cases in which both sequences correspond to full-insert cDNAs. For the statistical analysis we grouped the 5Ј and 3Ј UTRs. The excess of retention in both UTRs is statistically significant: elite group p = 10 −49 , MGC (Mammalian Gene Collection) p = 10 −17 . Expectation was calculated based on the density of introns in the respective transcript section (IR + Ii/I, where IR is the total number of retention events, Ii is the number of introns in the respective transcript section and I is the total number of introns in the dataset).

Retained introns code for protein domains
One way to assess the biological significance of an intron retention event is to evaluate its contribution to the protein domains encoded by the respective fullinsert cDNA. One hundred forty-seven full-insert cDNAs from the elite group containing at least one retention event totally located in the coding region were searched against Pfam (Bateman et al. 2002). See Materials and Methods for further details. Figure 3 illustrates the different types of contribution by a given retained intron to the coding of a protein domain. Only two domains were encoded entirely by two different retained introns, out of 147. Retained introns are less efficient in coding for protein domains than exons (30 exons, out of 830, encode an entire domain in the same set of full-insert sequences).
There is no statistical difference between the numbers obtained for retained introns and exons (Fisher's exact test, p = 0.21). However, retained introns are more efficient than nonretained introns present in the same set of full-insert cDNAs (no nonretained intron, out of 785, encodes an entire domain). This difference is statistically significant (Fisher's exact test, p < 0.026). Furthermore, if the retained intron was joined to the two flanking exons (exon/retained intron/exon), the number of domains entirely encoded by such segments increased to 15 (a hit was considered significant when at least 30% of the domain was encoded by the retained intron). For the 13 cases in which the intron partially encodes a domain, we tested if the exclusion of the sequence corresponding to the retained intron affected the ability to detect the correct domain. We performed a manual inspection of the protein domains encoded by a modified version of the intron retention variant (without the amino acids encoded by the retained intron). In six of these sequences, the domain encoded with the contribution of the retained intron was not found and in six others the domain was found with scores considerably lower and with e-values increased by orders of magnitude. In only one case could we not detect a substantial change in the scores after deleting the intronic sequence. This observation suggests that the amino acids encoded by the retained introns are a relevant part of a structural unit.
In addition, we found 24 other cases where the domain was longer than the "exon/retained intron/exon" unit (i.e., the domain was only partially coded by this unit) with a significant contribution of the intronic sequence (at least 30% of the domain). Thus, 39 out of 147 retained introns (26%) participated in the coding of a protein domain.
As an additional control, we tested if the retained introns located in the UTRs would participate in the coding of protein domains. None of the 762 amino acid sequences (254 UTR-retained introns translated in three phases) showed a complete or partial match against a domain in Pfam (Fisher's exact test, p < 10 −17 when compared to 39 cases in 147 retained introns in the CDS).

Selection on intron sequences involved in retention events
To evaluate whether selection for coding potential was acting on intron sequences involved in retention events, we analyzed those introns from our elite group located within the CDS (287 introns as shown in Table 2). We first compared the frequencies of all four bases in the set of introns involved in retention events to the set containing all nonretained introns. Because it has been shown that shorter introns have a higher GC content in humans (Lander et al. 2001), we divided both data sets into length categories. A significant difference ( 2 = 2276, 6 degrees of freedom, p < 10 −50 ) is observed between retained (higher GC content) and nonretained introns (Fig. 4). Figure 4 also shows that the GC distribution for retained introns is similar to the GC distribution for their flanking exons (p = 0.23, 2 = 7.9, 6 degrees of freedom). The distribution observed for both exons and nonretained introns is in accordance with Lander et al. (2001).
The frequency of stop codons in retained introns is lower than expected. Based on the nucleotide frequencies observed for retained introns located within the CDS, the expected number of stop codons is 1064 in 26,941 codons (see Materials and Methods). The observed number of stop codons (651) was significantly smaller than expected ( 2 = 103, 1 degree of freedom, p < 10 −25 ). We found 88 full-insert cDNAs, out of 287, in which the retention event generated a putative truncated protein due to the presence of a premature stop codon within the retained intron (available at http://www.compbio.ludwig.org.br/∼pgalante/IR). The identification of 88 cases in which the retention generates a putative truncated protein allowed us to examine the GC content within the retained intron for sequences upstream of and downstream from the premature stop codon. It is expected that selection would be relaxed on intronic sequences downstream from the premature stop codon. The GC content for sequences upstream of and downstream from the premature stop codon is 58% and 49%, respectively.

Codon usage in retained introns, exons, and nonretained introns
Codon usage in retained introns, exons, and nonretained introns was evaluated as described in Materials and Methods. Two strategies were used to compare the different data sets. The first strategy compares the whole codon usage table (61 codons) of the three data sets (results are labeled A in Fig. 5). The second strategy is based on the calculation of the average number of amino acids using the same codons in the three different data sets (results are labeled B in Fig. 5). Both approaches show that retained introns are more similar to exons than to nonretained introns.

Conservation of intron retention in mouse cDNA sequences
To assess if events of intron retention found in the human transcriptome are conserved in the mouse transcriptome we aligned human full-insert cDNAs presenting intron retention (elite group with at least one mouse cDNA hit) to mouse full-insert cDNAs and ESTs using BLAST 2.0  with an e-value threshold of e −10 . We found that 38%-57% of all retained introns (depending on the UTR/CDS localization) present a mouse hit (Table 3, first line). The sequence identity of orthologous retained introns is 84%, significantly higher than the average identity found for orthologous nonretained introns (60%) and similar to the average degree of identity found for orthologous exons (87%). These numbers are similar to those observed by . Interestingly, we found only three cases of retained intron in mouse with a premature stop codon, suggesting that the putative retention in mouse preserved the open reading frame.
To verify if the mouse cDNA also corresponds to an intron retention variant, we searched for at least one other cDNA that could define the corresponding intron at the genome level. For many examples, the region containing the retained intron was not covered by any other cDNA sequence, thus not allowing the identification of an intron in that region (second line of Table 3). An intron retention in human is deemed "conserved" in mouse if a mouse hit covers the retained intron in human and if there is also evidence that the retained sequence is also an intron in the mouse genome. We found 10 cases, out of 46, of conserved intron retention. This number is certainly an underestimate  because the prototype may not have been sequenced and all cDNA sequences covering that region may correspond to the variant showing the retention, as also pointed out by Thanaraj et al. (2003).

DISCUSSION
One could expect that the splicing process has an intrinsic error rate that is difficult to estimate. It is reasonable to suppose that at least a fraction of all splicing variants found in the cDNA databases correspond to artifacts caused by such errors. This is even more pronounced nowadays as databases are reaching a critical in-depth sampling. Thus, an important question refers to the amount of artifactual variants identified by the use of cDNA databases. Although valid for all types of alternative splicing, this seems to be more critical for intron retention, because unprocessed or partially processed messages would be easily detected as variants. Kan et al. (2002) found an intron retention rate of around 35%. They then used some statistics to infer the reliability of a given splicing variant. They found, for instance, that less than 5% of all genes exhibited intron retention at a 95% confidence interval (p < 0.05). We have found that 14.8% of all genes sampled here exhibited putative events of intron retention. This discrepancy may be due to the fact that we have not applied any statistics in our data set and/or to differences in the data sets used. On the other hand, we considered only sequences showing evidence of splicing to exclude totally unprocessed transcripts, and most of our analyses were performed on an elite group. Interestingly, if we consider the elite group only, we find evidence of intron retention in 4.6% of the whole data set used (640 genes with intron retention divided by 13,841 genes with at least two full-insert cDNAs).
An important aspect of our analysis is the apparent nonrandom distribution of intron retention events regarding some specific features. The rate of intron retention is more pronounced in the UTRs as clearly shown in Table 2. There is a threefold and 14-fold excess (above the expectation) of events in the 5Ј and 3Ј UTRs, respectively. Much of the bias is probably due to NMD that would trigger the degradation of transcripts bearing premature stop codons. Although this would explain the lower number of events within the CDS, it does not explain the different rates of intron retention between the 5Ј and 3Ј UTRs. This bias toward the 3Ј UTR could be explained by the direction of RNA synthesis, because splicing seems to be coupled to transcription (for review, see  and therefore the last intron would have less time available to be spliced out. However, for 90% of the retained introns, we were able to find at least one exon/exon boundary downstream from the retained intron, suggesting that the last intron was spliced as efficiently as the first introns in the transcript. One possible explanation implying biological significance is an effect of the retained intron on mRNA stability. Cis-acting elements present in the 3Ј UTR are known to alter the stability of mRNAs (Bashirulla et al. 2001), making them either more or less stable. The presence of such elements in the retained introns would affect the stability of that variant transcript.  observed that retained introns are significantly shorter than nonretained introns. We observed the same trend in our data set (see Supplemental Fig. 1 at http://www.compbio.ludwig.org.br/∼pgalante/IR/). However, the interpretation that this feature is related to biological significance is not straightforward because the retention of very long introns is unlikely to occur. Furthermore, other factors like cloning efficiency could contribute to the biased length distribution observed for retained introns.
Some analyzed features show that the retained introns are considerably more similar to exons than to other introns. The GC content of such sequences is higher than that of nonretained introns in the same size range and is similar to exons. Introns with a higher GC content might have a lower excision rate as demonstrated by Goodall and Filipowics (1991) in dicots. This phenomenon, however, does not seem to exist for other organisms, including vertebrates (Goodall and Filipowics 1989). The codon usage follows the same trend and is overall much more similar to that of exons than to introns. Likewise, retained introns encode more domains than introns but less than exons, showing that the intronic sequence is a relevant part of a larger functional unit, even when it is not coding the entire domain.
These aspects show that the pattern emerging from the data set of variants is not random and therefore can be We performed the analysis separately for CDS and UTR. The first row corresponds to the results of searching the retained introns in human against mouse cDNAs (ESTs and full-insert cDNA). The second row defines the cases that have more than one mouse cDNA in the region of the retained intron (these cases are informative for seeking an intron at the genome level). The third row corresponds to those cases defined in the second row where an intron is observed at the genome level. See text for more details.
interpreted as having some biological significance. We cannot rule out, however, that unknown properties of the splicing process would generate such a bias in the pattern of intron retention events. We have also found evidence that retained introns within the CDS are under selective pressure for coding potential. Not only is the GC content different from that of introns, as discussed above, but also the frequency of stop codons is considerably lower than random expectation. It should be noted, however, that the reduced frequency of stop codons can also be explained by NMD. Other evidence of selective pressure is that a fraction of introns being retained in human is also retained in mouse, presenting an average sequence identity comparable to orthologous exons (83%). As observed by Thanaraj et al. (2003) the presence of a same variant in two species is strongly affected by the cDNA representation of that gene in both transcriptomes. Therefore, it is reasonable to assume that the observed rate of cooccurrence is an underestimate.
As noted originally by , and recently stressed by , new variants are expected to correspond to a small fraction of all transcripts from a given gene. Mutations that affect splicing can allow the production of new proteins without the loss of the original one. If these new splicing variants are spurious products of the splicing reaction, filtering selection would be relaxed on the machinery that generates them, because they correspond to a small fraction of all transcripts and the original version is maintained. On the other hand, if the new variants have some biological significance, it is expected that selection will act on them to maintain the new biological function. Thus, the observed signals of selection support a biological role for the retained introns.
We have, however, identified 88 cases in which the intron retention event generates a putative truncated protein. For 40% of these cases the premature stop codon was still located in the 3Ј-most exon, which can explain why they were not filtered out by NMD. Among these 88 cases there are genes that seem to be involved in several types of syndromes such as Williams Beuren syndrome, chromosome 22 region, and the Betten, Spielmeyer-Vogt disease gene. Furthermore, there are several genes with retained introns clearly related to the tumorigenic process like p19A, kallikrein 3 and 4, TNF receptor, BCL2-like 11, and CDC2-like 10, among others. Further analyses are needed to explore the possibility of an association between the splicing variants and any feature of these diseases.
For those cases showing a putative premature stop codon, the GC content is higher in the intronic sequence upstream of the stop codon (58%) than in the downstream sequences (49%). This reinforces the notion that selection is acting on the retained introns.
It is important to stress that our results do not exclude the possibility that a fraction of intron retention events is indeed artifactual. The data indicate, however, that this fraction is not sufficiently large to eliminate signals that can be interpreted as products of events with biological significance. We make available subsets of retained introns that are more likely to be free of artifacts (www.compbio.ludwig. org.br/∼pgalante/IR). These include, for instance, a set of 385 cases of intron retention that are confirmed by at least two cDNA sequences. These sets might be useful for the identification of putative regulatory elements.
The definition of introns and exons may become fuzzy when we take into account intron retention. One could establish two requirements for a genomic sequence to be considered an exon: (1) It should be present in the mature mRNA, and (2) it should be flanked by sequences excised at the RNA level. Because the retained introns identified in this study (which seem to be under selection) violate the second requirement, they should be considered introns subjected to regulated alternative splicing. This raises the question of whether retained introns are portions of exons becoming introns or introns being incorporated into exons. The analysis performed here cannot answer this question. Expression data coupled to more complex interorganism comparisons could help to elucidate this issue.
Finally, regarding the mechanistic nature of the observed retention events, it is possible that the flanking exons or the retained introns themselves bear weak splice sites and/or do not contain proper splicing regulatory elements. An analysis of splice sites and regulatory elements currently underway will hopefully shed some light on this subject.

cDNA mapping and clustering
All human cDNAs available in dbEST ( July 2002) and mRNA sequences from known human genes from Uni-Gene release 153 (Schuler et al. 1996) were aligned to the masked human genome sequence (build 29, obtained from NCBI) by using pp-Blast , an implementation of MEGA-BLAST  for a parallel cluster. The parameters used in MEGABLAST were: -f T -J F -F F -W 24. The MEGA-BLAST output was parsed and a MySQL database was loaded with the mapping information. Spurious hits were excluded from the mapping database by using an additional set of alignment criteria. These include a minimum degree of identity for a cDNA/genome alignment set to 93% over at least 45% of the total EST length or 55% of the total length of the full-insert sequence. Furthermore, sequences mapping to more than one location on the genome were given a score for alignment quality. Higher score was associated with a higher identity over a longer alignment. Only the sequences with the highest scores were used in the analyses reported here. Clustering of cDNA sequences was based on their genomic coordinates as described by . Briefly, if two sequences shared at least partially the same gene structure they were joined into the same cluster. If no exon/intron boundary was defined, a sequence had to have at least a 100-bp overlap with another sequence at the genome level to be added to the respective cluster.

Identification of intron retention events
A set of 21,106 distinct cDNA clusters containing at least one known full-insert cDNA sequence was used in the present analysis. Intron retention events were characterized by the presence of at least one cDNA sequence in which the corresponding intron was not spliced out. We considered valid only those cases in which the sequence characterizing the retention also defined at least one more intron that was not retained. This excludes from our analysis artifacts derived from unprocessed messages.

Distribution of events along transcripts
We selected from the 21,106 full-insert cDNAs those sequences (16,951) containing an annotated CDS. We then mapped all retention events to these sequences. For this analysis, we have only used those cases in which both the prototype and the sequence containing the retained intron were full-insert cDNAs to avoid any positional bias that could be introduced by using 3Ј and 5Ј ESTs. A retention was considered to be located within a given transcript section if the 5Ј end of the intron was mapped to that section. In a few cases in which the presence of the intron generated a stop codon, the 3Ј end of the intron would be located in the 3Ј UTR whereas the 5Ј end would be located in the CDS.

Codon usage in nonretained introns, exons, and retained introns
Codon usage was determined for those cases belonging to the elite group. The exons and nonretained introns used as a control came also from the elite group to avoid any possible bias in the selection of genes. In the case of both retained and nonretained introns, the frame was defined by the frame of the upstream exon. To evaluate the similarity in codon usage in the three sets we made two distinct measurements. The first was a comparison of the entire codon usage table (61 codons, 60 degrees of freedom) among the three sets. We performed 1000 pairwise comparisons of the random distributions of 8216 codons (the number of codons in the set of retained introns) along the 61 possible codons to calculate the average 2 values (labeled A in Fig. 5) and standard deviation. Internal variations within each set (curved arrows in Fig. 5) were also calculated as an internal control. The second measurement refers to the number of amino acids that presented similar codon usage ( 2 tests at p Յ 0.05; labeled B in Fig. 5). As methionine and tryptophan have only one codon, they were excluded from this last analysis. Again 1000 random sets were built from each set and pairwise compared, grouped by amino acid. The numbers of amino acids with p values <0.05 in 2 tests were counted and presented in Figure 5.

Analysis of protein domains encoded by retained introns
Amino acid sequences of full-insert cDNAs, belonging to the elite group, with a retained intron located entirely in the CDS were submitted to domain search in Pfam 9.0 (Bateman et al. 2002) using the program hmm in a GeneMatcher hardware (Paracel Inc). Retained introns smaller than 99 bp and longer than 990 bp were excluded from this analysis. The final number of sequences submitted to Pfam was 147. To be considered a significant hit, the retained intron had to encode at least 30% of the protein domain and have an e-value <10 −1 . The coordinates of the domains were then converted to their cDNA positions. The cDNA coordinates of retained introns or exons (control set) entirely in the CDS were then compared to domain coordinates and grouped in those matching a whole domain or part of it. Domains were considered different in a given protein sequence when they presented distinct annotation and coordinates. For the comparison of exons, introns, and retained introns, domain search was performed with each exon or intron (of the full-insert cDNAs with intron retention) being a separate query. Only sequences with similar lengths were analyzed (<330 and >33 amino acids).

Experimental validation
The four genes tested by RT-PCR are: (1) BC015569-similar to SRp25 nuclear protein, (2) BC004239-jumping translocation breakpoint, (3) BC025673-similar to transcription elongation factor B (SIII), and (4) BC011903-similar to the small subunit of calpain 4. Reverse transcription was carried out using the Superscript First Strand Synthesis Kit, according to the manufacturer's instructions (Life Technologies). The cDNA was synthesized by incubating 2 µg of either brain, lung, prostate, or breast total RNA in 20 µL of 1× reverse transcriptase buffer containing 0.5 µg oligo dT primer, 0.5 mM dNTP, 5.0 mM MgCl 2 , 10 µM DTT, 80 U ribonuclease inhibitor, and 200 U Superscript II reverse transcriptase at 42°C for 1 h. The quality of total RNA was tested by PCR using MLH1 primers located at intronic sequences flanking exon 12 (forward: 5Ј-TGGTGTCTCTAGTTCTGG-3Ј, and reverse: 5Ј-CATTGTTGTAGTAGCTCTGC-3Ј), as an indicator of possible genomic DNA contamination. Three specific primers were designed for each intron retention event: two primers annealing to the flanking exons of the retained intron and one additional primer annealing to the retained intron itself. RT-PCR reactions were carried out in a 25-µL reaction mixture containing 1 µL of cDNA, 1× Taq DNA polymerase buffer, 1 µL of dNTP 2.5 mM, 6 pmoles of primers, 0.5 µL of MgCl 2 50 mM, and 1 unit Taq DNA polymerase (GIBCO/BRL). Standard PCR conditions were: 4 min at 94°C (initial denaturation), 45 sec at 94°C, 45 sec at 58°C, and 1 min at 72°C for 40 cycles and a final extension step of 10 min at 72°C. PCR products were analyzed both on 8% silver-stained polyacrylamide gels and on ethidium bromide agarose gel.
When evaluated globally, several features related to the distribution of NATs strongly suggest they have a prominent role in antisense regulation in gene expression [4][5][6][7]. For instance, expression of S-AS transcripts tends to be positively or negatively correlated and is more evolutionarily conserved than expected by chance [4,5,7]. Although experimental validation of a putative regulatory role has been achieved for a few models [8][9][10], it is still unknown whether antisense regulation is a rule or an exception in the human genome. NATs have been implicated in RNA and translational interference [11], genomic imprinting [12], transcriptional interference [13], Xinactivation [14], alternative splicing [10,15] and RNA editing [16]. Moreover, an accumulating body of evidence suggests that NATs might have a pivotal role in a range of human diseases [2].
NATs were initially identified in studies looking at individual genes. However, with the accumulation of whole genome and expressed sequences (mRNA and ESTs) in public databases, a significant number of NATs has been identified using computational analysis [17][18][19][20][21][22]. These studies showed a widespread occurrence of these transcripts in mammalian genomes. The first evidence that antisense transcription is a common feature of mammalian genomes came from analysis of reverse complementarity between all available mRNA sequences [17]. Subsequent studies, using larger collections of mRNA sequences, ESTs and genomic sequences, confirmed and extended these initial observations [18][19][20][21][22]. More recently, other sources of expression data, such as serial analysis of gene expression (SAGE) tags, were used to expand the catalog of NATs present in mammalian genomes [23,24]. At present, it is estimated that at least 15% and 20% of mouse and human transcripts, respectively, might form S-AS pairs [18,22], although a recent analysis [25] reported that 47% of human transcriptional units are involved in S-AS pairing (24.7% and 22.7% corresponding to S-AS pairs with exon and non-exon overlapping, respectively).
The major obstacle in using expressed sequence data for NAT identification is how to determine the correct orientation of the sequences, especially ESTs. Many ESTs were not directionally cloned and even well-known mRNA sequences were registered from both strands of cloned cDNAs or are incorrectly annotated. As done by others [18,22,23], we here established a set of stringent criteria, including the orientation of splicing sites, the presence of poly-A signal and tail as well as sequence annotation, to determine the correct orientation of each transcript relative to the genomic sequence and made a deep survey of NAT distribution in the human and mouse genomes. Using a set of computational and experimental procedures, we extensively explored expressed sequences and massively parallel signature sequencing (MPSS) data mapped onto the human and mouse genomes. Besides generating a catalog of known and new S-AS pairs, our analyses shed some light on functional and evolutionary aspects of S-AS pairs in mammalian genomes.

Overall distribution of S-AS pairs in human and mouse genomes
To identify transcripts that derive from opposite strands of the same locus, we used a modified version of an in-house knowledgebase previously described for humans [26][27][28]. This knowledgebase contains more than 6 million expressed sequences mapped onto the human genome sequence and clustered in approximately 111,000 groups. Furthermore, SAGE [29] and MPSS [30] tags were also annotated with all associated information, such as tag frequency, library source and tag-to-gene-assignment (using a strategy developed by us for SAGE Genie [31]). An equivalent knowledgebase was built for the mouse genome (for more details see Materials and methods).
We first designed software that searched the human and mouse genomes extracting gene information from transcripts mapped onto opposite strands of the same locus. Several parameters were used by the software to identify S-AS pairs, such as: sequence orientation given by the respective Gen-Bank entry; presence and orientation of splice site consensus; and presence of a poly-A tail (for more details see Materials and methods). We found 3,113 and 2,599 S-AS pairs in human and mouse genomes, respectively, containing at least one fullinsert cDNA (sequences annotated as 'mRNA' in GenBank and referred to here as such) in each orientation (Table 1). Furthermore, we also made use of EST data from both species. A critical issue when using ESTs is the orientation of the sequence, a feature not always available in the respective GenBank entries. We overcame this problem by simply using those ESTs that had a poly-A tail or spanned an intron and, therefore, disclosed their strand of origin by the orientation of a splicing consensus sequence (GT...AG rule). We found 6,964 and 5,492 additional S-AS pairs when EST data were incorporated into the analysis, totaling 10,077 and 8,091 pairs for human and mouse genomes, respectively ( Table 1). All of these pairs contained at least one mRNA since we did not analyze EST/EST pairs. It is important to note that we haven't considered in the present analysis non-polyadenylated transcripts and trans-NATs. Thus, the total number of NATs is likely to be even higher in both genomes. Data presented in Table 1 are split in cases where a single S-AS pair is present in a given locus (single bidirectional transcription) and in cases where more than one pair is present per locus (multiple bidirectional transcription). Additional data file 1 lists two representative GenBank entries for all S-AS pairs split by chromosome mapping in the two species. As previously observed [17], S-AS pairs are under-represented in the sex chromosomes of both species (Additional data file 2).
The above numbers confirm that S-AS pairs are much more frequent in mammalian genomes than originally estimated [4,17,18]. Our analyses suggest that at least 21,000 human and 16,000 mouse genes are involved in S-AS pairing. These numbers are more in agreement with those from [32]  analysis using tiling microarrays to evaluate gene expression of a fraction of the human genome. For the mouse genome, our numbers are in agreement with those reported by Katayama et al. [8]. A more recent analysis [25] also gives a similar estimate of S-AS pairs in both human and mouse genomes.
Could this high number of S-AS pairs be due to the stringency of our clustering strategy? If the same transcriptional unit is fragmented in close contigs due to 3' untranslated region (UTR) heterogeneity, the total number of clusters would be inflated, leading to an erroneous count of S-AS pairs. To evaluate this possibility, we relaxed our clustering parameters, requiring a minimum of 1 base-pair (bp) same strand overlap for clustering. Furthermore, we collapsed into a single cluster all pairs of clusters located in the same strand and less than 30 bp away from each other. Additional data file 3 shows the total number of clusters and S-AS pairs after this new clustering strategy was employed. As expected, both the total number of clusters and S-AS pairs decreased with the new clustering methodology. The total number of clusters decreased by 2% and 1% for human and mouse, respectively, while the total number of S-AS pairs decreased by 0.3% for both human and mouse. Thus, the small difference observed does not affect the conclusions on the genomic organization of S-AS pairs. For all further analyses, we decided to use the original dataset obtained with a more stringent clustering methodology.
We further explored the genomic organization of S-AS pairs using the subset of 3,113 human and 2,599 mouse pairs that contained mRNAs in both sense and antisense orientations. The genomic organization of S-AS pairs can be further divided into three subtypes based on their overlapping patterns: head-head (5'5'), tail-tail (3'3') or embedded (one gene contained entirely within the other) pairs (Table 2). For a schematic view of the genomic organization of S-AS pairs, see Additional data file 4. Embedded pairs are more frequent in both species, corresponding to 47.8% and 42.5% of all pairs in human and mouse, respectively. If we take into account the intron/exon organization of both genes, we observe that the most frequent overlap involves at least one exon-intron border. In spite of this, a significant amount of NATs maps completely within introns from the sense gene in both human and mouse (category 'Fully intronic' in Table 2). Interestingly, more than three-quarters of all S-AS pairs categorized as 'Fully intronic' fall within the embedded category for human and mouse. How unique is this distribution? Monte Carlo simulations, in which we randomly replaced NATs in relation to sense genes while keeping their 5'5'/embedded/3'3' orientation, show that the distribution of S-AS pairs is quite unique. All three categories of S-AS pairs deviate from a random distribution (chi-square = 11.5, df (degrees of freedom) = 2, p = 0.003 for embedded pairs; chi-square = 49, df = 2, p = 2.3 × 10 -11 for 5'5' pairs; chi-square = 132, df = 2, p = 2.1 × 10 -29 for 3'3' pairs). This peculiar distribution will be further discussed in the light of the expression analyses. Since these intronic NATs have been shown to be over-expressed in prostate tumors [33], our dataset should be further explored regarding differential expression in cancer. Due to their genomic distribution, any putative regulatory role of these intronic NATs would have to be restricted to the nucleus. Interestingly, Kiyosawa et al. [34] observed that a significant amount of NATs in mouse is poly-A negative and nuclear localized.
Another interesting observation is the higher frequency of intronless genes within the set of S-AS pairs (Table 3). About half (47%) of all mRNA/mRNA S-AS pairs in humans contains at least one intronless gene. This number is slightly lower for mouse (44%) ( Table 3). Interestingly, intronless genes are significantly enriched within the set of embedded pairs (chi-square = 95.9, p < 1.2 × 10 -22 for human and chisquare = 3.98 and p < 0.045 for mouse). For humans, 66% of all S-AS pairs containing at least one intronless gene are within the 'embedded' category; Sun et al. [5] found 43.4% of their S-AS pairs as 'embedded'. Furthermore, they found 35% of 3'3' pairs while we found only 25%. These differences are probably due to the fact that Sun et al. [5] included in their analyses pairs containing only ESTs.
All these results clearly show that subsets of S-AS pairs have distinct genomic organization, suggesting that they may play different biological roles in mammalian genomes. Below we will discuss these data in a functional/evolutionary context. Single bidirectional transcription corresponds to those loci in which only one S-AS pair is present. Multiple bidirectional transcription corresponds to those loci in which more than one S-AS pairs is present (at least one gene belongs to more than one S-AS pair).

Conservation of S-AS pairs between human and mouse
Using our set of human and mouse S-AS pairs, we measured the degree of conservation between S-AS pairs from human and mouse. Since the numbers reported so far are discrepant, ranging from a few hundred [5,6] to almost a thousand [25], we decided to use different strategies. We first used a strategy based on HomoloGene [35]. The number of S-AS pairs with both genes mapped to HomoloGene is 854 for human and 579 for mouse. Among these, 190 S-AS pairs are conserved between human and mouse. One problem with this type of analysis lies in its dependence on HomoloGene, which, for example, does not take into consideration genes that do not code for proteins. Therefore, we decided to implement a different strategy, in which we identified those pairs that had at least one conserved gene mapped by HomoloGene and tested each known gene's NAT for sequence level conservation. Using this strategy, we found an additional 546 cases, giving a total of 736 (190 + 546) conserved S-AS pairs between human and mouse. Finally, we also applied to our dataset the same strategy used by Engstrom et al. [25], in which they counted the number of human and mouse S-AS pairs that had exon overlap in corresponding positions in a BLASTZ alignment of the two genomes. We applied the same strategy to our dataset and found 1,136 and 1,144 corresponding S-AS pairs in human and mouse, respectively. As observed by Engstrom et al. [25] the numbers from human and mouse slightly differ because a small proportion of mouse pairs corresponded to several human pairs and vice versa. Additional data file 5 lists all S-AS pairs found by the three methodologies discussed above.
There is a predominance of 3'3' pairs in all sets of conserved S-AS pairs. For the first strategy solely based on Homolo-Gene, 67% of all pairs are 3'3' compared to 19% embedded and 14% 5'5'. For the dataset obtained using the strategy from Engstrom et al. [25], there is also a prevalence of 3'3'pairs (48%) compared to embedded (14%) and 5'5 (38%) pairs. We have also modified the method of Engstrom et al. [25] to take into account all S-AS pairs and not only those presenting exon-exon overlap. These data are shown in Additional Data File 6. We observed that S-AS pairs whose overlap is classified as 'Fully intronic' are less represented in the set of conserved S-AS pairs (18% in this set compared to 29% in the whole dataset of S-AS pairs). The same is true for S-AS pairs containing at least one intronless gene (26% in the set of conserved S-AS pairs compared to 47% in the whole dataset). These last results are in accordance with our previous observation that conserved S-AS pairs are enriched with 3'3' pairs. As seen in Tables 2 and 3, 3'3' pairs are poorly represented in the categories 'Fully intronic' (Table 2) and 'Intron/intronless' ( Table 3).

Discovery of new S-AS pairs in human and mouse genomes using MPSS data
Large-scale expression profiling tools have been used to discover and analyze the co-expression of S-AS pairs [5,23,34]. Quéré et al. [23], for instance, recently explored the SAGE  . Our strategy was based on the generation of virtual tags from each genome by simply searching the respective genome sequence for DpnII sites. Since these sites are palindromes, we extract, for each one, two virtual tags (13 and 16 nucleotide long tags for human and mouse, respectively), both immediately downstream of the restriction site but in opposite orientations (see Materials and methods for more details). In this way, we could evaluate the expression of transcriptional units present in both strands of DNA. We obtained 5,580,158 and 8,645,994 virtual tags for the human and mouse genomes, respectively. This set of virtual tags was then compared to a list of tags observed in the MPSS libraries. As true for any study using mapped tags, our analysis misses those cases in which a tag maps exactly at an exon/exon border at the cDNA level.
We first evaluated the number of cDNA-based S-AS pairs (shown in Table 1) that were further confirmed by the presence of an MPSS tag. Data for this analysis are presented as Additional data file 8. Roughly, 84% and 51% of all cDNAbased S-AS pairs were confirmed by MPSS data for human and mouse, respectively.
Since we were interested in finding new antisense transcripts, we searched for tags found in the MPSS libraries that were mapped on the opposite strand of both introns and exons of known genes. For this analysis we excluded those genes that were already part of S-AS pairs as described above. For humans, 4,308 genes have at least one MPSS tag derived from the antisense strand (Table 4). For 1,221 human genes there were two or more distinct MPSS tags in the antisense orientation. Another interesting observation is the larger number of MPSS tags antisense to exonic regions of the sense genes. Unexpectedly, we found a much smaller number of antisense tags for mouse (Table 4). Although the number of mouse libraries is larger (81 mouse and 41 human libraries), the number of unique tags is significantly smaller (56,061 for mouse and 340,820 for human). The assignment of these unique tags to known genes shows a smaller representation of known genes in the mouse dataset (51% against 66% for human). It is unlikely, however, that these differences can explain the dramatic difference shown in Table 4. Further analyses are needed to solve this apparent discrepancy.
To experimentally validate the existence of these novel human NAT candidates we used the GLGI (Generation of Longer cDNA fragments from SAGE for Gene Identification)-MPSS technique [38] to convert 96 antisense MPSS tags into their corresponding 3' cDNA fragments. A sense primer corresponding to the antisense MPSS tag was used for GLGI-MPSS amplification as described in Materials and methods. A predominant band was obtained for most of the GLGI-MPSS reactions ( Figure 1). Amplified fragments were purified, cloned, sequenced and aligned to the human genome sequence. We were able to generate a specific 3' cDNA fragment for 46 (50.5%) out of 91 novel antisense candidates. Of these 46, the poly-A tail of 19 aligned with stretches of As in the human genome sequence (this finding will be discussed further). The existence of three of these antisense transcripts, out of three that were tested, was further confirmed by orientation-specific RT-PCR (data not shown).
Among the 49.5% (91 -46 = 45) of candidates that were not considered to be validated, we found 25 that were amplified in the GLGI-MPSS experiment but whose exon-intron organ- Exonic and intronic refer to the genome organization of the sense gene. For instance, the category 'One exonic tag' corresponds to those genes with only one antisense tag complementary to its exonic region. All identified tags are found at a frequency ≥3 tags per million (see Materials and methods).
ization was identical to the sense gene. Although antisense sequences like these have already been observed [39], we did not consider them as validated antisense transcripts. Orientation-specific RT-PCR confirmed the existence of one transcript, out of two that were tested.

Alternative polyadenylation as a major factor in defining S-AS pairs
Dahary et al. [6] observed that S-AS overlap usually involves transcripts generated by alternative polyadenylation. This observation had already been reported by us and others [40]. We decided to test if these preliminary observations would survive a more quantitative analysis. We found that the S-AS overlap is predominantly due to alternative polyadenylation variants. Roughly, 51% of all S-AS pairs (274 out of 533 3'3' pairs) overlap due to the existence of at least one variant. This number is certainly underestimated since many variants are still not represented in the sequence databases. The above observation raises the exciting possibility that antisense regulation is associated with the regulation of alternative polyadenylation. It is expected that the presence of overlapping genes imposes constraints on their evolution since any mutation will be evaluated by natural selection according to its effect in both genes. Thus, in principle, overlapping genes should impose a negative effect on the fitness of a subject. Alternative polyadenylation has the potential to relax such negative selection since the overlapping is dependent on a post-transcriptional modification.
If alternative polyadenylation is a significant factor in defining S-AS pairs, we would expect a lower rate of alternative polyadenylation in chromosome X, which has the smallest density of S-AS pairs. Indeed, only 20% of all messages from the X chromosome show at least two polyadenylation variants, compared to 27.5%, on average, for the autosomes (chisquare = 34.91, df = 1, p < 0.0001).

A fraction of S-AS pairs is generated through internal priming and retroposition events
During the validation of new NATs identified using the MPSS data, we noticed that a significant fraction of GLGI amplicons (19 out of 46 validated fragments) had their 3' ends aligning GLGI-MPSS amplification Figure 1 GLGI-MPSS amplification. GLGI amplifications for 96 MPSS antisense tags were analyzed on agarose gels stained with ethidium bromide. Note that some lanes show only a single amplified band whereas others have more than one band and sometimes a smear. A 100 bp ladder (M) was used as molecular weight marker.
M 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29  to stretches of As in the human genome. This motivated us to search for similar cases in the set of cDNA-based S-AS pairs identified in this study. We found that 18% and 26% of all S-AS pairs have at least one gene with its 3' end aligning with a stretch of A's in the human and mouse genomes, respectively. This number is certainly inflated by ESTs since it decreases to 11.7% for human and 12.6% for mouse when only mRNA/ mRNA S-AS pairs are considered. Two possibilities could RT-PCR analysis for the internal priming (IP) candidates in fetal liver, colon and lung total RNA Figure 2 RT-PCR analysis for the internal priming (IP) candidates in fetal liver, colon and lung total RNA. RT-PCR was conducted in DNA-free RNA previously treated with DNAse (lanes 1 and 2) and in untreated RNA, which was, therefore, contaminated with genomic DNA (gDNA; lanes 3 and 4) for each candidate in the corresponding tissue. As a control, RT-PCR was conducted in the presence (lanes 1 and 3) and absence (lanes 2 and 4) of reverse transcriptase. gDNA was used as a positive control of the PCR reaction (lane 5) and no template as a negative control (lane 6). For fetal liver, in 3 IP candidates (5, 8 and 11) the PCR products (152 bp, 153 bp and 160 bp, respectively) were observed in the treated RNA when RT was added (lane 1) or in untreated RNA independent of the RT (lanes 3 and 4). For colon, in 1 IP candidate (9) the PCR product (158 bp) was observed in the treated RNA when RT was added (lane 1) or in untreated RNA independent of the RT (lanes 3 and 4). For the remaining IP candidates (1, 2, 4, 6, 7, 10 and 12), the PCR products ( account for this observation. First, a fraction of all antisense transcripts would be artifacts due to genomic priming with contaminant genomic DNA during cDNA library construction. An alternative is the possibility that antisense genes were constructed during evolution by retroposition events. Both possibilities are in agreement with the observation that antisense genes are depleted of introns. An experimental strategy was developed to evaluate the likelihood of genomic priming as a factor generating artifactual antisense cDNAs. A total of 11 mRNA candidates derived from cDNA libraries from fetal liver, colon and lung with a high proportion of sequences that had their 3' ends aligning to stretches of As in the human genome were selected for experimental validation by RT-PCR. cDNA samples used in these experiments were reverse transcribed from fetal liver, colon and lung total RNA treated or not with DNAse. As can be seen in Figure 2, specific amplifications could not be achieved for 7 (63.6%) out of the 11 selected candidates when cDNA sam-ples used as templates for PCR amplification were prepared from DNA-free RNA. On the other hand, when untreated RNA was used for cDNA synthesis, all candidates could be amplified, suggesting that a significant proportion of these internal priming sequences were indeed generated from contaminant genomic DNA.
Some other features support the artifactual origin of these antisense transcripts. First, cDNAs containing a stretch of As at their 3' genomic end have much less polyadenylation signals than genes in general (17% compared to 85%). Furthermore, these genes have a much narrower and rarer expression pattern when analyzed by SAGE and MPSS than genes in general (data not shown). These observations suggest that a significant fraction of all antisense genes are actually artifacts, due to genomic priming during library construction.
Retroposition generates intronless copies of existing genes through reverse transcription of mature mRNAs followed by integration of the resulting cDNA into the genome (for a review, see Long et al. [41]). Eventually, the cDNA copy can be involved in homologous recombination with the original source gene as has been suggested for yeast [42]. Retroposition was thought to generate non-functional copies of functional genes. However, several groups have shown that retroposition has generated a significant amount of new functional genes in several species [43][44][45]. Recently, Marques et al. [43] found almost 4,000 retrocopies of functional genes in the human genome. More recently, the same group reported that more than 1,000 of these retrocopies are transcribed, of which at least 120 have evolved as bona fide genes [46].
Retrocopies usually have a poly-A tail at their 3' end because of the insertion of this post-transcriptional modification together with the remaining cDNA. Thus, retroposition can explain the high incidence of antisense transcripts with a poly-A tail at their 3' end. To evaluate the contribution of retrocopies to the formation of S-AS pairs we compared the loci identified by Marques et al. [43] as retrocopies with the list of S-AS pairs identified in this study. Out of 413 retrocopies represented in the cDNA databases, 138 were involved in S-AS pairs (70 mRNA/mRNA and 68 mRNA/EST pairs). For the 70 mRNA/mRNA pairs, 78% were classified as embedded. This is in agreement with our previous observation that embedded pairs are enriched with intronless genes. Thus, retroposition seems to significantly contribute to the origin of embedded S-AS pairs.

Expression patterns within S-AS pairs
A critical issue to effectively evaluate the role of antisense transcripts in regulating distinct cellular phenomena is related to the expression pattern of both sense and antisense transcripts belonging to the same S-AS pair. Several reports have been published based on large-scale gene-expression analyses [5,19,23,47,48]. Similar to Wang et al. [48], we here used MPSS libraries available for human to explore this issue.
Expression pattern (in a set of 31 tissues covered by MPSS) of genes belonging to all three types of S-AS pairs (3'3', 5'5'and embedded) Figure 3 Expression pattern (in a set of 31 tissues covered by MPSS) of genes belonging to all three types of S-AS pairs (3'3', 5'5'and embedded). (a) Categories are as follows: 'no expression', for S-AS pairs whose expression was not detected (see Materials and methods for details); 'single-gene expression', for S-AS pairs in which expression is observed for only one gene in the pair; 'co-expression', for pairs in which expression is seen for both genes in the pair. (b) Rate of differential expression for the set of co-expressed S-AS pairs. Ratio of sense/antisense genes in the pair is shown on the x-axis. Tag to gene assignment was performed as previously described [31,49]. To ensure the MPSS sequences were unambiguously matched to the assigned transcript, we removed tags mapped to more than one locus. Frequencies for all tags assigned to genes in an S-AS pair were collected from all MPSS libraries. Figure 3 shows the expression pattern of S-AS pairs for all MPSS libraries for human. We divided the dataset into the following categories as before: 3'3', 5'5' or embedded. Several features are evident. The rate of co-expression in our dataset was 35.1% compared to 44.9% observed by Chen et al. [4]. The differences are probably due to experiment design in both reports (for example, differences in the dataset and in the way the rate was calculated). Second, the rate of coexpression is significantly higher for 3'3' pairs when compared to the frequency of the embedded pairs (50.3%, chisquare = 134, df = 1, p = 5.4 × 10 -31 ). This supports a previous conclusion from Sun et al. [5] that 3'3' S-AS pairs are significantly more co-expressed than other pairs and, therefore, are more prone to be involved in antisense regulation. It is important to mention that 5'5' pairs are also enriched in coexpressed pairs when compared to embedded pairs (chisquare = 23.5, df = 1, p = 1.2 × 10 -6 ). We observed no statistical difference among the three categories regarding differential expression of both genes in a pair.

Influence of antisense transcripts in the splicing of sense transcripts
It is quite clear nowadays that a significant fraction of all human genes undergo regulated alternative splicing, producing more than one mature mRNA from a gene (Galante et al. [27] and references therein). Although several regulatory elements in cis and trans have been identified (for a review see Pagani and Baralle [50]), it is reasonable to say that we are far from a complete understanding of how constitutive and alter-native splicing are regulated. One possible regulatory mechanism involves antisense sequences. Since the late 1980s, it is known that antisense RNA can inhibit splicing of a pre-mRNA in vitro [15]. A few years later, Munroe and Lazar [51] observed that NATs could inhibit the splicing of a message derived from the other DNA strand, more specifically the ErbAα gene. More recently, Yan et al. [52] characterized a new human gene, called SAF, which is transcribed from the opposite strand of the FAS gene. Over-expression of SAF altered the splicing pattern of FAS in a regulated way, suggesting that SAF controls the splicing of FAS. With the growing amount of genomic loci presenting both sense and antisense transcripts, a general role for S-AS pairing in splicing regulation has been proposed [47]. However, no systematic large-scale analysis has been reported so far investigating this issue for mammals. We made use of the human dataset described in this report to tackle this problem.
We first tested whether the rate of alternative splicing in the sense gene would be affected by the existence of an antisense transcript. It is expected that the effect of S-AS pairing on splicing would be restricted to those exon-intron borders located in the region involved in pairing. We therefore restricted the analysis to those exon-intron borders spanning the region involved in an S-AS pairing. Our strategy was to compare the number of splicing variants for those borders against all other exon-intron borders (those without an antisense transcript) in the same genes. To make the analysis more informative we split the borders into four categories (terminal donor, internal donor, internal acceptor and terminal acceptor). For both internal donor and acceptor sites, the presence of an antisense transcript slightly increased the rate of alternative splicing (Table 5; 4% and 3% increases, respectively). For the terminal sites, the presence of a NAT had the opposite effect (5% and 6% decrease for donor and acceptor, respectively). Table 5 also shows that these differences are predominantly due to intron retention. On the other hand, NATs located within the introns and exons (but not spanning the border) have no major effect on the splicing of the respective borders. The observed differences between borders with or without NATs is statistically significant (chi-square = 31.2, df = 1, p = 2.3 × 10 -8 for donor sites; and chi-square = 23, df = 1, p = 1.6 × 10 -6 for acceptor sites).
Recently, Wiemann et al. [53] reported a new variant of IL4L1 that contains the first two exons of an upstream gene, NUP62. This chimeric transcript was expressed in a tissue and cellspecific manner. The authors speculated that cell type specific alternative splicing was involved in the generation of this chimeric transcript. We speculate that NATs could be involved in the generation of this type of chimeric cDNA. The same antisense message pairing with both sense messages would form a double-stranded RNA that could induce the spliceosome to skip the paired region and join the two sense messages, a process very similar to the one proposed for trans-splicing in mammals [54]. Interestingly, we found five examples in our dataset of S-AS pairs in which the genomic organization of both sense and antisense genes suggest a process like this. Additional data file 9 illustrates one of these cases. It can be seen that two transcripts represented by cDNAs AK095876 and AK000438 join messages from genes SERF2 and HYPK. The antisense transcript is represented by cDNA AK097682. Additional data file 10 lists all other putative cases of chimeric transcripts. The fact that both sense genes share a common antisense transcript raises the possibility that antisense transcripts can mediate trans-splicing of the sense genes, thereby generating the chimeric transcript.

On the evolution of S-AS pairs: functional implications
It is reasonable to assume that a fraction of all S-AS pairs reached this genome organization solely by chance. However, evidence presented here and elsewhere suggest that this fraction is probably small [6,55,56]. For example, Dahary et al. [6] concluded that antisense transcription had a significant effect on vertebrate genome evolution since the genomic organization of S-AS pairs is much more conserved than the organization of genes in general. However, how did this organization come to be? In principle, S-AS genomic organization should carry a negative effect on the overall fitness of a subject. For each gene in an S-AS pair, its evolution is constrained not only by features of its own sequence but also by functional features encoded by the other gene in the pair. The fact that we observed a significant amount of S-AS pairs in mammalian genomes suggests that there are advantages inherent to this organization to counter-balance the negative effects. The proposed role of NATs in gene regulation is certainly advantageous. We propose here two evolutionary scenarios, not mutually exclusive, that would speed up the generation of S-AS pairs. In one scenario, alternative polyadenylation has a fundamental role. Sun et al. [5] observed a preferential targeting of 3' UTRs for NATs. Our observation that 51% of 3'3' S-AS pairs overlap because of polyadenylation variants suggests that selection has favored cases where overlapping occurs only in a time and spatially regulated manner.
In a second scenario, retroposition generates NATs, which lack introns and may even show a polyadenylation tail integrated into the genome. We observe here that retroposition contributed significantly to the origin of S-AS pairs, especially those classified as embedded. What would be the selective advantages of retrocopies as NATs? Chen et al. [56] observed that antisense genes have shorter introns when compared to genes in general. They speculated that this feature was advantageous during evolution since NATs need to be "rapid responsers" to execute their regulatory activities. Although transcription is a slow process in eukaryotes, another bottleneck in the expression of a gene is splicing. Furthermore, Nott et al. [57] observed that the presence of introns in a gene affects gene expression by enhancing mRNA accumulation. Thus, the argument from Chen et al. [56] gets stronger with the data reported here and by Nott et al. [57] since intronless antisense genes would be transcribed even faster; their transcripts would simply skip splicing and the half-life of the respective messages would be shorter. All key features for genes involved in regulatory activities.
An important issue is the conservation of S-AS pairs between human and mouse. Although we found more than a thousand conserved pairs, this number is still small compared to the whole set of S-AS pairs in both species. Several factors, however, suggest that the number reported here is an underestimate. First, as discussed by Engstrom et al. [25], sequence conservation might not be of primary importance for antisense regulation. Furthermore, it is likely that many truly conserved pairs were not detected because transcript sequences have not been discovered yet. This is more critical in the face of our findings that a significant proportion of 3'3' S-AS pairs depend on alternative polyadenylation for an overlap. It is also quite likely that some S-AS pairs are lineage-specific. For instance, our finding that retroposition contributes to the origin of many S-AS pairs could explain the appearance of lineage-specific S-AS pairs, assuming that the retroposition event occurred after the divergence between human and mouse.
These two evolutionary scenarios (alternative polyadenylation and retroposition) might produce S-AS pairs with different functional implications. The expression and evolutionary conservation analyses presented here, together with evidence from others [5,19,23,47,48] suggest that 3'3' overlap achieved by polyadenylation variants was used throughout evolution to regulate gene expression. Those pairs generated through retroposition may be involved in some other types of regulation, such as alternative splicing.

Conclusion
This is the deepest survey so far of S-AS pairs in the human and mouse genomes. We made use of all cDNAs available in the public domain together with 122 MPSS libraries for human and mouse. The major findings of the present report include: as many as 10,077 and 8,091 S-AS pairs were identified for human and mouse respectively; using MPSS data, we found 4,308 and 216 new putative S-AS loci in human and mouse, respectively; a small fraction of all S-AS pairs are artifacts caused by genomic priming during cDNA library construction; a significant amount of S-AS pairs is due to retroposition events of one of the genes in the pair; quantitative analyses suggest that the presence of an antisense gene, complementary to an exon-intron border of the sense gene, increases the rate of retention of the respective intron. Furthermore, we propose an evolutionary model in which alternative polyadenylation and retroposition are important forces in the generation of S-AS pairs.
Taken together, these results offer, up to now, the vastest catalog of S-AS pairs in human and mouse genomes.

Mapping cDNAs and MPSS tags onto the human and mouse genomes
We used a modified protocol similar to the one described previously to identify transcription clusters in the human and mouse genomes [27,28]. Briefly, genome sequence (NCBI build no. 35 for human and NCBI build no. 33 for mouse), EST collections (5,992,459 sequences for human and 4,246,824 sequences for mouse) and mRNA sequences (186,358 for human and 120,058 for mouse) were downloaded from UCSC [58]. All cDNAs were mapped to the respective genome sequence using BLAT (default parameters) [59]. The best hit for each cDNA in the genome was identified, followed by a pairwise alignment using Sim4 [60]. Only transcripts presenting identity ≥94%, coverage ≥50% and all splice sites in the same orientations were used.
Correct orientation of ESTs was determined by the presence of a poly-A tail (a stretch of 8 As at the 3' end) and/or a splicing donor (GT) and acceptor (AG) sites. All mRNAs were considered in the 'sense' orientation (oriented from 5' end to 3' end). All cDNAs mapped and reliably orientated were assembled into clusters. One cluster contains cDNAs presenting the same orientation and sharing at least one exon-intron boundary or a minimum of 30 nucleotides of overlap (only for those sequences without a common exon/intron organization).
For the mapping of MPSS data, we first extracted 'virtual' tags for both human and mouse genomes by simply finding all DpnII sites and extracting a 13 (human) or 16 (mouse) nucleotide long sequence immediately downstream of the restriction site in both orientations. These 'virtual' tags present only once in the respective genomes were further used and matched against the 'real' tags found in 41 and 81 MPSS libraries for human and mouse, respectively. Only MPSS tags classified as 'reliable' (present in more than one sequencing run) and 'significant' (tags per million >3) were considered as trusted signatures.

Identification of S-AS pairs
S-AS pairs were identified as those cases in which two clusters, in opposite orientations, overlap at the genome level. For the correct orientation of all mapped cDNAs, we took into consideration several parameters, including: sequence annotation as available in the respective GenBank entry; splice junctions; and poly-A tails and poly-T heads. We excluded from our analyses all cDNAs that presented conflicting orientations as defined by the three criteria above. If only two clusters overlap in the opposite orientation, they were classified as a single bidirectional S-AS pair. If a given cluster overlaps with more than one antisense cluster, they were classified as multiple bidirectional S-AS pairs. S-AS pairs were also classified according to their genomic pattern. Parameters evaluated included: pattern of S-AS overlap (exonic, intronic and exonic/intronic); spanning of introns by the components of a pair as defined by their alignment onto the genome; and chromosome localization and relative orientation within the S-AS pairs (tail-tail, head-head and embedded).

Conservation between human and mouse S-AS pairs
We used three strategies to evaluate the degree of conservation between human and mouse S-AS pairs. First, all pairs were searched against the dataset from HomoloGene [35] and those pairs conserved in both species were counted. In our second strategy, we selected those S-AS pairs in which at least one gene was conserved according to HomoloGene. We then used Needle, an alignment algorithm [61], to test sequence conservation between the respective antisense genes. We classified as conserved those global alignments with identity >30%. Finally, we also used the strategy from Engstrom et al. [25]. We used the net alignment between human and mouse genomes (retrieved from the UCSC Genome Browser database) to define the corresponding (synthenic) regions. We considered a human S-AS pair to be conserved in mouse if it had an exon region aligning (>20 bp) to an exon region from a mouse pair.

Investigation of the expression pattern of S-AS transcripts
We evaluated the expression pattern of S-AS pairs at the whole genome level based on their expression profiles obtained from MPSS libraries (available at [36]). The procedure was previously described by us for SAGE and MPSS [27,31,49]. The tag to gene assignment was done by scanning and extracting virtual tags (13 nucleotide-long sequences present downstream to the 3'-most DpnII restriction sites of each mRNA sequence). To accurately represent the 3' end of a transcript, only mRNA sequences containing a poly-A tail were used. All tags mapped to two or more different genes Galante et al. R40.13  script. PCR amplifications were performed using 1 μl of the first-strand cDNA as a template in a final volume of 25 μl and 1× buffer, 1.5 mM MgCl 2 , 0.1 mM dNTP, 0.4 μM of gene specific primers and 1 U of Platinum Taq DNA polymerase (Invitrogen). The following cycling conditions were used for amplification: initial denaturation of 95°C for 2 minutes; 94°C for 40 s; reaction-specific annealing temperature for 40 s and 72°C for 1 minute for 35 cycles; followed by a final extension step at 72°C for 7 minutes. All PCR products were resolved on 8% polyacrylamide gels. Controls for the absence of self-priming during cDNA synthesis were done with reverse transcriptase in the absence of primers, and controls for the absence of DNA were done by incubation with primers but with no reverse transcriptase.

Availability
To make our dataset fully accessible to the community we have set up a worldwide web portal [62] containing all raw data generated in this study and a series of tools to explore the data.

Additional data files
The following additional data are available with the online version of this paper. Additional data file 1 is a list of representative GenBank entries for all S-AS pairs in both human and mouse. Additional data file 2 is a table showing the total number of S-AS pairs by chromosome for both human and mouse. Additional data file 3 shows the number of clusters and S-AS pairs when a less stringent clustering methodology is applied. Additional data file 4 shows a schematic view of all possible genomic organizations of S-AS pairs. Additional data file 5 lists all S-AS pairs conserved between human and mouse using the three strategies described in the text. Additional data file 6 shows the fraction of S-AS pairs conserved between human and mouse that are classified as 'Fully intronic' and the fraction of conserved S-AS pairs that contain at least one intronless gene. Additional data file 7 is a list of all MPSS libraries used in this study. Additional data file 8 presents the number of cDNA-based pairs that were further confirmed by the MPSS data. Additional data file 9 is a figure illustrating chimeric transcripts joining two adjacent genes (SERF2 and HYPK) with a NAT located between them. Additional file 10 lists all cases of chimeric transcripts identified in our dataset.
Additional data file 1 Representative GenBank entries for all S-AS pairs in both human and mouse Representative GenBank entries for all S-AS pairs in both human and mouse. Click here for file Additional data file 2 Total number of S-AS pairs by chromosome for both human and mouse Total number of S-AS pairs by chromosome for both human and mouse. Click here for file Additional data file 3 Number of clusters and S-AS pairs when a less stringent clustering methodology is applied Number of clusters and S-AS pairs when a less stringent clustering methodology is applied. Click here for file Additional data file 4 All possible genomic organizations of S-AS pairs All possible genomic organizations of S-AS pairs. Click here for file Additional data file 5 All S-AS pairs conserved between human and mouse using the three strategies described in the text All S-AS pairs conserved between human and mouse using the three strategies described in the text. Click here for file Additional data file 6 The fraction of S-AS pairs conserved between human and mouse that are classified as 'Fully intronic' and the fraction of conserved S-AS pairs that contain at least one intronless gene The fraction of S-AS pairs conserved between human and mouse that are classified as 'Fully intronic' and the fraction of conserved S-AS pairs that contain at least one intronless gene. Click here for file Additional data file 7 All MPSS libraries used in this study All MPSS libraries used in this study. Click here for file Additional data file 8 The number of cDNA-based pairs that were further confirmed by the MPSS data The number of cDNA-based pairs that were further confirmed by the MPSS data. Click here for file Additional data file 9 Chimeric transcripts joining two adjacent genes (SERF2 and HYPK) with a NAT located between them Chimeric transcripts joining two adjacent genes (SERF2 and HYPK) with a NAT located between them. Click here for file Additional data file 10 All cases of chimeric transcripts identified in our dataset All cases of chimeric transcripts identified in our dataset. Click here for file