From First Base: The Sequence of the Tip of the X Chromosome of Drosophila melanogaster, a Comparison of Two Sequencing Strategies

Linking the Genetic Map of the X Chromosome to a Molecular Framework

A decade ago, the founding members of the EDGP argued the case for constructing an accurate physical map of the genome of D. melanogaster linked to the genetic map (Sidén-Kiamos et al. 1990). To this end, cosmid clones were selected by hybridization with PCR-amplified DNA microdissected from each of the 100 individual divisions of the major polytene chromosome arms. A physical map was generated by determining overlaps between the cosmids based on the shared fragments generated by restriction endonuclease digestion (Sulston et al. 1988). The localization of cosmids was verified by in situ hybridization to the polytene chromosomes and by determining STSs of cosmid end sequences (Louis et al. 1997). This physical map, and the cosmid library on which it was based, are available as a public resource (http://www.hgmp.mrc.ac.uk/Biology/descriptions/drosophila.html).

A physical map was also constructed by the BDGP (Kimmerly et al. 1996) based on segments of DNA cloned in a P1 phage vector that were aligned using PCR based STS content mapping. However, it was clear that both the cosmid and P1 maps would be an incomplete resource for sequencing the genome. Moreover, although the YAC map of Ajioka et al. (1991) does give good coverage, in our hands YAC clones were impractical for DNA sequencing purposes. We therefore undertook to build another map based on BAC clones because these vectors can, in principle, accommodate larger inserts of DNA. The generation of these BAC clones, that give an approximately 10-fold coverage of the genome, will be described in detail elsewhere. The library is available as a public resource (http://www.hgmp.mrc.ac.uk/Biology/descriptions/dros_bac.html). Clones from both this and a BAC library of partial EcoRI digestion products of DNA constructed for the BDGP (Hoskins et al. 2000) were physically ordered and linked by hybridization with a total of 647 hybridization probes each of 40 nucleotides in length corresponding to sequences distributed along the length of theX chromosome. The resulting maps, whose full description will also be provided elsewhere, allowed us to determine a minimal tiling path of clones for sequencing purposes. We selected such a minimal tiling path extending through polytene divisions 4–10, and determined the sequence of these clones at ∼1.5-fold coverage (http://edgp.ebi.ac.uk/cgi-bin/progress.pl). This provided a skeletal sequence scan of ∼6 Mb of the chromosome that was made available to the Celera/BDGP/Baylor shotgun sequencing project for use as an assembly scaffold.

The accurate sequencing of polytene divisions 1–3 was initiated on a minimal tiling path of cosmid clones, subsequently extended using the BAC clones to fill gaps in the cosmid map. The clones selected for sequencing are presented in Figure 1A, and the assembled nonredundant sequence can be directly accessed athttp://edgp.ebi.ac.uk/cgi-bin/progress.pl, which links to the EMBL-Bank deposits.

Figure 1.

Physical maps of the interval 1A–3C. (A) Minimal tiling pattern of clones sequenced in divisions 1A–3C. BACR clones are indicated in red; BACN and BACH clones are indicated in green; cosmid clones are indicated in blue; redundant clones sequenced are indicated in pink; a few small regions were sequenced from other clones, these are indicated in yellow. The BACR, BACN, and BACH clones are from the same strain as that sequenced by the BDGP and Celera; the cosmids are from a different strain (see Methods). Scale divisions are 10 Kb. (B) Genes, transposable elements, and P-elementinsertions in divisions 1A–3C. Known genes are shown in red; genes with significant protein similarities to nondrosophilid proteins are shown in blue; predicted genes with EST hits are shown in yellow; predicted genes with no EST hits are shown in green; predicted genes with protein motif matches are shown in pink. Transposable elements are shown in orange within the sequence coordinate line. The sites ofP-element and EP-element insertions are indicated by gray triangles. The large square brackets from 2100 to 2480 Kb embrace the zeste-white region (Figure 2). Scale divisions are 10 Kb (bold) and 1 Kb (regular).

Open in new tabDownload PowerPointLink to figure

General Features of Gene Content

As explained in Methods, we have used two general classes of computational method to predict genes in this chromosome region: similarity-based methods and ab initio methods. Together these two approaches have enabled us to predict 277 protein-coding genes overall, of which 94 (33.9%) had been sequenced previously by the community (Table1; Figure 1B). A total of 25 genes (9%) were predicted solely by ab initio methods, a lower fraction than in the Adh region (19%). A possible reason for this difference is that we used a stricter criterion for accepting a gene predicted only by an ab initio method than did Ashburner et al. (1999). Of the predicted genes, 205 have matches with ESTs from the BDGP (Rubin et al. 2000b) and NIH (Andrews et al. 2000) projects. The fraction of previously knownDrosophila genes that had EST matches (77.1%) is the same as that of the genes predicted by sequence similarity (77.2%), and is very similar to the proportion of matches from the Adh region (71%). Assuming that the criteria used to predict genes are adequate, these figures provide a good indication for the proportion ofDrosophila genes currently represented in EST collections. Presumably the shortfall reflects mainly that the cDNAs used to generate the ESTs have been derived from a restricted number of developmental stages. The value of ESTs in confirming gene identity and splicing patterns provides a strong argument to extend the generation of EST data to other developmental stages and tissues (Andrews et al. 2000; Rubin et al. 2000b). Based on the analysis of EST hits, we identified nine genes that are alternatively spliced in their coding regions, and thus able to direct the synthesis of two or more different proteins (Table 1, asterisks). It is striking that of the 183 newly predicted genes, 55% have significant similarities with sequences in other organisms thus indicating the extent of conserved function.

The average size of the coding regions of the genes predicted in the tip of the X chromosome is 1.8 Kb, with 2.7 introns per gene. The gene with the highest number of introns isEG:BACR25B3.1 (26 introns in the coding region). The average size of the introns is 475 bp, with the shortest being 26 bp (EG:63B12.3) and the longest being 34,401 bp (sidekick [sdk], EG:BACR19J1.1). The calculated average number of introns per gene in this chromosomal region is consistent with previous studies that have indicated the majority of Drosophila genes contain one or two small introns located near their 5′ ends (although exon and intron numbers will have been underestimated as ab initio gene prediction methods will not predict untranslated exons). There are, however, some exceptionally large genes. These include sdk, which encodes an immunoglobulin-C2 domain protein, and is required to prevent the“mystery cell” of the developing eye disc differentiating as a photoreceptor (Nguyen et al. 1997). This gene, sequenced previously as a cDNA, covers 60 Kb and includes at least 14 exons. Another very large gene is futsch (EG:49E4.1), covering 18 Kb and encoding a protein of 5327 amino acids predicted to encode a microtubule-associated protein, on the basis of its similarity with human MAP1B (SWISS-PROT:P46821), which is only half the size. RecentlyHummel et al. (2000) have shown that futsch encodes the well-known Drosophila neural antigen 22C10. Four other genes have large transcription units: Appl, 35.1 Kb; br, 27.7 Kb: EG25B3.1, 20.0 Kb; and csw, 17.4 Kb. The overall GC content of this collection of genes from the tip of theX chromosome is significantly lower (45.5%) than the overall GC content of the genes in the Joint Sequence (56.1%).

One of the surprising results of the analysis of the Adhregion sequence (Ashburner et al. 1999) was the number of genes predicted to be included within the introns of other genes (8%). These were most frequently, but not exclusively, arranged as anti-parallel transcription units. The present analysis of the tip of theX permits a comparison with another segment of genomic DNA. We predict four nested genes. This corresponds to 1.4 % of all of the genes we identify. This is probably an underestimate, because ab initio gene prediction programs do not predict genes within genes.

One group of duplicated genes worthy of specific mention in this region are the cytochrome P450s, small monooxygenases often involved in the metabolism of xenobiotic compounds. Eighty-seven genes encoding these microsomal or mitochondrial enzymes had been identified in the essentially complete Joint Sequence of D. melanogaster (Nelson 2000). Only two (l(2)35Fb in the Adh region [Ashburner et al. 1999] and disembodied [Chávez et al. 2000]) have been associated with a mutant phenotype, although polymorphisms at others have implicated them in differential resistance to DDT and other compounds (Berge et al. 1998). One characteristic of the genes encoding these proteins is that they often occur in small clusters, indicating an expansion of the gene family by duplication. In region 1–3 we have identified five cytochrome P450-encoding genes (Cyp4g1, Cyp4d1, Cyp4d2, Cyp4ae1, and Cyp4d14); of these, the latter three are in tandem within about 7.5 Kb at 2E1 andCyp4d1 is some 12 Kb distal at 2D6. The Cyp4g1 (at 1B4) gene appears to be more abundantly transcribed than any other P450 gene in D. melanogaster, at least judging from the large number of its EST sequences (59; Nelson 2000).

We have analyzed all of the known or predicted proteins by several methods, most extensively by BLASTP against data sets derived from SWISS-PROT and TrEMBL sorted by taxonomic origin (seeAshburner et al. 1999). We have also analyzed all of the protein sequences by various methods to detect protein motifs, and domains. Overall, 71% of the known or predicted proteins have aBLASTP match with an expectation of 10⁻⁷ or less when compared with nondrosophilid protein sequences. Similarly, 137 contain at least one known motif or domain (other than the PROSITE Nuclear Localization Signal profile) as determined by matches against InterPro (http://www.ebi.ac.uk/interpro/). These numbers are, of course, both preliminary and transitory. All of these data have been communicated to FlyBase and can be found in the supplementary data (see Methods). We have chosen only to present the PFAM hits in Table 1, as an indication of the data obtained.

As we have discussed previously (Benos et al. 2000), examples of 12 different transposable elements were identified within the region analyzed: 412, roo, Doc, FB, jockey, mgd1, Tirant, S-element, 1360, Burdock, blastopia, and yoyo. It is possible that more transposable elements may be present in the region; however, we have not identified them molecularly.

Chromosomal Regions of Particular Interest

The zeste-white Region

The discovery of polytene chromosomes in the larvae ofDrosophila in the early 1930's was a major event in the history of genetics. These chromosomes are characterized by a nonperiodic pattern of darkly staining bands and lightly staining interbands, reflecting differences in the degree of DNA packing. These patterns are both colinear with the genetic map, as proven by Bridges (1937) and extraordinarily stable; they can be recognized in species that have diverged many millions of years ago. The detailed maps of Bridges (see Lefevre 1976; Sorsa 1988) enumerated 5072 polytene chromosome bands (and, hence, interbands). Bridges suggested, somewhat tentatively, that there may be a one-to-one correspondence between these bands and genes, a hypothesis that became known as the “one band/one gene hypothesis”. A prediction of this hypothesis was that Drosophila had ∼5000 genes. This idea was apparently supported by estimates of the number of vital loci on the X chromosome, ∼1000 or ∼5000 for the genome as a whole (Lea 1955; Lefevre and Watkins 1986). Further apparent confirmation of the one band/one gene hypothesis came from a number of attempts to “saturate” small regions of the genome with mutations, and hence estimate the number of genes in that region (Alikhanian 1937). Most famous of these experiments was that of Judd and students (Judd et al. 1972; Young and Judd 1978) who studied a small region of the distal X chromosome between bands 3A2 and 3C2. By saturation mutagenesis in this 16-band region, Judd and colleagues, and subsequent studies (e.g., Lim and Synder 1974) defined 20 genes, of which 15 were vital. A number of other studies also concluded that the ratio between gene and band number was about one (Zhimulev 1999). It is now clear that, although the number of vital loci inDrosophila is indeed ∼5000, the use of lethal mutations to define genes results in a substantial underestimate; only about one-third of genes are vital.

The complete sequence of the tip of the X chromosome now gives us the chance to review the important study of Judd and colleagues with a molecular perspective (see also Judd's own recent historical review, Judd 1998). The region between the genesgiant and white studied by Judd et al. (thezeste-white region) is 360 Kb in length and is predicted to contain 45 genes (Fig. 2). It is indeed remarkable that conventional genetic analyses had identified 20 of these. Of these 20, 12 can be placed directly on the genetic map, by virtue of identity of sequence; the remaining eight genes, known only from lethal mutations, have not been sequenced independently.

Figure 2.

The zeste-white interval. The top is a reproduction of Figure 5 from Judd et al. (1972) showing the polytene chromosome region 3A–3C and the complementation groups discovered by mutational analysis. Below this projections are made onto the interval 2100 Kb to 2480 Kb of the EDGP sequence showing the correspondence between the genetic analysis and the genes known or predicted in this region from sequence analysis.

Open in new tabDownload PowerPointLink to figure

Unraveling the famous zeste-white region in the ultimate detail of its complete DNA sequence leaves major questions concerning the chromomeric structure of polytene chromosomes unanswered, of course. The banding pattern is attributable to aperiodicities in the packing ratio of the DNA, associated with proteins, in chromatin. Does this pattern have any functional significance whatsoever? No answer to this question can yet be given. How is the banding pattern determined? At one level the answer to this is obvious, by the DNA sequence. We have already described an inverted repeat sequence in the chromosomal DNA flanking the broad complex that could account for the unusual chromosomal banding pattern of this region (Benos et al. 2000). However, more subtle aspects of DNA sequence may define the domains of the majority of polytene chromosome bands, and the full answer to this problem will require considerable further analysis.

P-element Insertions

The majority of P-element screens to have been carried out to date have been performed on the autosomes. Spradling and colleagues (1999) have described their attempts to consolidate a number of suchP-element collections, including a large collection of lethal P-element insertions on the second chromosome (Török et al. 1993). Similarly, the EDGP have described a collection of lethal insertions on chromosome 3 (Deak et al. 1997). We have begun to generate a comparable collection ofP-element insertion mutants on the X chromosome in anticipation of their value for functional genomics. The initial group of mutants corresponds to ∼500 lethal insertions that have been mapped by hybridization of P-element probes to polytene chromosomes in situ. The characterization of this collection will be presented elsewhere. We have localized the insertion sites for 64P-element-induced lethal mutations that map to divisions 1–3, and determined the gene(s) whose function is likely to be affected by each insertion (Table 2). We have carried out a similar computational analysis on a collection of random EP-element insertions sequenced by the BDGP (Rørth et al. 1998). Forty-seven of these had been mapped to divisions 1–3 by in situ hybridization; this is a density of one element per 55 Kb, about twice that found for EP-elements in the Adh region (1/108 Kb). This difference in density is not due to the existence of major hotspots for insertion of EP-elements on theX chromosome tip, nor to a higher proportion of the insertions on the X tip being outwith genes (in both regions ∼47% of EP-element insertions are within genes).

From a total of 111 P-element insertions that we have located within the region analyzed, 41% fall in regions in which they are expected to affect the expression of genes already known, whereas 50% are expected to affect the expression of predicted genes. These expectations are based on the positions of the P-elementinsertion either within transcribed regions or within 5 Kb 5′ to these. Some insertions might affect two different genes, one on either side of the insertion (Table 2). Only 13 elements or clusters of elements map more distantly, 7–33 Kb 5′ to the nearest known or predicted gene (footnotes in Table 2; of these, five elements or groups were selected as lethal, but may or may not cause the lethality).

Comparison with the Joint Sequence

The determination of the sequence and gene annotation of chromosomal divisions 1–3 was completed and submitted to the EMBL-Bank by February 7, 2000, six weeks before the publication and release of the annotated Joint Sequence of the D. melanogaster genome in March 2000 (Adams et al. 2000). Although preexisting gene features were taken into account during the analysis of the Joint Sequence, these are essentially independent annotation experiments that can be compared. Moreover, direct comparison of the nucleotide sequence determined by the EDGP with the Joint Sequence, allows one to assess some of the strengths and weaknesses of the two different sequencing strategies. We have compared both individual gene predictions and the overall sequence between these two studies.

Overall Sequence Comparison

The Joint Sequence for region 1A–3C is found on nine GenBank entries (Fig. 3). We have compared it to the contiguous EDGP sequence using the MUMmer program ofDelcher et al. (1999) (Fig. 3). At the nucleotide level, the differences between our sequence and that of the Joint Sequence in this region are of two types: small indels and large (1 Kb or more) blocks of difference. Thirty large blocks of sequence are present in only one of the sequences. Ten of these blocks occur at identical nucleotide positions in both the Joint and EDGP sequences (Fig. 3). The null hypothesis is that these pairs of blocks are independent. As will be shown below, this is probably not true for all. Excluding these, the difference between the two studies at the nucleotide level is 3.03% (n = 2,568,355 common nucleotides). This figure may seem high, but over half (56%) of the EDGP sequence was from clones derived from a very different strain from that used for the Joint Sequence. We have partitioned this difference into that seen in known or predicted coding exons, known or predicted introns, and other sequences; the figures are 0.90%, 2.29%, and 3.98%, respectively.

Figure 3.

Sequence comparisons. A comparison of EDGP sequence of the tip of theX chromosome with that of the Drosophila Joint Sequence in the same region. The comparison was made using theMUMmer program (see Methods). The GenBank accession numbers corresponding to the Joint Sequence are shown on the left (AE003417–AE003425); note that this is part of a unitig (Myers et al. 2000). The blocks indicate regions of ⋝1 Kb present in one sequence but not the other. The position and length of each block of sequence ⋝1 Kb that is unique to one sequence is shown; each GenBank accession is numbered to the left of the unitig, the corresponding base position within the EDGP sequence is shown in italics to theright of the unitig. The EDGP sequence is numbered continuously. The length of each block of unique sequence is in parentheses. The nature of these sequence segments is shown in thecenter (note that a segment may include sequences in addition to those identified here). The segments corresponding to transposable elements are indicated in orange; those corresponding to known genes are red; a gray “neck” depicts a sequence interrupted by a large block of n's of length N (nN). The Greek superscripts (^α,^β,^γ,^ɛ,^δ) refer to the class of sequence difference (see text). Note that there are an additional nine transposable elements in the EDGP sequence that are not seen to differ in the Joint Sequence.

Open in new tabDownload PowerPointLink to figure

Most of the 30 blocks of sequence that appear to be absent from one or other sequence are either regions that have not been elucidated fully in the Joint Sequence, or correspond to transposable elements of variable location and/or length. In particular, 17 blocks in one or the other sequence correspond to recognizable transposable elements of variable length and/or location (^α in Fig. 3). These include two roo elements of different length found at the same position (nucleotide 572,960) in both sequences; five rooelements of variable location; and 10 single occurrences of other transposable element families at unique locations (BEL, 412, FB4, I, 412-like and mgd1 in the Joint Sequence, andDoc, Tirant, Burdock, and FB in the EDGP Sequence). It should be noted that two of the long runs ofn in the Joint Sequence correspond to transposable elements in the EDGP Sequence (see below). The 17 differences in transposable elements are not surprising, as the majority of the two sequences were derived from two quite different fruitfly strains. In the EDGP sequence we have identified 18 transposable elements or fragments of elements and at least 7 of these differ in position in the Joint Sequence.

Ten of the 30 blocks are long gaps in the Joint Sequence (^β, ^γ, ^ɛ in Fig. 3), represented in the GenBank accessions by long runs of n, with a total estimated length of 39,938 nucleotides. For four of the 10 gaps (^β), the length of the gap in the Joint Sequence is considerably larger than the corresponding region in the EDGP sequence; for example the run of 4722 n's at position 1,245,921 corresponds to 102 bp in the EDGP sequence. We presume the reason for this is that the gap in the Joint Sequence represents a transposable element. Indeed, two gaps (^ɛ) are caused by transposable elements: The 6353-bp gap at 2,294,896 corresponds to a 6062-bp Burdock element in the EDGP sequence, and the 8060-bp gap at 2,511,915 corresponds to a roo element in the EDGP sequence. Of the four remaining gaps (^γ), two are complex (at 237,007 bp and 556,147 bp) and cannot be explained simply; one corresponds to the ph-d/ph-p gene duplication (see below), and the final gap, at 2,011,597 bp will be discussed below.

The remaining three long blocks (^δ in Fig. 3) of the 30 that differ between the two sequences are informative, and will be discussed more fully. Two are only found in the EDGP sequence and are clearly the result of misassemblies in the Joint Sequence. The first of these is just 3′ to the Actn gene and is 4.7-Kb long; the probable explanation for it is that the Joint Sequence has failed to properly assemble a duplicated sequence that includes a partial duplication of the predicted gene EG:133E12.4. This duplication was first indicated by the matches of EST sequences (e.g., EMBL accession no. AA202518, EMBL accession no. AA696909) to both an exon of EG:133E12.4 and to a region between this gene andActn. The duplication is 4777 bp in length and the two copies are only mismatched over a 77-bp internal gap (1.5% mismatch). The second is in the region of the duplicate gene pair ph-d andph-p; the Joint Sequence has an incorrect model forph-p. That this region includes a long tandem repeat is known from the work of Deatrick et al. (1991).

The third region, at 2,011,597, is more complex. There is an 18.5 Kb region (of which 7.1 Kb are n's) in the Joint Sequence absent from the EDGP sequence; this sequence is not in the shotgun sequence of either relevant EDGP clone, cosmid 82C7, or BACH48C10. In addition, there is a 10.3-Kb sequence at the junction of these clones in the EDGP sequence that is absent from the Joint Sequence. Finally, 11 Kb of cosmid 82C7 is in the opposite orientation when compared to BACH48C10; note that the cosmid and BAC DNAs are from different strains (see Methods).

These three major sequence differences could be caused by polymorphisms; all occur within regions of EDGP cosmid sequence. However we consider that the hypothesis of misassembly, at least for the Actn and ph-d/ph-p region differences, is the more likely. The current “finishing” of the Joint Sequence by the BDGP should settle these problems.

Repeated regions are well known to present a problem to the software used to build long contiguous regions of sequence, and there is evidence of this in at least two regions of the Joint Sequence. It is interesting that in both cases the assembler appears to have had difficulties with tandem near repeats of quite long regions. Using statistical criteria, the software that assembled the Joint Sequence was able to identify and filter out the highly repetitive sequences, based on their higher than expected representation (Myers et al. 2000). However, the low copy repetitive sequences (such as the tandemly duplicated regions in these two cases) are difficult to identify by these methods. If this comparison of the X tip is typical of the genome as a whole, then it indicates some 90 misassemblies in the euchromatic sequence of the Joint Sequence.

The differences revealed by this comparison of the genomic sequence from the two projects includes both differences in sequencing method (clone-based in the case of the EDGP, and shotgun in the case of the Joint Sequence) and differences in strain from which the DNA was derived. Even when the sequenced DNA is from the same strain, but isolated some years apart, there are differences in sequence and transposable elements. For example, Myers et al. (2000) compared theAdh region sequenced by the BDGP using predominantly P1 clones (Ashburner et al. 1999) with that from the Joint Sequence. Although the differences are smaller than in the comparison made in this study, they are qualitatively very similar.

There are clear differences in gene predictions between the EDGP and Joint Sequence projects, both in the existence of genes and in the precise models of genes predicted in common. Again this is not too surprising, given that the Joint Sequence was annotated very largely by automatic methods, whereas the EDGP had the luxury of time to make a more careful study of each gene model. These differences point out that we have a long way to go before the annotation of eukaryotic sequences can be left entirely in the hands of computer programs (Ashburner 2000;Lewis et al. 2000).

This analysis has, for obvious reasons, concentrated on the differences between the two available sequences of this chromosome region. This must not obscure the fact that in general the two analyses are in remarkable agreement, and point to the overall utility of the “complete” genomic sequence now available for D. melanogaster.

Clone Libraries and Map Construction

DNA from two strains has been sequenced. About 44% of the sequence is from BAC clones derived from the same strain as that sequenced by the BDGP and by Celera; in contrast, the cosmid clones sequenced were from a different strain (Fig. 1). The relationship between these strains cannot be determined. Both strains were free ofP-elements.

The cosmid library used for the construction of the Xchromosome physical map was derived from a wild-type (Canton-S) strain and described in detail by Sidén-Kiamos et al. (1990). It has an estimated average insert size of 35 Kb and contains ∼18,000 clones providing a fourfold coverage of the genome. The library is available on high density double spotted filters from the MRC HGMP Resource Centre (http://www.hgmp.mrc.ac.uk/Biology/Bio.html).

Three BAC clone libraries were used; each was constructed from DNA from the y²; cn bw sp isogenic strain. Two BAC libraries were made at CEPH (Centre d'Etude du Polymorphisme Humaine). One (BACN clones) was prepared with NdeII inserts and the other (BACH clones) with HindIII inserts, both in the vector pBeloBACII. These two libraries were made with pools of size-fractioned DNA that gave mean insert sizes of up to 90 Kb. The 23,400 clones gave ∼10-fold coverage of the genome. The third library was ofEcoRI digested DNA (BACR clones) and was constructed in the vector pBACe.3.6 by Aaron Mammoser and Kazutoyo Oseogawa at the Roswell Park Cancer Institute (Buffalo, NY) in collaboration with the BDGP (Hoskins et al. 2000). This library gave an ∼17-fold coverage of the genome with an average insert size of 165 Kb.

Sequencing

Cosmids and BACs were sequenced by a two-stage approach involving random sequencing of sub-clones followed by directed sequencing to resolve problems. DNA from cosmids and BACs was sonicated and fragments of 1.4–2 Kb were cloned into either M13 or pUC18 vectors. Clones were sequenced using dye-terminator chemistry and loaded on ABI373 or ABI377 automated sequencing machines. Sequence base calling and contig assembly was accomplished using Phred/Phrap software (Ewing and Green 1998; Ewing et al. 1998) and editing took place in either Consed (Gordon et al. 1998) or Gap4(Bonfield et al. 1995). Gaps were filled using a combination of custom primer walking and PCR.

Cosmid and BAC DNAs were nebulized and end repaired. Following agarose gel purification, fragments of ∼1500 nucleotides were ligated to linearized vector (pTZ19R or pCR-BluntII) and cloned in the KK2186 strain of Escherichia coli. Bacterial clones were picked at random and cultured overnight. Plasmid DNAs were prepared by an alkaline lysis method and purified using the QIAprep 96 Turbo Miniprep kit (QIAGEN). Insert DNA were sequenced from both ends using universal primers. Cycle sequencing was performed with labeled terminators using AmpliTaq and the Big Dye Terminator Cycle Sequencing Ready Reaction kit (Applied Biosystems).

The Heidelberg group employed the RANDI strategy that combines the advantages of RANdom and DIrected approaches. It involves systematic simultaneous sequencing on both strands from clones of combined libraries without cloning gaps. The random library fragments were generated by separate partial digestion with two four-cutter restriction enzymes (Tsp, Sau3A), gel-purified and ligated into plasmid vector. In parallel, BAC or cosmid DNA was completely digested with EcoRI (or HindIII) and fragments were isolated from agarose gel and inserted into the pUC vector. Their sequences served as a “scaffold” in the assembly of the complete sequence of the BAC genomic insert and also as templates for primer walking in the finishing stage. Cycle sequencing of plasmid DNA was performed with the AmpliTaqFS core kit (Applied Biosystems), using forward and reverse primers labeled with FITC or CY5. An MJ Research PT-200 cycler was used for 25 cycles (97°C, 15 sec; 55°C, 30 sec; 68°C, 30 sec). Reactions were loaded off-gel on the 72-clone porous-membrane combs, applied to 60-cm long polyacrylamide gels (4.5% Hydrolink Long Ranger gel solution, FMC) and analyzed on the ARAKIS sequencing system with array detectors, developed at EMBL (Erfle et al. 1997). This system allows simultaneous on-line sequencing of both strands (doublex sequencing), with the two sequencing products obtained in a single sequencing reaction, each labeled with a different fluorescent dye (Wiemann et al. 1995). Up to 2000 bases are thus obtained simultaneously in one sequencing reaction, which represents an efficient system for identifying large numbers of long sequences in one run. Raw sequencing data were evaluated, analyzed, and the consensus sequence assembled, using the software packages (LaneTracker and GeneSkipper) developed at EMBL. Remaining sequencing gaps were covered by primer walking (Voss et al. 1993). Direct cosmid or BAC DNA sequencing was carried out essentially as described elsewhere (Benes et al. 1997).

P-element Stocks and Mapping

A large-scale screen for insertions of the enhancer trap vectorP{lacW} (Bier et al. 1989) in essential Xchromosome genes has been performed in H. Jäckle's laboratory (Peter et al., in prep.). Females homozygous for a male sterile insertion of the P{lacW} element in chromosome 2 were crossed en masse to w/Y; wg Sp/CyO; P{ry ⁺=delta2–3}(99B)males. In the next generation five homozygous FM6 females were mated to two w/Y; P{lacW}/CyO; P{ry ⁺=delta2–3}(99B)/+ males. F2 daughters in which the CyO andP{lacW} chromosomes had cosegregated were individually mated to Fm7c/Y males. Lines that produced only FM6sons in the F3 generation were kept as candidates for a lethal insertion. If these re-tested, then the lethal insertion was kept in stock balanced with FM7c.

P{lacW} insertion sites were mapped by either plasmid rescue or inverse PCR. DNA from adult flies was isolated using a QIAGEN column, digested overnight with an appropriate restriction enzyme, and then ligated under conditions favoring intramolecular joining. For plasmid rescue, E. coli cells were electroporated with the DNA and plated for the selection of ampicillin resistant colonies. These were used to inoculate small scale overnight cultures from which plasmid DNA was then isolated. Cycle sequencing was performed with a primer complementary to the 31-bp inverted repeat of the P-element on an ABI373 DNA sequencer using dye terminator technology. In the case of inverse PCR, we followed essentially the protocol from the BDGP. We used their primers Plac1 and Plac4 for the amplification of 5′ sequences and primers Pry4 and Plw3–1 for the amplification of 3′ sequences, respectively. Sequencing was done as before with primer SP1 for 5′ and primer SP6 for 3′ analysis.

Sequence Analysis

Sequences were analyzed by the EDGP on a clone-by-clone basis; i.e., only fully sequenced clones (cosmids or BACs) were included. The overall analysis scheme is similar to that adopted by other genome projects (e.g., C. elegans Sequencing Consortium 1998).

tRNA genes were identified by tRNAscan-SE program, v. 1.0 (Lowe and Eddy 1997). Candidate protein coding genes were predicted independently by GENEFINDER version 0.84 (P. Green, unpubl.) and the publicly available Genscan version 1.0 (Burge and Karlin 1997). These two programs employ fundamentally different algorithms and complemented each other on gene discovery. GENSCAN andGENEFINDER had been trained on a vertebrate gene set and aDrosophila-specific set (compiled by G. Helt, pers. comm.), respectively. We measured the accuracy of prediction of the two programs with already known Drosophila genes and we found them to be comparable. However, each of them performed better on a different set of genes. As expected, Drosophila-trained GENEFINDER showed a preference for genes with fewer exons and smaller introns when compared to the vertebrate-trainedGENSCAN.

Additional supporting evidence for the predicted genes, as well as indications of their function, was obtained by similarity searches against SWISS-PROT and TrEMBL protein databases (Bairoch and Apweiler 2000), Drosophila nucleic acid sequences (derived from EMBL-Bank), and Drosophila EST sets, generated by the BDGP (Rubin et al. 2000b) and by Andrews et al. (2000). (Note that the annotation of version 1 of the Joint Sequence did not use the entire BDGP EST data set; in particular 4,654 3′ ESTs, out of a total of 86,121, were not used [S. Lewis, pers. comm.]). EST alignments were also used to fine-tune the intron/exon boundaries of the predicted genes. Simple repetitive sequences were filtered out byTANDEM, INVERTED, andQUICKTANDEM programs (R. Durbin, pers. comm..) whereas repeats of higher complexity were screened out using similarity searches against Drosophila repetitive and transposable element databases (see below). For protein and nucleotide database searches we used BLASTX and BLASTN, v. 1.4.9. (Altschul et al. 1990), respectively.

Finally, protein domains/motifs of the predicted genes were identified by PPSEARCH and HMMER (v. 2.1.1) programs, scanning the PROSITE and PFAM databases, respectively. PROSITE output was further filtered using the EMOTIF program (Nevill-Manning et al. 1998).

All data generated by the automatic computational analysis described above were parsed into an ACeDB-based database (http://www.acedb.org/), XDrosDB, tailored to the needs of the EDGP. The combined data were manually examined/analyzed using ACeDB software. During this analysis we disregarded genes with a GENEFINDER score <50, if there was no other supporting evidence for them (i.e., protein similarity and/or EST matches). This cutoff is stricter than the one used by the BDGP (cutoff = 20) for the analysis of the Adhregion (Ashburner et al. 1999); and, presumably, increases the number of rejected genes (false negatives). However, we chose to set it this high to avoid overpredicting genes (false positives).

During the initial phase of our work, we, in collaboration with the BDGP, created and subsequently curated three datasets. One consisted of 1332 D. melanogaster coding sequences from genes that have been previously studied genetically and/or biochemically. This is a nonredundant set, i.e., only one copy of each gene is included in it. In case a gene appears in multiple entries in the public databases (e.g., alternatively transcribed, submitted from more than one laboratory, etc.), we manually selected one copy (usually the best documented or longest open reading frame). We used this dataset to test the accuracy of the two chosen gene prediction programs (GENEFINDER, GENSCAN), as well as a source for hexanucleotides score calculation (GENEFINDER). This dataset has been subsequently expanded/updated to include genes identified by Drosophila genome projects (EDGP, BDGP, and Celera), with the help of Leyla Bayraktaroglou (FlyBase at Harvard). Both the original and expanded versions, together with information about their history, can be found at:ftp://ftp.ebi.ac.uk/pub/databases/edgp/sequence_sets/ or fromhttp://fruitfly.berkeley.edu/.

Similarly, a nonredundant collection of 47 D. melanogastertransposable elements and another consisting of 96 miscellaneous repetitive sequences were also assembled during the initial phase of our project. These datasets were used to identify complex repetitive regions, as described previously. They are also available from the same ftp site or from the BDGP site.

For clarity, we use the term “Joint Sequence” to refer to v1.0 of the complete sequence of the genome of D. melanogaster (Adams et al. 2000) released on March 24, 2000 by Celera. Comparisons of predicted, or known, protein sequences from the EDGP project with those from the Joint Sequence were done by CLUSTALW using the protein sequences of release 1.0 of the Joint Sequence (http://www.fruitfly.org/sequence/sequence_db/aa_gadfly.dros of March 21, 2000). These comparisons were then analyzed by hand. The comparison of the entire sequence of the X chromosome tip with the sequence of the same region from the Joint Sequence was done using theMUMmer program (Delcher et al. 1999), which aligns long genomic regions by finding corresponding maximal unique matches. Nine separate alignments were done using the following GenBank accession nos.: AE003417, AE003418, AE003419, AE003420, AE003421, AE003422,AE003423, AE003424, and AE003425, each being matched against the entire EDGP sequence. The resulting alignments were analyzed by hand to find regions where the discrepancies between the sequences were large. Figure 3 was drawn by hand and is a graphic depiction of the alignment produced by MUMmer. Large segments absent from one of the sequences have been highlighted.

The results presented in this study were obtained by or before February 7, 2000. However, if we had repeated the same analysis today we would have assigned function (by protein similarity) to 23 more of the predicted genes (raising the percentage of the genes with significant protein similarities to 66% of the 206 newly identified genes).

Supplementary data are available fromftp://ebi.ac.uk/pub/databases/edgp/EDGP-GenomeResearch_suppdata_2001.