EST Alignments andClustering

Aligning ESTs by inter se sequence comparison allows the grouping of clones derived from the same transcript and is the most common method used to estimate the number of unique transcripts represented in an EST collection. We used this approach in our first EST project to produce the DGCr1; Rubin et al. 2000). 5′ ESTs were grouped using BLAST (Altschul et al. 1997) and then assembled using Phrap (Ewing and Green 1998; Ewing et al. 1998). The 3′-end sequence was then determined from the clone that extended the most 5′ in each cluster. Comparison of the 3′ ESTs was used as a final check for redundancy because 5′ ESTs from incomplete cDNA clones often do not overlap with ESTs from full-length clones derived from the same transcript.

Sim4 is a software tool designed for aligning ESTs or cDNAs with genomic sequence allowing for gaps at the positions of introns. Using Sim4 to scan the entire genome for sequences matching each EST proved impractical because of the long compute time required. We first used BLAST to align the ESTs with 100-kb genomic segments and then further characterized the transcript structure of each EST using Sim4. We successfully aligned 213,238 ESTs to the genome using a high stringency cutoff that required >90% of the EST sequence to be aligned.

We used an iterative approach to cluster ESTs that was based on their alignments with the genomic sequence. If any exons in an EST alignment overlapped exons within an existing cluster, the EST was incorporated into that cluster and the intron–exon structure(s) of the cluster updated. Because one gene could nest in the intron of another gene, EST alignments were grouped on the basis of overlap between exons and not introns. Figure ​Figure11 shows a user interface we developed to view the structures of EST clusters. We subdivided EST clusters such that each subcluster contained only ESTs whose alignments had the same intron–exon structure. For 3′ ESTs, we placed that EST in the same subcluster as the 5′ EST from the same cDNA clone. For the EST cluster shown in Figure ​Figure1A,1A, there are six subclusters (numbered SC2 to SC7) that have been color coded to indicate the number of their constituent ESTs. These subclusters were further combined into merged subclusters that most likely reflect the structure(s) of alternatively spliced transcripts (Fig. ​(Fig.1B).1B). However, the data are insufficient to establish whether all the ESTs that make up a merged subcluster are in fact derived from a single splice form. Thus the gene models indicated by the merged subclusters, whereas consistent with the available data, are not proven and do not represent all possible models. Because our initial goal was to identify a single full-length cDNA clone for each gene, we used a very conservative process in selecting clones for the DGCr2, which is described in detail below. We selected only a single clone from each cluster for inclusion, irrespective of the number of subclusters. In the case illustrated in Figure ​Figure1B,1B, we selected the clone corresponding to the EST extending the farthest 5′ in subcluster 2 (SC2). Our analysis identified cDNAs derived from apparent alternatively spliced transcripts for >20% of the genes with ESTs. These data and clones will be used in the future to isolate and further characterize cDNAs representing alternative splice forms.

Open in a separate window

Figure 1

Graphical display of expressed sequence tag (EST) clusters. (A) is an example of six different subclusters aligning to the Curated Gene CG8977 and numbered SC2-SC7. The subclusters are color coded with respect to the number of EST members as shown. (B) is a gene model based on the merged subclusters illustrating two possible splice variants and numbered SC2-as and SC4-as. SC2-as is a merge of SC2, -3, -5, and -7 and SC4-as is a merge of SC4 and SC6.

This clustering approach identified 16,744 clusters from the 213,238 successful EST alignments. Annotations of the Release 2.0 Drosophila genome sequence are composed of gene predictions and previously characterized genes, designated as Curated Genes (CGs). Of the 16,744 clusters, 12,154 overlap 9664 different CGs. The remaining 4590 clusters do not directly overlap a CG. We expect these nonoverlapping clusters for two reasons. First, many 5′ ESTs are composed entirely of 5′-UTR sequences, whereas the majority of CGs result from gene-prediction algorithms that only predict open reading frames (ORFs) and thus do not include the sequences of the UTRs. For example, 13% of full-length sequences from our gene collection have a 5′ UTR >500 bp in length. The average sequence trace from the ABI-PRISM DNA sequencer (Applied Biosystems) is ~500 bp in length after quality trimming, so we would expect that some of these clusters would not overlap with a gene prediction. Second, 36% of the nonoverlapping clusters are >20 kb upstream of a CG and likely represent new genes that were not detected by gene-finding programs. Verifying that these nonoverlapping clusters represent authentic genes will require additional sequence and gene-expression information.

Genomic contamination or reverse transcription of unspliced mRNA precursors has the potential to corrupt gene collections. To evaluate our EST clusters for such artifacts, we examined separately the position, relative to CGs, of EST clusters whose alignments indicate that they contain spliced ESTs and those containing only unspliced ESTs. Of the 12,154 clusters that overlap a CG, 72% contain spliced ESTs. Those 1455 clusters that do not overlap a CG but lie within 5 kb upstream of a CG contained spliced ESTs 28% of the time, whereas only 18% of the 3135 clusters that lie >5 kb upstream of a CG contained spliced ESTs. To minimize the number of cDNAs derived from unspliced transcripts that were selected, we made the conservative decision to select ESTs only from clusters that contained spliced ESTs, unless the cluster directly overlapped a CG, in which case, we also accepted clusters comprised only of unspliced ESTs.

Selecting DGC2

In the absence of direct overlap between an EST cluster and the closest downstream CG, it proved difficult to determine unambiguously if they represent the same gene. We developed the following rule-based approach to select clones for the DGCr2 in such cases. We defined a “CG region” as the genomic region extending from the 3′ end of a CG to the 3′ end of the adjacent upstream CG (Fig. ​(Fig.2A).2A). In cases in which the intergenic distance was >5 kb, we arbitrarily limited the CG region to the span extending from the 3′ end of the CG to 5 kb upstream of its 5′ end and considered the remainder of the intergenic space as a distinct CG region (Fig. ​(Fig.2B).2B). We made the further assumption that a CG region contains no more than one gene. The types of alignments observed for DGCr2 clusters and their frequencies are depicted in Figure ​Figure2,2, panels C through G. In 310 cases, clusters overlapped more than one CG as in Figure ​Figure2F.2F. Such cases were evidence that original gene predictions erroneously split a single gene into two CGs.

Open in a separate window

Figure 2

CG regions. The observed percentage of each alignment type for clusters representing DGCr2 is indicated in the right column. CG regions are represented as open boxes. CGs are shown numbered in gray boxes (exons) connected by lines (introns). The genomic sequence is shown as a hatched line. Aligned EST clusters are shown in black, and the ESTs chosen for DGC2 are gray with a black border. All CGs and alignments are shown representing one strand of the genome proceeding 5′ to 3′ as indicated in (A). See text for a full description of CG regions.

After defining the CG regions, the rules for choosing clusters and the clone in each cluster for inclusion in DGCr2 are straightforward. CG regions that already had a DGCr1 clone aligned within them were omitted from the analysis. For the CG regions that had no DGCr1 clone aligned, we examined each cluster within the region progressing from 5′ to 3′. For regions that contained both a cluster directly overlapping the CG and a nonoverlapping cluster more 5′, we chose the most 5′ cluster if it contained a spliced EST. In cases in which the more 5′ cluster contained only unspliced ESTs, we selected the EST from the cluster that overlapped the CG. After a cluster was identified, we selected the EST alignment that extended the most 5′ within the cluster. Using these criteria, we identified 5061 new clones for the DGC. Preferentially picking the DGCr2 clones from clusters that directly overlap CGs might result in our choosing a significant percentage of clones that were not full length; however, this does not appear to be the case. We compiled a test set of 328 CGs represented by a DGCr2 clone and for which full-length cDNA sequences previously existed in GenBank. We asked whether the DGCr2 clone extended as far 5′ as the GenBank sequence and whether the DGCr2 clone contains the complete ORF. For the 53 cDNAs in our test set <1 kb in length, the corresponding DGC cDNA was as long or longer than the test set cDNA 77% of the time and contained the complete ORF 100% of the time. For the 101 cDNAs in the test set between 1 kb and 2 kb, these numbers were 71% and 99%; for those 74 cDNAs between 2 and 3 kb, these numbers were 66% and 88%; and for those 100 cDNAs >3 kb, these numbers were 59% and 91%. This analysis shows that larger genes are more likely to be represented as an incomplete cDNA. However, we conclude that our selection criteria did not significantly bias the DGC toward incomplete clones.

We were also concerned that the high initial gene discovery rate of the AT library, which had the consequence that 1000 (20%) of the new clones chosen for the DGCr2 are AT clones, was caused by testis-specific transcription start sites of genes also expressed in other tissues. However, we found that 768 of the 1000 AT clones in DGCr2 represented clusters consisting of only AT ESTs. Furthermore, 663 (86%) of these AT-only clusters did not contain ESTs from other tissues in the same CG region. Moreover, preliminary evidence from GeneChip Array (Affymetrix) experiments comparing RNA from male and female adults indicates that transcripts corresponding to 30% of these AT-only clusters are present at >25-fold higher levels in adult males (P. Spellman and G.M. Rubin, pers. comm.). Although a more rigorous analysis of transcript structures and expression patterns will be necessary to address this question more directly, our initial analysis indicates that many of the testis clones in the DGC are abundantly expressed in testis and may be testis-specific genes.

Table ​Table22 shows the number of clones from each library included in DGCr1 and DGCr2. Clones from the AT, RE, and RH libraries account for 76% of DGCr2, indicating the value of these new libraries in identifying clones for the DGC. Sequencing more clones from older libraries as well as realigning all ESTs from the original project also proved valuable. A total of 1205 (24%) DGCr2 clones are from the original EST collection; of those 1205 clones, 328 (27%) are new ESTs from the SD and GM libraries, and the other 882 (73%) are ESTs from the original sequencing effort. Thus, 17% of the new DGCr2 clones already existed but were not included in DGCr1. Many of these clones represent replacements for clones that were originally selected for the DGCr1 but which subsequently failed one of the quality control tests used to validate this collection (Rubin et al. 2000). Moreover, in selecting clones for the DGCr1 we generally did not include CGs represented by a single EST. Genome alignments allowed us to identify which of these singleton clusters contained a spliced EST and gave us confidence that the EST represented a legitimate transcript. In total, DGCr2 represents 4274 ESTs that overlap a single CG, 113 that overlap two adjacent CGs, 198 that are within 5 kb upstream of a CG, and 476 that are spliced but align >5 kb upstream of a CG. If the 10,910 clones in our combined DGCr1 plus DGCr2 collections do indeed represent unique genes, the DGC would consist of ~81% of genes predicted to exist in D. melanogaster.

End Sequencing and Analysis

Random cDNA clones were end sequenced from the 5′ end using ABI Big Dye II Dye terminator chemistry. End sequencing of 3′ ends was performed using ABI d-rhodamine dye terminator chemistry on the AT, GM, and SD clones and ABI Big Dye II Dye terminator chemistry for the RE and RH clones. Primers for each of the cDNA cloning vectors are described in detail on the BDGP web page (http://www.fruitfly.org/about/methods/index.html). Reactions were run on ABI 3700 capillary sequence machines and chromatograms were evaluated using Phred (Ewing and Green 1998; Ewing et al. 1998). Reads were vector trimmed using Crossmatch (P. Green, http://bozeman.mbt.washing-ton.edu/phrap.docs/phrap.html). The ESTs were quality trimmed by evaluating chromatogram trace peaks using a Java program written by S. Lewis that implements an algorithm described previously (Hillier et al. 1996). The high quality cutoff was defined at the point in which either of the following conditions was true: The ratio between peak height and valley fell below a threshold or the ratio of the areas beneath the highest peak and the next highest fell below a threshold. This point was extended to the location where 4 consecutive bases with scores <q15 occurred. This marked the end of the “high quality” sequence. By further evaluating the Phred quality scores, the reads were then trimmed for submission at the base before a run of 6 bases with scores <q10. This marked the end of the submitted sequence. Reads were also screened for Escherichia coli transposons using the Genbank data set GB.vector. Reads that were >150 bp of contiguous high quality sequence were submitted to the National Center for Biotechnology Information (NCBI) and have accession numbers BF485518-BF503517, BF503521-BF506780, BG631888-BG631996, BG633696-BG637540, BG640063-BG641469, BI141709-BI142246, BI161485-BI173971, BI212109-BI216987, BI227448-BI233322, BI234009-BI243989, BI351612-BI354228, BI354231-BI355901, BI355935-BI358751, BI361285-BI376197, BI481532-BI487261, BI563331-BI593695, BI604243-BI620155, BI620158-BI635012, BI635064-BI638027, and BI638030-BI642053. Before high quality trimming, the average number of bases with phred quality ≥q20 for ESTs submitted to Genbank was 554. Clones corresponding to all ESTs, including those that make up DGCr2, are available from Research Genetics (http://www.resgen.com/products/DEST.php3). We are in the process of colony purifying and rearraying the DGCr2 clone set and expect it to be available for distribution by July of 2002.

We thank R. Hoskins, A. Huang, and S. Misra for critically reading and improving the manuscript. We thank H. Guarin for excellent technical assistance and the entire staff of the BDGP sequencing center. This work was funded by grants from the Department of Energy (grant no. DE-FG03-98 ER62625 and DE-FG03-99ER62739 and National Institutes of Health (grant no. P50 HG00750) to G.M. Rubin. This work was also supported by the Riken Genome Exploration Research grant to Yoshihide Hayashizaki.

The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked “advertisement” in accordance with 18 USC section 1734 solely to indicate this fact.