J. Craig Venter,^1* Mark D. Adams,¹ Eugene W. Myers,¹ Peter W. Li,¹ Richard J. Mural,¹ Granger G. Sutton,¹ Hamilton O. Smith,¹ Mark Yandell,¹ Cheryl A. Evans,¹ Robert A. Holt,¹ Jeannine D. Gocayne,¹ Peter Amanatides,¹ Richard M. Ballew,¹ Daniel H. Huson,¹ Jennifer Russo Wortman,¹ Qing Zhang,¹ Chinnappa D. Kodira,¹ Xiangqun H. Zheng,¹ Lin Chen,¹ Marian Skupski,¹ Gangadharan Subramanian,¹ Paul D. Thomas,¹ Jinghui Zhang,¹ George L. Gabor Miklos,² Catherine Nelson,³ Samuel Broder,¹ Andrew G. Clark,⁴ Joe Nadeau,⁵ Victor A. McKusick,⁶ Norton Zinder,⁷ Arnold J. Levine,⁷ Richard J. Roberts,⁸ Mel Simon,⁹ Carolyn Slayman,¹⁰ Michael Hunkapiller,¹¹ Randall Bolanos,¹ Arthur Delcher,¹ Ian Dew,¹ Daniel Fasulo,¹ Michael Flanigan,¹ Liliana Florea,¹ Aaron Halpern,¹ Sridhar Hannenhalli,¹ Saul Kravitz,¹ Samuel Levy,¹ Clark Mobarry,¹ Knut Reinert,¹ Karin Remington,¹ Jane Abu-Threideh,¹ Ellen Beasley,¹ Kendra Biddick,¹ Vivien Bonazzi,¹ Rhonda Brandon,¹ Michele Cargill,¹ Ishwar Chandramouliswaran,¹ Rosane Charlab,¹ Kabir Chaturvedi,¹ Zuoming Deng,¹ Valentina Di Francesco,¹ Patrick Dunn,¹ Karen Eilbeck,¹ Carlos Evangelista,¹ Andrei E. Gabrielian,¹ Weiniu Gan,¹ Wangmao Ge,¹ Fangcheng Gong,¹ Zhiping Gu,¹ Ping Guan,¹ Thomas J. Heiman,¹ Maureen E. Higgins,¹ Rui-Ru Ji,¹ Zhaoxi Ke,¹ Karen A. Ketchum,¹ Zhongwu Lai,¹ Yiding Lei,¹ Zhenya Li,¹ Jiayin Li,¹ Yong Liang,¹ Xiaoying Lin,¹ Fu Lu,¹ Gennady V. Merkulov,¹ Natalia Milshina,¹ Helen M. Moore,¹ Ashwinikumar K Naik,¹ Vaibhav A. Narayan,¹ Beena Neelam,¹ Deborah Nusskern,¹ Douglas B. Rusch,¹ Steven Salzberg,¹² Wei Shao,¹ Bixiong Shue,¹ Jingtao Sun,¹ Zhen Yuan Wang,¹ Aihui Wang,¹ Xin Wang,¹ Jian Wang,¹ Ming-Hui Wei,¹ Ron Wides,¹³ Chunlin Xiao,¹ Chunhua Yan,¹ Alison Yao,¹ Jane Ye,¹ Ming Zhan,¹ Weiqing Zhang,¹ Hongyu Zhang,¹ Qi Zhao,¹ Liansheng Zheng,¹ Fei Zhong,¹ Wenyan Zhong,¹ Shiaoping C. Zhu,¹ Shaying Zhao,¹² Dennis Gilbert,¹ Suzanna Baumhueter,¹ Gene Spier,¹ Christine Carter,¹ Anibal Cravchik,¹ Trevor Woodage,¹ Feroze Ali,¹ Huijin An,¹ Aderonke Awe,¹ Danita Baldwin,¹ Holly Baden,¹ Mary Barnstead,¹ Ian Barrow,¹ Karen Beeson,¹ Dana Busam,¹ Amy Carver,¹ Angela Center,¹ Ming Lai Cheng,¹ Liz Curry,¹ Steve Danaher,¹ Lionel Davenport,¹ Raymond Desilets,¹ Susanne Dietz,¹ Kristina Dodson,¹ Lisa Doup,¹ Steven Ferriera,¹ Neha Garg,¹ Andres Gluecksmann,¹ Brit Hart,¹ Jason Haynes,¹ Charles Haynes,¹ Cheryl Heiner,¹ Suzanne Hladun,¹ Damon Hostin,¹ Jarrett Houck,¹ Timothy Howland,¹ Chinyere Ibegwam,¹ Jeffery Johnson,¹ Francis Kalush,¹ Lesley Kline,¹ Shashi Koduru,¹ Amy Love,¹ Felecia Mann,¹ David May,¹ Steven McCawley,¹ Tina McIntosh,¹ Ivy McMullen,¹ Mee Moy,¹ Linda Moy,¹ Brian Murphy,¹ Keith Nelson,¹ Cynthia Pfannkoch,¹ Eric Pratts,¹ Vinita Puri,¹ Hina Qureshi,¹ Matthew Reardon,¹ Robert Rodriguez,¹ Yu-Hui Rogers,¹ Deanna Romblad,¹ Bob Ruhfel,¹ Richard Scott,¹ Cynthia Sitter,¹ Michelle Smallwood,¹ Erin Stewart,¹ Renee Strong,¹ Ellen Suh,¹ Reginald Thomas,¹ Ni Ni Tint,¹ Sukyee Tse,¹ Claire Vech,¹ Gary Wang,¹ Jeremy Wetter,¹ Sherita Williams,¹ Monica Williams,¹ Sandra Windsor,¹ Emily Winn-Deen,¹ Keriellen Wolfe,¹ Jayshree Zaveri,¹ Karena Zaveri,¹ Josep F. Abril,¹⁴ Roderic Guigó,¹⁴ Michael J. Campbell,¹ Kimmen V. Sjolander,¹ Brian Karlak,¹ Anish Kejariwal,¹ Huaiyu Mi,¹ Betty Lazareva,¹ Thomas Hatton,¹ Apurva Narechania,¹ Karen Diemer,¹ Anushya Muruganujan,¹ Nan Guo,¹ Shinji Sato,¹ Vineet Bafna,¹ Sorin Istrail,¹ Ross Lippert,¹ Russell Schwartz,¹ Brian Walenz,¹ Shibu Yooseph,¹ David Allen,¹ Anand Basu,¹ James Baxendale,¹ Louis Blick,¹ Marcelo Caminha,¹ John Carnes-Stine,¹ Parris Caulk,¹ Yen-Hui Chiang,¹ My Coyne,¹ Carl Dahlke,¹ Anne Deslattes Mays,¹ Maria Dombroski,¹ Michael Donnelly,¹ Dale Ely,¹ Shiva Esparham,¹ Carl Fosler,¹ Harold Gire,¹ Stephen Glanowski,¹ Kenneth Glasser,¹ Anna Glodek,¹ Mark Gorokhov,¹ Ken Graham,¹ Barry Gropman,¹ Michael Harris,¹ Jeremy Heil,¹ Scott Henderson,¹ Jeffrey Hoover,¹ Donald Jennings,¹ Catherine Jordan,¹ James Jordan,¹ John Kasha,¹ Leonid Kagan,¹ Cheryl Kraft,¹ Alexander Levitsky,¹ Mark Lewis,¹ Xiangjun Liu,¹ John Lopez,¹ Daniel Ma,¹ William Majoros,¹ Joe McDaniel,¹ Sean Murphy,¹ Matthew Newman,¹ Trung Nguyen,¹ Ngoc Nguyen,¹ Marc Nodell,¹ Sue Pan,¹ Jim Peck,¹ Marshall Peterson,¹ William Rowe,¹ Robert Sanders,¹ John Scott,¹ Michael Simpson,¹ Thomas Smith,¹ Arlan Sprague,¹ Timothy Stockwell,¹ Russell Turner,¹ Eli Venter,¹ Mei Wang,¹ Meiyuan Wen,¹ David Wu,¹ Mitchell Wu,¹ Ashley Xia,¹ Ali Zandieh,¹ Xiaohong Zhu¹

A 2.91-billion base pair (bp) consensus sequence of the euchromatic portion of the human genome was generated by the whole-genome shotgun sequencing method. The 14.8-billion bp DNA sequence was generated over 9 months from 27,271,853 high-quality sequence reads (5.11-fold coverage of the genome) from both ends of plasmid clones made from the DNA of five individuals. Two assembly strategies--a whole-genome assembly and a regional chromosome assembly--were used, each combining sequence data from Celera and the publicly funded genome effort. The public data were shredded into 550-bp segments to create a 2.9-fold coverage of those genome regions that had been sequenced, without including biases inherent in the cloning and assembly procedure used by the publicly funded group. This brought the effective coverage in the assemblies to eightfold, reducing the number and size of gaps in the final assembly over what would be obtained with 5.11-fold coverage. The two assembly strategies yielded very similar results that largely agree with independent mapping data. The assemblies effectively cover the euchromatic regions of the human chromosomes. More than 90% of the genome is in scaffold assemblies of 100,000 bp or more, and 25% of the genome is in scaffolds of 10 million bp or larger. Analysis of the genome sequence revealed 26,588 protein-encoding transcripts for which there was strong corroborating evidence and an additional ~12,000 computationally derived genes with mouse matches or other weak supporting evidence. Although gene-dense clusters are obvious, almost half the genes are dispersed in low G+C sequence separated by large tracts of apparently noncoding sequence. Only 1.1% of the genome is spanned by exons, whereas 24% is in introns, with 75% of the genome being intergenic DNA. Duplications of segmental blocks, ranging in size up to chromosomal lengths, are abundant throughout the genome and reveal a complex evolutionary history. Comparative genomic analysis indicates vertebrate expansions of genes associated with neuronal function, with tissue-specific developmental regulation, and with the hemostasis and immune systems. DNA sequence comparisons between the consensus sequence and publicly funded genome data provided locations of 2.1 million single-nucleotide polymorphisms (SNPs). A random pair of human haploid genomes differed at a rate of 1 bp per 1250 on average, but there was marked heterogeneity in the level of polymorphism across the genome. Less than 1% of all SNPs resulted in variation in proteins, but the task of determining which SNPs have functional consequences remains an open challenge.

¹ Celera Genomics, 45 West Gude Drive, Rockville, MD 20850, USA.
² GenetixXpress, 78 Pacific Road, Palm Beach, Sydney 2108, Australia.
³ Berkeley Drosophila Genome Project, University of California, Berkeley, CA 94720, USA.
⁴ Department of Biology, Penn State University, 208 Mueller Lab, University Park, PA 16802, USA.
⁵ Department of Genetics, Case Western Reserve University School of Medicine, BRB-630, 10900 Euclid Avenue, Cleveland, OH 44106, USA.
⁶ Johns Hopkins University School of Medicine, Johns Hopkins Hospital, 600 North Wolfe Street, Blalock 1007, Baltimore, MD 21287-4922, USA.
⁷ Rockefeller University, 1230 York Avenue, New York, NY 10021-6399, USA.
⁸ New England BioLabs, 32 Tozer Road, Beverly, MA 01915, USA.
⁹ Division of Biology, 147-75, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA.
¹⁰ Yale University School of Medicine, 333 Cedar Street, P.O. Box 208000, New Haven, CT 06520-8000, USA.
¹¹ Applied Biosystems, 850 Lincoln Centre Drive, Foster City, CA 94404, USA.
¹² The Institute for Genomic Research, 9712 Medical Center Drive, Rockville, MD 20850, USA.
¹³ Faculty of Life Sciences, Bar-Ilan University, Ramat-Gan, 52900 Israel.
¹⁴ Grup de Recerca en Informàtica Mèdica, Institut Municipal d'Investigació Mèdica, Universitat Pompeu Fabra, 08003-Barcelona, Catalonia, Spain.
^*   To whom correspondence should be addressed. E-mail: humangenome{at}celera.com

The modern history of DNA sequencing began in 1977, when Sanger reported his method for determining the order of nucleotides of DNA using chain-terminating nucleotide analogs (3). In the same year, the first human gene was isolated and sequenced (4). In 1986, Hood and co-workers (5) described an improvement in the Sanger sequencing method that included attaching fluorescent dyes to the nucleotides, which permitted them to be sequentially read by a computer. The first automated DNA sequencer, developed by Applied Biosystems in California in 1987, was shown to be successful when the sequences of two genes were obtained with this new technology (6). From early sequencing of human genomic regions (7), it became clear that cDNA sequences (which are reverse-transcribed from RNA) would be essential to annotate and validate gene predictions in the human genome. These studies were the basis in part for the development of the expressed sequence tag (EST) method of gene identification (8), which is a random selection, very high throughput sequencing approach to characterize cDNA libraries. The EST method led to the rapid discovery and mapping of human genes (9). The increasing numbers of human EST sequences necessitated the development of new computer algorithms to analyze large amounts of sequence data, and in 1993 at The Institute for Genomic Research (TIGR), an algorithm was developed that permitted assembly and analysis of hundreds of thousands of ESTs. This algorithm permitted characterization and annotation of human genes on the basis of 30,000 EST assemblies (10).

The complete 49-kbp bacteriophage lambda genome sequence was determined by a shotgun restriction digest method in 1982 (11). When considering methods for sequencing the smallpox virus genome in 1991 (12), a whole-genome shotgun sequencing method was discussed and subsequently rejected owing to the lack of appropriate software tools for genome assembly. However, in 1994, when a microbial genome-sequencing project was contemplated at TIGR, a whole-genome shotgun sequencing approach was considered possible with the TIGR EST assembly algorithm. In 1995, the 1.8-Mbp Haemophilus influenzae genome was completed by a whole-genome shotgun sequencing method (13). The experience with several subsequent genome-sequencing efforts established the broad applicability of this approach (14, 15).

A key feature of the sequencing approach used for these megabase-size and larger genomes was the use of paired-end sequences (also called mate pairs), derived from subclone libraries with distinct insert sizes and cloning characteristics. Paired-end sequences are sequences 500 to 600 bp in length from both ends of double-stranded DNA clones of prescribed lengths. The success of using end sequences from long segments (18 to 20 kbp) of DNA cloned into bacteriophage lambda in assembly of the microbial genomes led to the suggestion (16) of an approach to simultaneously map and sequence the human genome by means of end sequences from 150-kbp bacterial artificial chromosomes (BACs) (17, 18). The end sequences spanned by known distances provide long-range continuity across the genome. A modification of the BAC end-sequencing (BES) method was applied successfully to complete chromosome 2 from the Arabidopsis thaliana genome (19).

In 1997, Weber and Myers (20) proposed whole-genome shotgun sequencing of the human genome. Their proposal was not well received (21). However, by early 1998, as less than 5% of the genome had been sequenced, it was clear that the rate of progress in human genome sequencing worldwide was very slow (22), and the prospects for finishing the genome by the 2005 goal were uncertain.

In early 1998, PE Biosystems (now Applied Biosystems) developed an automated, high-throughput capillary DNA sequencer, subsequently called the ABI PRISM 3700 DNA Analyzer. Discussions between PE Biosystems and TIGR scientists resulted in a plan to undertake the sequencing of the human genome with the 3700 DNA Analyzer and the whole-genome shotgun sequencing techniques developed at TIGR (23). Many of the principles of operation of a genome-sequencing facility were established in the TIGR facility (24). However, the facility envisioned for Celera would have a capacity roughly 50 times that of TIGR, and thus new developments were required for sample preparation and tracking and for whole-genome assembly. Some argued that the required 150-fold scale-up from the H. influenzae genome to the human genome with its complex repeat sequences was not feasible (25). The Drosophila melanogaster genome was thus chosen as a test case for whole-genome assembly on a large and complex eukaryotic genome. In collaboration with Gerald Rubin and the Berkeley Drosophila Genome Project, the nucleotide sequence of the 120-Mbp euchromatic portion of the Drosophila genome was determined over a 1-year period (26-28). The Drosophila genome-sequencing effort resulted in two key findings: (i) that the assembly algorithms could generate chromosome assemblies with highly accurate order and orientation with substantially less than 10-fold coverage, and (ii) that undertaking multiple interim assemblies in place of one comprehensive final assembly was not of value.

These findings, together with the dramatic changes in the public genome effort subsequent to the formation of Celera (29), led to a modified whole-genome shotgun sequencing approach to the human genome. We initially proposed to do 10-fold sequence coverage of the genome over a 3-year period and to make interim assembled sequence data available quarterly. The modifications included a plan to perform random shotgun sequencing to ~5-fold coverage and to use the unordered and unoriented BAC sequence fragments and subassemblies published in GenBank by the publicly funded genome effort (30) to accelerate the project. We also abandoned the quarterly announcements in the absence of interim assemblies to report.

Although this strategy provided a reasonable result very early that was consistent with a whole-genome shotgun assembly with eightfold coverage, the human genome sequence is not as finished as the Drosophila genome was with an effective 13-fold coverage. However, it became clear that even with this reduced coverage strategy, Celera could generate an accurately ordered and oriented scaffold sequence of the human genome in less than 1 year. Human genome sequencing was initiated 8 September 1999 and completed 17 June 2000. The first assembly was completed 25 June 2000, and the assembly reported here was completed 1 October 2000. Here we describe the whole-genome random shotgun sequencing effort applied to the human genome. We developed two different assembly approaches for assembling the ~3 billion bp that make up the 23 pairs of chromosomes of the Homo sapiens genome. Any GenBank-derived data were shredded to remove potential bias to the final sequence from chimeric clones, foreign DNA contamination, or misassembled contigs. Insofar as a correctly and accurately assembled genome sequence with faithful order and orientation of contigs is essential for an accurate analysis of the human genetic code, we have devoted a considerable portion of this manuscript to the documentation of the quality of our reconstruction of the genome. We also describe our preliminary analysis of the human genetic code on the basis of computational methods. Figure 1 (see fold-out chart associated with this issue; files for each chromosome can be found in Web fig. 1 on Science Online at www.sciencemag.org/cgi/content/full/291/5507/1304/DC2) provides a graphical overview of the genome and the features encoded in it. The detailed manual curation and interpretation of the genome are just beginning.

To aid the reader in locating specific analytical sections, we have divided the paper into seven broad sections. A summary of the major results appears at the beginning of each section.

1 Sources of DNA and Sequencing Methods

2  Genome Assembly Strategy and Characterization

3 Gene Prediction and Annotation

4 Genome Structure

5 Genome Evolution

6  A Genome-Wide Examination of Sequence Variations

7  An Overview of the Predicted Protein- Coding Genes in the Human Genome

8 Conclusions

1 Sources of DNA and Sequencing Methods

Various policies of the United States and the World Medical Association, specifically the Declaration of Helsinki, offer recommendations for conducting experiments with human subjects. We convened an Institutional Review Board (IRB) (31) that helped us establish the protocol for obtaining and using human DNA and the informed consent process used to enroll research volunteers for the DNA-sequencing studies reported here. We adopted several steps and procedures to protect the privacy rights and confidentiality of the research subjects (donors). These included a two-stage consent process, a secure random alphanumeric coding system for specimens and records, circumscribed contact with the subjects by researchers, and options for off-site contact of donors. In addition, Celera applied for and received a Certificate of Confidentiality from the Department of Health and Human Services. This Certificate authorized Celera to protect the privacy of the individuals who volunteered to be donors as provided in Section 301(d) of the Public Health Service Act 42 U.S.C. 241(d).

Celera and the IRB believed that the initial version of a completed human genome should be a composite derived from multiple donors of diverse ethnic backgrounds Prospective donors were asked, on a voluntary basis, to self-designate an ethnogeographic category (e.g., African-American, Chinese, Hispanic, Caucasian, etc.). We enrolled 21 donors (32).

Three basic items of information from each donor were recorded and linked by confidential code to the donated sample: age, sex, and self-designated ethnogeographic group. From females, ~130 ml of whole, heparinized blood was collected. From males, ~130 ml of whole, heparinized blood was collected, as well as five specimens of semen, collected over a 6-week period. Permanent lymphoblastoid cell lines were created by Epstein-Barr virus immortalization. DNA from five subjects was selected for genomic DNA sequencing: two males and three females--one African-American, one Asian-Chinese, one Hispanic-Mexican, and two Caucasians (see Web fig. 2 on Science Online at www.sciencemag.org/cgi/content/291/5507/1304/DC1). The decision of whose DNA to sequence was based on a complex mix of factors, including the goal of achieving diversity as well as technical issues such as the quality of the DNA libraries and availability of immortalized cell lines.

1.1 Library construction and sequencing

Table 1. Celera-generated data input into assembly. Individual Number of reads for different insert libraries Total number of base pairs 2 kbp 10 kbp 50 kbp Total No. of sequencing reads A 0 0 2,767,357 2,767,357 1,502,674,851 B 11,736,757 7,467,755 66,930 19,271,442 10,464,393,006 C 853,819 881,290 0 1,735,109 942,164,187 D 952,523 1,046,815 0 1,999,338 1,085,640,534 F 0 1,498,607 0 1,498,607 813,743,601 Total 13,543,099 10,894,467 2,834,287 27,271,853 14,808,616,179 Fold sequence coverage A 0 0 0.52 0.52 (2.9-Gb genome) B 2.20 1.40 0.01 3.61 C 0.16 1.17 0 0.32 D 0.18 0.20 0 0.37 F 0 0.28 0 0.28 Total 2.54 2.04 0.53 5.11 Fold clone coverage A 0 0 18.39 18.39 B 2.96 11.26 0.44 14.67 C 0.22 1.33 0 1.54 D 0.24 1.58 0 1.82 F 0 2.26 0 2.26 Total 3.42 16.43 18.84 38.68 Insert size* (mean) Average 1,951 bp 10,800 bp 50,715 bp Insert size* (SD) Average 6.10% 8.10% 14.90% % Mates Average 74.50 80.80 75.60 * Insert size and SD are calculated from assembly of mates on contigs. % Mates is based on laboratory tracking of sequencing runs.

In designing the DNA-sequencing process, we focused on developing a simple system that could be implemented in a robust and reproducible manner and monitored effectively (Fig. 2) (34).

Current sequencing protocols are based on the dideoxy sequencing method (35), which typically yields only 500 to 750 bp of sequence per reaction. This limitation on read length has made monumental gains in throughput a prerequisite for the analysis of large eukaryotic genomes. We accomplished this at the Celera facility, which occupies about 30,000 square feet of laboratory space and produces sequence data continuously at a rate of 175,000 total reads per day. The DNA-sequencing facility is supported by a high-performance computational facility (36).

The process for DNA sequencing was modular by design and automated. Intermodule sample backlogs allowed four principal modules to operate independently: (i) library transformation, plating, and colony picking; (ii) DNA template preparation; (iii) dideoxy sequencing reaction set-up and purification; and (iv) sequence determination with the ABI PRISM 3700 DNA Analyzer. Because the inputs and outputs of each module have been carefully matched and sample backlogs are continuously managed, sequencing has proceeded without a single day's interruption since the initiation of the Drosophila project in May 1999. The ABI 3700 is a fully automated capillary array sequencer and as such can be operated with a minimal amount of hands-on time, currently estimated at about 15 min per day. The capillary system also facilitates correct associations of sequencing traces with samples through the elimination of manual sample loading and lane-tracking errors associated with slab gels. About 65 production staff were hired and trained, and were rotated on a regular basis through the four production modules. A central laboratory information management system (LIMS) tracked all sample plates by unique bar code identifiers. The facility was supported by a quality control team that performed raw material and in-process testing and a quality assurance group with responsibilities including document control, validation, and auditing of the facility. Critical to the success of the scale-up was the validation of all software and instrumentation before implementation, and production-scale testing of any process changes.

1.2 Trace processing

1.3 Quality assessment and control

2 Genome Assembly Strategy and Characterization

Shotgun sequence assembly is a classic example of an inverse problem: given a set of reads randomly sampled from a target sequence, reconstruct the order and the position of those reads in the target. Genome assembly algorithms developed for Drosophila have now been extended to assemble the ~25-fold larger human genome. Celera assemblies consist of a set of contigs that are ordered and oriented into scaffolds that are then mapped to chromosomal locations by using known markers. The contigs consist of a collection of overlapping sequence reads that provide a consensus reconstruction for a contiguous interval of the genome. Mate pairs are a central component of the assembly strategy. They are used to produce scaffolds in which the size of gaps between consecutive contigs is known with reasonable precision. This is accomplished by observing that a pair of reads, one of which is in one contig, and the other of which is in another, implies an orientation and distance between the two contigs (Fig. 3). Finally, our assemblies did not incorporate all reads into the final set of reported scaffolds. This set of unincorporated reads is termed "chaff," and typically consisted of reads from within highly repetitive regions, data from other organisms introduced through various routes as found in many genome projects, and data of poor quality or with untrimmed vector.

2.1 Assembly data sets

The second data set was from the publicly funded Human Genome Project (PFP) and is primarily derived from BAC clones (30). The BAC data input to the assemblies came from a download of GenBank on 1 September 2000 (Table 2) totaling 4443.3 Mbp of sequence. The data for each BAC is deposited at one of four levels of completion. Phase 0 data are a set of generally unassembled sequencing reads from a very light shotgun of the BAC, typically less than 1×. Phase 1 data are unordered assemblies of contigs, which we call BAC contigs or bactigs. Phase 2 data are ordered assemblies of bactigs. Phase 3 data are complete BAC sequences. In the past 2 years the PFP has focused on a product of lower quality and completeness, but on a faster time-course, by concentrating on the production of Phase 1 data from a 3× to 4× light-shotgun of each BAC clone.

Table 2. GenBank data input into assembly. Center Statistics Completion phase sequence 0 1 and 2 3 Whitehead Institute/ Number of accession records 2,825 6,533 363 MIT Center for Number of contigs 243,786 138,023 363 Genome Research, Total base pairs 194,490,158 1,083,848,245 48,829,358 USA Total vector masked (bp) 1,553,597 875,618 2,202 Total contaminant masked (bp) 13,654,482 4,417,055 98,028 Average contig length (bp) 798 7,853 134,516 Washington University, Number of accession records 19 3,232 1,300 USA Number of contigs 2,127 61,812 1,300 Total base pairs 1,195,732 561,171,788 164,214,395 Total vector masked (bp) 21,604 270,942 8,287 Total contaminant masked (bp) 22,469 1,476,141 469,487 Average contig length (bp) 562 9,079 126,319 Baylor College of Number of accession records 0 1,626 363 Medicine, USA Number of contigs 0 44,861 363 Total base pairs 0 265,547,066 49,017,104 Total vector masked (bp) 0 218,769 4,960 Total contaminant masked (bp) 0 1,784,700 485,137 Average contig length (bp) 0 5,919 135,033 Production Sequencing Number of accession records 135 2,043 754 Facility, DOE Joint Number of contigs 7,052 34,938 754 Genome Institute, Total base pairs 8,680,214 294,249,631 60,975,328 USA Total vector masked (bp) 22,644 162,651 7,274 Total contaminant masked (bp) 665,818 4,642,372 118,387 Average contig length (bp) 1,231 8,422 80,867 The Institute of Physical Number of accession records 0 1,149 300 and Chemical Number of contigs 0 25,772 300 Research (RIKEN), Total base pairs 0 182,812,275 20,093,926 Japan Total vector masked (bp) 0 203,792 2,371 Total contaminant masked (bp) 0 308,426 27,781 Average contig length (bp) 0 7,093 66,978 Sanger Centre, UK Number of accession records 0 4,538 2,599 Number of contigs 0 74,324 2,599 Total base pairs 0 689,059,692 246,118,000 Total vector masked (bp) 0 427,326 25,054 Total contaminant masked (bp) 0 2,066,305 374,561 Average contig length (bp) 0 9,271 94,697 Others* Number of accession records 42 1,894 3,458 Number of contigs 5,978 29,898 3,458 Total base pairs 5,564,879 283,358,877 246,474,157 Total vector masked (bp) 57,448 279,477 32,136 Total contaminant masked (bp) 575,366 1,616,665 1,791,849 Average contig length (bp) 931 9,478 71,277 All centers combined Number of accession records 3,021 21,015 9,137 Number of contigs 258,943 409,628 9,137 Total base pairs 209,930,983 3,360,047,574 835,722,268 Total vector masked (bp) 1,655,293 2,438,575 82,284 Total contaminant masked (bp) 14,918,135 16,311,664 3,365,230 Average contig length (bp) 811 8,203 91,466 * Other centers contributing at least 0.1% of the sequence include: Chinese National Human Genome Center; Genomanalyse Gesellschaft fuer Biotechnologische Forschung mbH; Genome Therapeutics Corporation; GENOSCOPE; Chinese Academy of Sciences; Institute of Molecular Biotechnology; Keio University School of Medicine; Lawrence Livermore National Laboratory; Cold Spring Harbor Laboratory; Los Alamos National Laboratory; Max-Planck Institut fuer Molekulare, Genetik; Japan Science and Technology Corporation; Stanford University; The Institute for Genomic Research; The Institute of Physical and Chemical Research, Gene Bank; The University of Oklahoma; University of Texas Southwestern Medical Center, University of Washington. The 4,405,700,825 bases contributed by all centers were shredded into faux reads resulting in 2.96× coverage of the genome.

We screened the bactig sequences for contaminants by using the BLAST algorithm against three data sets: (i) vector sequences in Univec core (38), filtered for a 25-bp match at 98% sequence identity at the ends of the sequence and a 30-bp match internal to the sequence; (ii) the nonhuman portion of the High Throughput Genomic (HTG) Seqences division of GenBank (39), filtered at 200 bp at 98%; and (iii) the nonredundant nucleotide sequences from GenBank without primate and human virus entries, filtered at 200 bp at 98%. Whenever 25 bp or more of vector was found within 50 bp of the end of a contig, the tip up to the matching vector was excised. Under these criteria we removed 2.6 Mbp of possible contaminant and vector from the Phase 3 data, 61.0 Mbp from the Phase 1 and 2 data, and 16.1 Mbp from the Phase 0 data (Table 2). This left us with a total of 4363.7 Mbp of PFP sequence data 20% finished, 75% rough-draft (Phase 1 and 2), and 5% single sequencing reads (Phase 0). An additional 104,018 BAC end-sequence mate pairs were also downloaded and included in the data sets for both assembly processes (18).

2.2 Assembly strategies

For the whole-genome assembly, the PFP data was first disassembled or "shredded" into a synthetic shotgun data set of 550-bp reads that form a perfect 2× covering of the bactigs. This resulted in 16.05 million "faux" reads that were sufficient to cover the genome 2.96× because of redundancy in the BAC data set, without incorporating the biases inherent in the PFP assembly process. The combined data set of 43.32 million reads (8×), and all associated mate-pair information, were then subjected to our whole-genome assembly algorithm to produce a reconstruction of the genome. Neither the location of a BAC in the genome nor its assembly of bactigs was used in this process. Bactigs were shredded into reads because we found strong evidence that 2.13% of them were misassembled (40). Furthermore, BAC location information was ignored because some BACs were not correctly placed on the PFP physical map and because we found strong evidence that at least 2.2% of the BACs contained sequence data that were not part of the given BAC (41), possibly as a result of sample-tracking errors (see below). In short, we performed a true, ab initio whole-genome assembly in which we took the expedient of deriving additional sequence coverage, but not mate pairs, assembled bactigs, or genome locality, from some externally generated data.

In the compartmentalized shotgun assembly (CSA), Celera and PFP data were partitioned into the largest possible chromosomal segments or "components" that could be determined with confidence, and then shotgun assembly was applied to each partitioned subset wherein the bactig data were again shredded into faux reads to ensure an independent ab initio assembly of the component. By subsetting the data in this way, the overall computational effort was reduced and the effect of interchromosomal duplications was ameliorated. This also resulted in a reconstruction of the genome that was relatively independent of the whole-genome assembly results so that the two assemblies could be compared for consistency. The quality of the partitioning into components was crucial so that different genome regions were not mixed together. We constructed components from (i) the longest scaffolds of the sequence from each BAC and (ii) assembled scaffolds of data unique to Celera's data set. The BAC assemblies were obtained by a combining assembler that used the bactigs and the 5× Celera data mapped to those bactigs as input. This effort was undertaken as an interim step solely because the more accurate and complete the scaffold for a given sequence stretch, the more accurately one can tile these scaffolds into contiguous components on the basis of sequence overlap and mate-pair information. We further visually inspected and curated the scaffold tiling of the components to further increase its accuracy. For the final CSA assembly, all but the partitioning was ignored, and an independent, ab initio reconstruction of the sequence in each component was obtained by applying our whole-genome assembly algorithm to the partitioned, relevant Celera data and the shredded, faux reads of the partitioned, relevant bactig data.

2.3 Whole-genome assembly

The WGA assembler consists of a pipeline composed of five principal stages: Screener, Overlapper, Unitigger, Scaffolder, and Repeat Resolver, respectively. The Screener finds and marks all microsatellite repeats with less than a 6-bp element, and screens out all known interspersed repeat elements, including Alu, Line, and ribosomal DNA. Marked regions get searched for overlaps, whereas screened regions do not get searched, but can be part of an overlap that involves unscreened matching segments.

The Overlapper compares every read against every other read in search of complete end-to-end overlaps of at least 40 bp and with no more than 6% differences in the match. Because all data are scrupulously vector-trimmed, the Overlapper can insist on complete overlap matches. Computing the set of all overlaps took roughly 10,000 CPU hours with a suite of four-processor Alpha SMPs with 4 gigabytes of RAM. This took 4 to 5 days in elapsed time with 40 such machines operating in parallel.

Every overlap computed above is statistically a 1-in-10¹⁷ event and thus not a coincidental event. What makes assembly combinatorially difficult is that while many overlaps are actually sampled from overlapping regions of the genome, and thus imply that the sequence reads should be assembled together, even more overlaps are actually from two distinct copies of a low-copy repeated element not screened above, thus constituting an error if put together. We call the former "true overlaps" and the latter "repeat-induced overlaps." The assembler must avoid choosing repeat-induced overlaps, especially early in the process.

We achieve this objective in the Unitigger. We first find all assemblies of reads that appear to be uncontested with respect to all other reads. We call the contigs formed from these subassemblies unitigs (for quely assembled con). Formally, these unitigs are the uncontested interval subgraphs of the graph of all overlaps (42). Unfortunately, although empirically many of these assemblies are correct (and thus involve only true overlaps), some are in fact collections of reads from several copies of a repetitive element that have been overcollapsed into a single subassembly. However, the overcollapsed unitigs are easily identified because their average coverage depth is too high to be consistent with the overall level of sequence coverage. We developed a simple statistical discriminator that gives the logarithm of the odds ratio that a unitig is composed of unique DNA or of a repeat consisting of two or more copies. The discriminator, set to a sufficiently stringent threshold, identifies a subset of the unitigs that we are certain are correct. In addition, a second, less stringent threshold identifies a subset of remaining unitigs very likely to be correctly assembled, of which we select those that will consistently scaffold (see below), and thus are again almost certain to be correct. We call the union of these two sets U-unitigs. Empirically, we found from a 6× simulated shotgun of human chromosome 22 that we get U-unitigs covering 98% of the stretches of unique DNA that are >2 kbp long. We are further able to identify the boundary of the start of a repetitive element at the ends of a U-unitig and leverage this so that U-unitigs span more than 93% of all singly interspersed Alu elements and other 100-to 400-bp repetitive segments.

The result of running the Unitigger was thus a set of correctly assembled subcontigs covering an estimated 73.6% of the human genome. The Scaffolder then proceeded to use mate-pair information to link these together into scaffolds. When there are two or more mate pairs that imply that a given pair of U-unitigs are at a certain distance and orientation with respect to each other, the probability of this being wrong is again roughly 1 in 10¹⁰, assuming that mate pairs are false less than 2% of the time. Thus, one can with high confidence link together all U-unitigs that are linked by at least two 2- or 10-kbp mate pairs producing intermediate-sized scaffolds that are then recursively linked together by confirming 50-kbp mate pairs and BAC end sequences. This process yielded scaffolds that are on the order of megabase pairs in size with gaps between their contigs that generally correspond to repetitive elements and occasionally to small sequencing gaps. These scaffolds reconstruct the majority of the unique sequence within a genome.

For the Drosophila assembly, we engaged in a three-stage repeat resolution strategy where each stage was progressively more aggressive and thus more likely to make a mistake. For the human assembly, we continued to use the first "Rocks" substage where all unitigs with a good, but not definitive, discriminator score are placed in a scaffold gap. This was done with the condition that two or more mate pairs with one of their reads already in the scaffold unambiguously place the unitig in the given gap. We estimate the probability of inserting a unitig into an incorrect gap with this strategy to be less than 10⁷ based on a probabilistic analysis.

We revised the ensuing "Stones" substage of the human assembly, making it more like the mechanism suggested in our earlier work (43). For each gap, every read R that is placed in the gap by virtue of its mated pair M being in a contig of the scaffold and implying R's placement is collected. Celera's mate-pairing information is correct more than 99% of the time. Thus, almost every, but not all, of the reads in the set belong in the gap, and when a read does not belong it rarely agrees with the remainder of the reads. Therefore, we simply assemble this set of reads within the gap, eliminating any reads that conflict with the assembly. This operation proved much more reliable than the one it replaced for the Drosophila assembly; in the assembly of a simulated shotgun data set of human chromosome 22, all stones were placed correctly.

The final method of resolving gaps is to fill them with assembled BAC data that cover the gap. We call this external gap "walking." We did not include the very aggressive "Pebbles" substage described in our Drosophila work, which made enough mistakes so as to produce repeat reconstructions for long interspersed elements whose quality was only 99.62% correct. We decided that for the human genome it was philosophically better not to introduce a step that was certain to produce less than 99.99% accuracy. The cost was a somewhat larger number of gaps of somewhat larger size.

At the final stage of the assembly process, and also at several intermediate points, a consensus sequence of every contig is produced. Our algorithm is driven by the principle of maximum parsimony, with quality-value-weighted measures for evaluating each base. The net effect is a Bayesian estimate of the correct base to report at each position. Consensus generation uses Celera data whenever it is present. In the event that no Celera data cover a given region, the BAC data sequence is used.

A key element of achieving a WGA of the human genome was to parallelize the Overlapper and the central consensus sequence-constructing subroutines. In addition, memory was a real issue--a straightforward application of the software we had built for Drosophila would have required a computer with a 600-gigabyte RAM. By making the Overlapper and Unitigger incremental, we were able to achieve the same computation with a maximum of instantaneous usage of 28 gigabytes of RAM. Moreover, the incremental nature of the first three stages allowed us to continually update the state of this part of the computation as data were delivered and then perform a 7-day run to complete Scaffolding and Repeat Resolution whenever desired. For our assembly operations, the total compute infrastructure consists of 10 four-processor SMPs with 4 gigabytes of memory per cluster (Compaq's ES40, Regatta) and a 16-processor NUMA machine with 64 gigabytes of memory (Compaq's GS160, Wildfire). The total compute for a run of the assembler was roughly 20,000 CPU hours.

The assembly of Celera's data, together with the shredded bactig data, produced a set of scaffolds totaling 2.848 Gbp in span and consisting of 2.586 Gbp of sequence. The chaff, or set of reads not incorporated in the assembly, numbered 11.27 million (26%), which is consistent with our experience for Drosophila. More than 84% of the genome was covered by scaffolds >100 kbp long, and these averaged 91% sequence and 9% gaps with a total of 2.297 Gbp of sequence. There were a total of 93,857 gaps among the 1637 scaffolds >100 kbp. The average scaffold size was 1.5 Mbp, the average contig size was 24.06 kbp, and the average gap size was 2.43 kbp, where the distribution of each was essentially exponential. More than 50% of all gaps were less than 500 bp long, >62% of all gaps were less than 1 kbp long, and no gap was >100 kbp long. Similarly, more than 65% of the sequence is in contigs >30 kbp, more than 31% is in contigs >100 kbp, and the largest contig was 1.22 Mbp long. Table 3 gives detailed summary statistics for the structure of this assembly with a direct comparison to the compartmentalized shotgun assembly.

Table 3. Scaffold statistics for whole-genome and compartmentalized shotgun assemblies. Scaffold size All >30 kbp >100 kbp >500 kbp >1000 kbp Compartmentalized shotgun assembly No. of bp in scaffolds 2,905,568,203 2,748,892,430 2,700,489,906 2,489,357,260 2,248,689,128 (including intrascaffold gaps) No. of bp in contigs 2,653,979,733 2,524,251,302 2,491,538,372 2,320,648,201 2,106,521,902 No. of scaffolds 53,591 2,845 1,935 1,060 721 No. of contigs 170,033 112,207 107,199 93,138 82,009 No. of gaps 116,442 109,362 105,264 92,078 81,288 No. of gaps 1 kbp 72,091 69,175 67,289 59,915 53,354 Average scaffold size (bp) 54,217 966,219 1,395,602 2,348,450 3,118,848 Average contig size (bp) 15,609 22,496 23,242 24,916 25,686 Average intrascaffold gap size (bp) 2,161 2,054 1,985 1,832 1,749 Largest contig (bp) 1,988,321 1,988,321 1,988,321 1,988,321 1,988,321 % of total contigs 100 95 94 87 79 Whole-genome assembly No. of bp in scaffolds (including intrascaffold gaps) 2,847,890,390 2,574,792,618 2,525,334,447 2,328,535,466 2,140,943,032 No. of bp in contigs 2,586,634,108 2,334,343,339 2,297,678,935 2,143,002,184 1,983,305,432 No. of scaffolds 118,968 2,507 1,637 818 554 No. of contigs 221,036 99,189 95,494 84,641 76,285 No. of gaps 102,068 96,682 93,857 83,823 75,731 No. of gaps 1 kbp 62,356 60,343 59,156 54,079 49,592 Average scaffold size (bp) 23,938 1,027,041 1,542,660 2,846,620 3,864,518 Average contig size (bp) 11,702 23,534 24,061 25,319 25,999 Average intrascaffold gap size (bp) 2,560 2,487 2,426 2,213 2,082 Largest contig (bp) 1,224,073 1,224,073 1,224,073 1,224,073 1,224,073 % of total contigs 100 90 89 83 77

2.4 Compartmentalized shotgun assembly

The first phase of the CSA strategy was to separate Celera reads into those that matched the BAC contigs for a particular PFP BAC entry, and those that did not match any public data. Such matches must be guaranteed to properly place a Celera read, so all reads were first masked against a library of common repetitive elements, and only matches of at least 40 bp to unmasked portions of the read constituted a hit. Of Celera's 27.27 million reads, 20.76 million matched a bactig and another 0.62 million reads, which did not have any matches, were nonetheless identified as belonging in the region of the bactig's BAC because their mate matched the bactig. Of the remaining reads, 2.92 million were completely screened out and so could not be matched, but the other 2.97 million reads had unmasked sequence totaling 1.189 Gbp that were not found in the GenBank data set. Because the Celera data are 5.11× redundant, we estimate that 240 Mbp of unique Celera sequence is not in the GenBank data set.

In the next step of the CSA process, a combining assembler took the relevant 5× Celera reads and bactigs for a BAC entry, and produced an assembly of the combined data for that locale. These high-quality sequence reconstructions were a transient result whose utility was simply to provide more reliable information for the purposes of their tiling into sets of overlapping and adjacent scaffold sequences in the next step. In outline, the combining assembler first examines the set of matching Celera reads to determine if there are excessive pileups indicative of unscreened repetitive elements. Wherever these occur, reads in the repeat region whose mates have not been mapped to consistent positions are removed. Then all sets of mate pairs that consistently imply the same relative position of two bactigs are bundled into a link and weighted according to the number of mates in the bundle. A "greedy" strategy then attempts to order the bactigs by selecting bundles of mate-pairs in order of their weight. A selected mate-pair bundle can tie together two formative scaffolds. It is incorporated to form a single scaffold only if it is consistent with the majority of links between contigs of the scaffold. Once scaffolding is complete, gaps are filled by the "Stones" strategy described above for the WGA assembler.

The GenBank data for the Phase 1 and 2 BACs consisted of an average of 19.8 bactigs per BAC of average size 8099 bp. Application of the combining assembler resulted in individual Celera BAC assemblies being put together into an average of 1.83 scaffolds (median of 1 scaffold) consisting of an average of 8.57 contigs of average size 18,973 bp. In addition to defining order and orientation of the sequence fragments, there were 57% fewer gaps in the combined result. For Phase 0 data, the average GenBank entry consisted of 91.52 reads of average length 784 bp. Application of the combining assembler resulted in an average of 54.8 scaffolds consisting of an average of 58.1 contigs of average size 873 bp. Basically, some small amount of assembly took place, but not enough Celera data were matched to truly assemble the 0.5× to 1× data set represented by the typical Phase 0 BACs. The combining assembler was also applied to the Phase 3 BACs for SNP identification, confirmation of assembly, and localization of the Celera reads. The phase 0 data suggest that a combined whole-genome shotgun data set and 1× light-shotgun of BACs will not yield good assembly of BAC regions; at least 3× light-shotgun of each BAC is needed.

The 5.89 million Celera fragments not matching the GenBank data were assembled with our whole-genome assembler. The assembly resulted in a set of scaffolds totaling 442 Mbp in span and consisting of 326 Mbp of sequence. More than 20% of the scaffolds were >5 kbp long, and these averaged 63% sequence and 27% gaps with a total of 302 Mbp of sequence. All scaffolds >5 kbp were forwarded along with all scaffolds produced by the combining assembler to the subsequent tiling phase.

At this stage, we typically had one or two scaffolds for every BAC region constituting at least 95% of the relevant sequence, and a collection of disjoint Celera-unique scaffolds. The next step in developing the genome components was to determine the order and overlap tiling of these BAC and Celera-unique scaffolds across the genome. For this, we used Celera's 50-kbp mate-pairs information, and BAC-end pairs (18) and sequence tagged site (STS) markers (44) to provide long-range guidance and chromosome separation. Given the relatively manageable number of scaffolds, we chose not to produce this tiling in a fully automated manner, but to compute an initial tiling with a good heuristic and then use human curators to resolve discrepancies or missed join opportunities. To this end, we developed a graphical user interface that displayed the graph of tiling overlaps and the evidence for each. A human curator could then explore the implication of mapped STS data, dot-plots of sequence overlap, and a visual display of the mate-pair evidence supporting a given choice. The result of this process was a collection of "components," where each component was a tiled set of BAC and Celera-unique scaffolds that had been curator-approved. The process resulted in 3845 components with an estimated span of 2.922 Gbp.

In order to generate the final CSA, we assembled each component with the WGA algorithm. As was done in the WGA process, the bactig data were shredded into a synthetic 2× shotgun data set in order to give the assembler the freedom to independently assemble the data. By using faux reads rather than bactigs, the assembly algorithm could correct errors in the assembly of bactigs and remove chimeric content in a PFP data entry. Chimeric or contaminating sequence (from another part of the genome) would not be incorporated into the reassembly of the component because it did not belong there. In effect, the previous steps in the CSA process served only to bring together Celera fragments and PFP data relevant to a large contiguous segment of the genome, wherein we applied the assembler used for WGA to produce an ab initio assembly of the region.

WGA assembly of the components resulted in a set of scaffolds totaling 2.906 Gbp in span and consisting of 2.654 Gbp of sequence. The chaff, or set of reads not incorporated into the assembly, numbered 6.17 million, or 22%. More than 90.0% of the genome was covered by scaffolds spanning >100 kbp long, and these averaged 92.2% sequence and 7.8% gaps with a total of 2.492 Gbp of sequence. There were a total of 105,264 gaps among the 107,199 contigs that belong to the 1940 scaffolds spanning >100 kbp. The average scaffold size was 1.4 Mbp, the average contig size was 23.24 kbp, and the average gap size was 2.0 kbp where each distribution of sizes was exponential. As such, averages tend to be underrepresentative of the majority of the data. Figure 5 shows a histogram of the bases in scaffolds of various size ranges. Consider also that more than 49% of all gaps were <500 bp long, more than 62% of all gaps were <1 kbp, and all gaps are <100 kbp long. Similarly, more than 73% of the sequence is in contigs > 30 kbp, more than 49% is in contigs >100 kbp, and the largest contig was 1.99 Mbp long. Table 3 provides summary statistics for the structure of this assembly with a direct comparison to the WGA assembly.

2.5 Comparison of the WGA and CSA scaffolds

From this tabulation, we estimated the amount of unique sequence in each assembly in two ways. The first was to determine the number of bases of each assembly that were not covered by a matching segment in the other assembly. Some 82.5 Mbp of the WGA (3.95%) was not covered by the CSA, whereas 204.5 Mbp (8.26%) of the CSA was not covered by the WGA. This estimate did not require any consistency of the assemblies or any uniqueness of the matching segments. Thus, another analysis was conducted in which matches of less than 1 kbp between a pair of scaffolds were excluded unless they were confirmed by other matches having a consistent order and orientation. This gives some measure of consistent coverage: 1.982 Gbp (95.00%) of the WGA is covered by the CSA, and 2.169 Gbp (87.69%) of the CSA is covered by the WGA by this more stringent measure.

The comparison of WGA to CSA also permitted evaluation of scaffolds for structural inconsistencies. We looked for instances in which a large section of a scaffold from one assembly matched only one scaffold from the other assembly, but failed to match over the full length of the overlap implied by the matching segments. An initial set of candidates was identified automatically, and then each candidate was inspected by hand. From this process, we identified 31 instances in which the assemblies appear to disagree in a nonlocal fashion. These cases are being further evaluated to determine which assembly is in error and why.

In addition, we evaluated local inconsistencies of order or orientation. The following results exclude cases in which one contig in one assembly corresponds to more than one overlapping contig in the other assembly (as long as the order and orientation of the latter agrees with the positions they match in the former). Most of these small rearrangements involved segments on the order of hundreds of base pairs and rarely >1 kbp. We found a total of 295 kbp (0.012%) in the CSA assemblies that were locally inconsistent with the WGA assemblies, whereas 2.108 Mbp (0.11%) in the WGA assembly were inconsistent with the CSA assembly.

The CSA assembly was a few percentage points better in terms of coverage and slightly more consistent than the WGA, because it was in effect performing a few thousand shotgun assemblies of megabase-sized problems, whereas the WGA is performing a shotgun assembly of a gigabase-sized problem. When one considers the increase of two-and-a-half orders of magnitude in problem size, the information loss between the two is remarkably small. Because CSA was logistically easier to deliver and the better of the two results available at the time when downstream analyses needed to be begun, all subsequent analysis was performed on this assembly.

2.6 Mapping scaffolds to the genome

In order to determine the effectiveness of the fingerprint maps and GM99 for mapping scaffolds, we first examined the reliability of these maps by comparison with large scaffolds. Only 1% of the STS markers on the 10 largest scaffolds (those >9 Mbp) were mapped on a different chromosome on GM99. Two percent of the STS markers disagreed in position by more than five framework bins. However, for the fingerprint maps, a 2% chromosome discrepancy was observed, and on average 23.8% of BAC locations in the scaffold sequence disagreed with fingerprint map placement by more than five BACs. When further examining the source of discrepancy, it was found that most of the discrepancy came from 4 of the 10 scaffolds, indicating this there is variation in the quality of either the map or the scaffolds. All four scaffolds were assembled, as well as the other six, as judged by clone coverage analysis, and showed the same low discrepancy rate to GM99, and thus we concluded that the fingerprint map global order in these cases was not reliable. Smaller scaffolds had a higher discordance rate with GM99 (4.21% of STSs were discordant by more than five framework bins), but a lower discordance rate with the fingerprint maps (11% of BACs disagreed with fingerprint maps by more than five BACs). This observation agrees with the clone coverage analysis (46) that Celera scaffold construction was better supported by long-range mate pairs in larger scaffolds than in small scaffolds.

We created two orderings of Celera scaffolds on the basis of the markers (BAC or STS) on these maps. Where the order of scaffolds agreed between GM99 and the WashU BAC map, we had a high degree of confidence that that order was correct; these scaffolds were termed "anchor scaffolds." Only scaffolds with a low overall discrepancy rate with both maps were considered anchor scaffolds. Scaffolds in GM99 bins were allowed to permute in their order to match WashU ordering, provided they did not violate their framework orders. Orientation of individual scaffolds was determined by the presence of multiple mapped markers with consistent order. Scaffolds with only one marker have insufficient information to assign orientation. We found 70.1% of the genome in anchored scaffolds, more than 99% of which are also oriented (Table 4). Because GM99 is of lower resolution than the WashU map, a number of scaffolds without STS matches could be ordered relative to the anchored scaffolds because they included sequence from the same or adjacent BACs on the WashU map. On the other hand, because of occasional WashU global ordering discrepancies, a number of scaffolds determined to be "unmappable" on the WashU map could be ordered relative to the anchored scaffolds with GM99. These scaffolds were termed "ordered scaffolds." We found that 13.9% of the assembly could be ordered by these additional methods, and thus 84.0% of the genome was ordered unambiguously.

Table 4. Summary of scaffold mapping. Scaffolds were mapped to the genome with different levels of confidence (anchored scaffolds have the highest confidence; unmapped scaffolds have the lowest). Anchored scaffolds were consistently ordered by the WashU BAC map and GM99. Ordered scaffolds were consistently ordered by at least one of the following: the WashU BAC map, GM99, or component tiling path. Bounded scaffolds had order conflicts between at least two of the external maps, but their placements were adjacent to a neighboring anchored or ordered scaffold. Unmapped scaffolds had, at most, a chromosome assignment. The scaffold subcategories are given below each category. Mapped scaffold category Number Length (bp) % Total length Anchored 1,526 1,860,676,676 70 Oriented 1,246 1,852,088,645 70 Unoriented 280 8,588,031 0.3 Ordered 2,001 369,235,857 14 Oriented 839 329,633,166 12 Unoriented 1,162 39,602,691 2 Bounded 38,241 368,753,463 14 Oriented 7,453 274,536,424 10 Unoriented 30,788 94,217,039 4 Unmapped 11,823 55,313,737 2 Known 281 2,505,844 0.1 chromosome Unknown chromosome 11,542 52,807,893 2

Next, all scaffolds that could be placed, but not ordered, between anchors were assigned to the interval between the anchored scaffolds and were deemed to be "bounded" between them. For example, small scaffolds having STS hits from the same GeneMap bin or hitting the same BAC cannot be ordered relative to each other, but can be assigned a placement boundary relative to other anchored or ordered scaffolds. The remaining scaffolds either had no localization information, conflicting information, or could only be assigned to a generic chromosome location. Using the above approaches, ~98% of the genome was anchored, ordered, or bounded.

Finally, we assigned a location for each scaffold placed on the chromosome by spreading out the scaffolds per chromosome. We assumed that the remaining unmapped scaffolds, constituting 2% of the genome, were distributed evenly across the genome. By dividing the sum of unmapped scaffold lengths with the sum of the number of mapped scaffolds, we arrived at an estimate of interscaffold gap of 1483 bp. This gap was used to separate all the scaffolds on each chromosome and to assign an offset in the chromosome.

During the scaffold-mapping effort, we encountered many problems that resulted in additional quality assessment and validation analysis. At least 978 (3% of 33,173) BACs were believed to have sequence data from more than one location in the genome (47). This is consistent with the bactig chimerism analysis reported above in the Assembly Strategies section. These BACs could not be assigned to unique positions within the CSA assembly and thus could not be used for ordering scaffolds. Likewise, it was not always possible to assign STSs to unique locations in the assembly because of genome duplications, repetitive elements, and pseudogenes.

Because of the time required for an exhaustive search for a perfect overlap, CSA generated 21,607 intrascaffold gaps where the mate-pair data suggested that the contigs should overlap, but no overlap was found. These gaps were defined as a fixed 50 bp in length and make up 18.6% of the total 116,442 gaps in the CSA assembly.

We chose not to use the order of exons implied in cDNA or EST data as a way of ordering scaffolds. The rationale for not using this data was that doing so would have biased certain regions of the assembly by rearranging scaffolds to fit the transcript data and made validation of both the assembly and gene definition processes more difficult.

2.7 Assembly and validation analysis

Completeness. Completeness is defined as the percentage of the euchromatic sequence represented in the assembly. This cannot be known with absolute certainty until the euchromatin sequence has been completed. However, it is possible to estimate completeness on the basis of (i) the estimated sizes of intrascaffold gaps; (ii) coverage of the two published chromosomes, 21 and 22 (48, 49); and (iii) analysis of the percentage of an independent set of random sequences (STS markers) contained in the assembly. The whole-genome libraries contain heterochromatic sequence and, although no attempt has been made to assemble it, there may be instances of unique sequence embedded in regions of heterochromatin as were observed in Drosophila (50, 51).

The sequences of human chromosomes 21 and 22 have been completed to high quality and published (48, 49). Although this sequence served as input to the assembler, the finished sequence was shredded into a shotgun data set so that the assembler had the opportunity to assemble it differently from the original sequence in the case of structural polymorphisms or assembly errors in the BAC data. In particular, the assembler must be able to resolve repetitive elements at the scale of components (generally multimegabase in size), and so this comparison reveals the level to which the assembler resolves repeats. In certain areas, the assembly structure differs from the published versions of chromosomes 21 and 22 (see below). The consequence of the flexibility to assemble "finished" sequence differently on the basis of Celera data resulted in an assembly with more segments than the chromosome 21 and 22 sequences. We examined the reasons why there are more gaps in the Celera sequence than in chromosomes 21 and 22 and expect that they may be typical of gaps in other regions of the genome. In the Celera assembly, there are 25 scaffolds, each containing at least 10 kb of sequence, that collectively span 94.3% of chromosome 21. Sixty-two scaffolds span 95.7% of chromosome 22. The total length of the gaps remaining in the Celera assembly for these two chromosomes is 3.4 Mbp. These gap sequences were analyzed by RepeatMasker and by searching against the entire genome assembly (52). About 50% of the gap sequence consisted of common repetitive elements identified by RepeatMasker; more than half of the remainder was lower copy number repeat elements.

A more global way of assessing completeness is to measure the content of an independent set of sequence data in the assembly. We compared 48,938 STS markers from Genemap99 (51) to the scaffolds. Because these markers were not used in the assembly processes, they provided a truly independent measure of completeness. ePCR (53) and BLAST (54) were used to locate STSs on the assembled genome. We found 44,524 (91%) of the STSs in the mapped genome. An additional 2648 markers (5.4%) were found by searching the unassembled data or "chaff." We identified 1283 STS markers (2.6%) not found in either Celera sequence or BAC data as of September 2000, raising the possibility that these markers may not be of human origin. If that were the case, the Celera assembled sequence would represent 93.4% of the human genome and the unassembled data 5.5%, for a total of 98.9% coverage. Similarly, we compared CSA against 36,678 TNG radiation hybrid markers (55a) using the same method. We found that 32,371 markers (88%) were located in the mapped CSA scaffolds, with 2055 markers (5.6%) found in the remainder. This gave a 94% coverage of the genome through another genome-wide survey.

Correctness. Correctness is defined as the structural and sequence accuracy of the assembly. Because the source sequences for the Celera data and the GenBank data are from different individuals, we could not directly compare the consensus sequence of the assembly against other finished sequence for determining sequencing accuracy at the nucleotide level, although this has been done for identifying polymorphisms as described in Section 6. The accuracy of the consensus sequence is at least 99.96% on the basis of a statistical estimate derived from the quality values of the underlying reads.

The structural consistency of the assembly can be measured by mate-pair analysis. In a correct assembly, every mated pair of sequencing reads should be located on the consensus sequence with the correct separation and orientation between the pairs. A pair is termed "valid" when the reads are in the correct orientation and the distance between them is within the mean ± 3 standard deviations of the distribution of insert sizes of the library from which the pair was sampled. A pair is termed "misoriented" when the reads are not correctly oriented, and is termed "misseparated" when the distance between the reads is not in the correct range but the reads are correctly oriented. The mean ± the standard deviation of each library used by the assembler was determined as described above. To validate these, we examined all reads mapped to the finished sequence of chromosome 21 (48) and determined how many incorrect mate pairs there were as a result of laboratory tracking errors and chimerism (two different segments of the genome cloned into the same plasmid), and how tight the distribution of insert sizes was for those that were correct (Table 5). The standard deviations for all Celera libraries were quite small, less than 15% of the insert length, with the exception of a few 50-kbp libraries. The 2- and 10-kbp libraries contained less than 2% invalid mate pairs, whereas the 50-kbp libraries were somewhat higher (~10%). Thus, although the mate-pair information was not perfect, its accuracy was such that measuring valid, misoriented, and misseparated pairs with respect to a given assembly was deemed to be a reliable instrument for validation purposes, especially when several mate pairs confirm or deny an ordering.

Table 5. Mate-pair validation. Celera fragment sequences were mapped to the published sequence of chromosome 21. Each mate pair uniquely mapped was evaluated for correct orientation and placement (number of mate pairs tested). If the two mates had incorrect relative orientation or placement, they were considered invalid (number of invalid mate pairs). Library type Library no. Chromosome 21 Genome Mean insert size (bp) SD (bp) SD/mean (%) No. of mate pairs tested No. of invalid mate pairs % invalid Mean insert size (bp) SD (bp) SD/ mean (%) 2 kbp 1 2,081 106 5.1 3,642  38 1.0 2,082 90 4.3 2 1,913 152 7.9 28,029 413 1.5 1,923 118 6.1 3 2,166 175 8.1 4,405  57 1.3 2,162 158 7.3 10 kbp 4 11,385 851 7.5 4,319  80 1.9 11,370 696 6.1 5 14,523 1,875 12.9 7,355 156 2.1 14,142 1,402 9.9 6 9,635 1,035 10.7 5,573 109 2.0 9,606 934 9.7 7 10,223 928 9.1 34,079 399 1.2 10,190 777 7.6 50 kbp 8 64,888 2,747 4.2 16   1 6.3 65,500 5,504 8.4 9 53,410 5,834 10.9 914 170 18.6 53,311 5,546 10.4 10 52,034 7,312 14.1 5,871 569 9.7 51,498 6,588 12.8 11 52,282 7,454 14.3 2,629 213 8.1 52,282 7,454 14.3 12 46,616 7,378 15.8 2,153 215 10.0 45,418 9,068 20.0 13 55,788 10,099 18.1 2,244 249 11.1 53,062 10,893 20.5 14 39,894 5,019 12.6 199   7 3.5 36,838 9,988 27.1 BES 15 48,931 9,813 20.1 144  10 6.9 47,845 4,774 10.0 16 48,130 4,232 8.8 195  14 7.2 47,924 4,581 9.6 17 106,027 27,778 26.2 330  16 4.8 152,000 26,600 17.5 18 160,575 54,973 34.2 155   8 5.2 161,750 27,000 16.7 19 164,155 19,453 11.9 642  44 6.9 176,500 19,500 11.05 Sum 102,894 2,768 2.7 (mean = 2.7)

The clone coverage of the genome was 39×, meaning that any given base pair was, on average, contained in 39 clones or, equivalently, spanned by 39 mate-paired reads. Areas of low clone coverage or areas with a high proportion of invalid mate pairs would indicate potential assembly problems. We computed the coverage of each base in the assembly by valid mate pairs (Table 6). In summary, for scaffolds >30 kbp in length, less than 1% of the Celera assembly was in regions of less than 3× clone coverage. Thus, more than 99% of the assembly, including order and orientation, is strongly supported by this measure alone.

Table 6. Genome-wide mate pair analysis of compartmentalized shotgun (CSA) and PFP assemblies.* Genome library CSA PFP % valid % mis-oriented % mis-separated % valid % mis-oriented % mis-separated 2 kbp 98.5 0.6 1.0 95.7 2.0 2.3 10 kbp 96.7 1.0 2.3 81.9 9.6 8.6 50 kbp 93.9 4.5 1.5 64.2 22.3 13.5 BES 94.1 2.1 3.8 62.0 19.3 18.8 Mean 97.4 1.0 1.6 87.3 6.8 5.9 * Data for individual chromosomes can be found in Web fig. 3 on Science Online at www.sciencemag.org/cgi/content/full/291/5507/1304/DC1. Mates are misseparated if their distance is >3 SD from the mean library size.

We examined the locations and number of all misoriented and misseparated mates. In addition to doing this analysis on the CSA assembly (as of 1 October 2000), we also performed a study of the PFP assembly as of 5 September 2000 (30, 55b). In this latter case, Celera mate pairs had to be mapped to the PFP assembly. To avoid mapping errors due to high-fidelity repeats, the only pairs mapped were those for which both reads matched at only one location with less than 6% differences. A threshold was set such that sets of five or more simultaneously invalid mate pairs indicated a potential breakpoint, where the construction of the two assemblies differed. The graphic comparison of the CSA chromosome 21 assembly with the published sequence (Fig. 6A) serves as a validation of this methodology. Blue tick marks in the panels indicate breakpoints. There were a similar (small) number of breakpoints on both chromosome sequences. The exception was 12 sets of scaffolds in the Celera assembly (a total of 3% of the chromosome length in 212 single-contig scaffolds) that were mapped to the wrong positions because they were too small to be mapped reliably. Figures 6 and 7 and Table 6 illustrate the mate-pair differences and breakpoints between the two assemblies. There was a higher percentage of misoriented and misseparated mate pairs in the large-insert libraries (50 kbp and BAC ends) than in the small-insert libraries in both assemblies (Table 6). The large-insert libraries are more likely to identify discrepancies simply because they span a larger segment of the genome. The graphic comparison between the two assemblies for chromosome 8 (Fig. 6, B and C) shows that there are many more breakpoints for the PFP assembly than for the Celera assembly. Figure 7 shows the breakpoint map (blue tick marks) for both assemblies of each chromosome in a side-by-side fashion. The order and orientation of Celera's assembly shows substantially fewer breakpoints except on the two finished chromosomes. Figure 7 also depicts large gaps (>10 kbp) in both assemblies as red tick marks. In the CSA assembly, the size of all gaps have been estimated on the basis of the mate-pair data. Breakpoints can be caused by structural polymorphisms, because the two assemblies were derived from different human genomes. They also reflect the unfinished nature of both genome assemblies.

3 Gene Prediction and Annotation

3.1 Automated gene annotation

The number of protein-coding genes in mammals has been controversial from the outset. Initial estimates based on reassociation data placed it between 30,000 to 40,000, whereas later estimates from the brain were >100,000 (56). More recent data from both the corporate and public sectors, based on extrapolations from EST, CpG island, and transcript density-based extrapolations, have not reduced this variance. The highest recent number of 142,634 genes emanates from a report from Incyte Pharmaceuticals, and is based on a combination of EST data and the association of ESTs with CpG islands (57). In stark contrast are three quite different, and much lower estimates: one of ~35,000 genes derived with genome-wide EST data and sampling procedures in conjunction with chromosome 22 data (58); another of 28,000 to 34,000 genes derived with a comparative methodology involving sequence conservation between humans and the puffer fish Tetraodon nigroviridis (59); and a figure of 35,000 genes, which was derived simply by extrapolating from the density of 770 known and predicted genes in the 67 Mbp of chromosomes 21 and 22, to the approximately 3-Gbp euchromatic genome.

The problem of computational identification of transcriptional units in genomic DNA sequence can be divided into two phases. The first is to partition the sequence into segments that are likely to correspond to individual genes. This is not trivial and is a weakness of most de novo gene-finding algorithms. It is also critical to determining the number of genes in the human gene inventory. The second challenge is to construct a gene model that reflects the probable structure of the transcript(s) encoded in the region. This can be done with reasonable accuracy when a full-length cDNA has been sequenced or a highly homologous protein sequence is known. De novo gene prediction, although less accurate, is the only way to find genes that are not represented by homologous proteins or ESTs. The following section describes the methods we have developed to address these problems for the prediction of protein-coding genes.

We have developed a rule-based expert system, called Otto, to identify and characterize genes in the human genome (60). Otto attempts to simulate in software the process that a human annotator uses to identify a gene and refine its structure. In the process of annotating a region of the genome, a human curator examines the evidence provided by the computational pipeline (described below) and examines how various types of evidence relate to one another. A curator puts different levels of confidence in different types of evidence and looks for certain patterns of evidence to support gene annotation. For example, a curator may examine homology to a number of ESTs and evaluate whether or not they can be connected into a longer, virtual mRNA. The curator would also evaluate the strength of the similarity and the contiguity of the match, in essence asking whether any ESTs cross splice-junctions and whether the edges of putative exons have consensus splice sites. This kind of manual annotation process was used to annotate the Drosophila genome.

The Otto system can promote observed evidence to a gene annotation in one of two ways. First, if the evidence includes a high-quality match to the sequence of a known gene [here defined as a human gene represented in a curated subset of the RefSeq database (61)], then Otto can promote this to a gene annotation. In the second method, Otto evaluates a broad spectrum of evidence and determines if this evidence is adequate to support promotion to a gene annotation. These processes are described below.

Initially, gene boundaries are predicted on the basis of examination of sets of overlapping protein and EST matches generated by a computational pipeline (62). This pipeline searches the scaffold sequences against protein, EST, and genome-sequence databases to define regions of sequence similarity and runs three de novo gene-prediction programs.

To identify likely gene boundaries, regions of the genome were partitioned by Otto on the basis of sequence matches identified by BLAST. Each of the database sequences matched in the region under analysis was compared by an algorithm that takes into account both coordinates of the matching sequence, as well as the sequence type (e.g., protein, EST, and so forth). The results were used to group the matches into bins of related sequences that may define a gene and identify gene boundaries. During this process, multiple hits to the same region were collapsed to a coherent set of data by tracking the coverage of a region. For example, if a group of bases was represented by multiple overlapping ESTs, the union of these regions matched by the set of ESTs on the scaffold was marked as being supported by EST evidence. This resulted in a series of "gene bins," each of which was believed to contain a single gene. One weakness of this initial implementation of the algorithm was in predicting gene boundaries in regions of tandemly duplicated genes. Gene clusters frequently resulted in homologous neighboring genes being joined together, resulting in an annotation that artificially concatenated these gene models.

Next, known genes (those with exact matches of a full-length cDNA sequence to the genome) were identified, and the region corresponding to the cDNA was annotated as a predicted transcript. A subset of the curated human gene set RefSeq from the National Center for Biotechnology Information (NCBI) was included as a data set searched in the computational pipeline. If a RefSeq transcript matched the genome assembly for at least 50% of its length at >92% identity, then the SIM4 (63) alignment of the RefSeq transcript to the region of the genome under analysis was promoted to the status of an Otto annotation. Because the genome sequence has gaps and sequence errors such as frameshifts, it was not always possible to predict a transcript that agrees precisely with the experimentally determined cDNA sequence. A total of 6538 genes in our inventory were identified and transcripts predicted in this way.

Regions that have a substantial amount of sequence similarity, but do not match known genes, were analyzed by that part of the Otto system that uses the sequence similarity information to predict a transcript. Here, Otto evaluates evidence generated by the computational pipeline, corresponding to conservation between mouse and human genomic DNA, similarity to human transcripts (ESTs and cDNAs), similarity to rodent transcripts (ESTs and cDNAs), and similarity of the translation of human genomic DNA to known proteins to predict potential genes in the human genome. The sequence from the region of genomic DNA contained in a gene bin was extracted, and the subsequences supported by any homology evidence were marked (plus 100 bases flanking these regions). The other bases in the region, those not covered by any homology evidence, were replaced by N's. This sequence segment, with high confidence regions represented by the consensus genomic sequence and the remainder represented by N's, was then evaluated by Genscan to see if a consistent gene model could be generated. This procedure simplified the gene-prediction task by first establishing the boundary for the gene (not a strength of most gene-finding algorithms), and by eliminating regions with no supporting evidence. If Genscan returned a plausible gene model, it was further evaluated before being promoted to an "Otto" annotation. The final Genscan predictions were often quite different from the prediction that Genscan returned on the same region of native genomic sequence. A weakness of using Genscan to refine the gene model is the loss of valid, small exons from the final annotation.

The next step in defining gene structures based on sequence similarity was to compare each predicted transcript with the homology-based evidence that was used in previous steps to evaluate the depth of evidence for each exon in the prediction. Internal exons were considered to be supported if they were covered by homology evidence to within ±10 bases of their edges. For first and last exons, the internal edge was required to be within 10 bases, but the external edge was allowed greater latitude to allow for 5' and 3' untranslated regions (UTRs). To be retained, a prediction for a multi-exon gene must have evidence such that the total number of "hits," as defined above, divided by the number of exons in the prediction must be >0.66 or must correspond to a RefSeq sequence. A single-exon gene must be covered by at least three supporting hits (±10 bases on each side), and these must cover the complete predicted open reading frame. For a single-exon gene, we also required that the Genscan prediction include both a start and a stop codon. Gene models that did not meet these criteria were disregarded, and those that passed were promoted to Otto predictions. Homology-based Otto predictions do not contain 3' and 5' untranslated sequence. Although three de novo gene-finding programs [GRAIL, Genscan, and FgenesH (63)] were run as part of the computational analysis, the results of these programs were not directly used in making the Otto predictions. Otto predicted 11,226 additional genes by means of sequence similarity.

3.2 Otto validation

Table 7. Sensitivity and specificity of Otto and Genscan. Sensitivity and specificity were calculated by first aligning the prediction to the published RefSeq transcript, tallying the number (N) of uniquely aligned RefSeq bases. Sensitivity is the ratio of N to the length of the published RefSeq transcript. Specificity is the ratio of N to the length of the prediction. All differences are significant (Tukey HSD; P < 0.001). Method Sensitivity Specificity Otto (RefSeq only)* 0.939 0.973 Otto (homology) 0.604 0.884 Genscan 0.501 0.633 * Refers to those annotations produced by Otto using only the Sim4-polished RefSeq alignment rather than an evidence-based Genscan prediction. Refers to those annotations produced by supplying all available evidence to Genscan.

Because Otto uses an evidence-based approach to reconstruct genes, the absence of experimental evidence for intervening exons may inadvertantly result in a set of exons that cannot be spliced together to give rise to a transcript. In such cases, Otto may "split genes" when in fact all the evidence should be combined into a single transcript. We also examined the tendency of these methods to incorrectly split gene predictions. These trends are shown in Fig. 8. Both RefSeq and homology-based predictions by Otto split known genes into fewer segments than Genscan alone.

3.3 Gene number

Table 8. Numbers of exons and transcripts supported by various types of evidence for Otto and de novo gene prediction methods. Highlighted cells indicate the gene sets analyzed in this paper (boldface, set of genes selected for protein analysis; italic, total set of accepted de novo predictions). Total Types of evidence No. of lines of evidence* Mouse Rodent Protein Human  1  2  3  4 Otto Number of transcripts 17,969 17,065 14,881 15,477 16,374 17,968 17,501 15,877 12,451 Number of exons 141,218 111,174 89,569 108,431 118,869 140,710  127,955 99,574 59,804 De novo Number of transcripts 58,032 14,463 5,094 8,043 9,220  21,350  8,619 4,947 1,904 Number of exons 319,935 48,594 19,344 26,264 40,104  79,148  31,130 17,508 6,520 No. of exons per Otto 7.84 5.77 6.01 6.99 7.24 7.81 7.19 6.00 4.28 transcript De novo 5.53 3.17 3.80 3.27 4.36  3.7 3.56 3.42 3.16 * Four kinds of evidence (conservation in 3× mouse genomic DNA, similarity to human EST or cDNA, similarity to rodent EST or cDNA, and similarity to known proteins) were considered to support gene predictions from the different methods. The use of evidence is quite liberal, requiring only a partial match to a single exon of predicted transcript. This number includes alternative splice forms of the 17,764 genes mentioned elsewhere in the text.

The sum of this number (21,350) and the number of Otto annotations (17,764), 39,114, is near the upper limit for the human gene complement. As seen in Table 8, if the requirement for other supporting evidence is made more stringent, this number drops rapidly so that demanding two types of evidence reduces the total gene number to 26,383 and demanding three types reduces it to ~23,000. Requiring that a prediction be supported by all four categories of evidence is too stringent because it would eliminate genes that encode novel proteins (members of currently undescribed protein families). No correction for pseudogenes has been made at this point in the analysis.

In a further attempt to identify genes that were not found by the autoannotation process or any of the de novo gene finders, we examined regions outside of gene predictions that were similar to the EST sequence, and where the EST matched the genomic sequence across a splice junction. After correcting for potential 3' UTRs of predicted genes, about 2500 such regions remained. Addition of a requirement for at least one of the following evidence types--homology to mouse genomic sequence fragments, rodent ESTs, or cDNAs--or similarity to a known protein reduced this number to 1010. Adding this to the numbers from the previous paragraph would give us estimates of about 40,000, 27,000, and 24,000 potential genes in the human genome, depending on the stringency of evidence considered. Table 8 illustrates the number of genes and presents the degree of confidence based on the supporting evidence. Transcripts encoded by a set of 26,383 genes were assembled for further analysis. This set includes the 6538 genes predicted by Otto on the basis of matches to known genes, 11,226 transcripts predicted by Otto based on homology evidence, and 8619 from the subset of transcripts from de novo gene-prediction programs that have two types of supporting evidence. The 26,383 genes are illustrated along chromosome diagrams in Fig. 1. These are a very preliminary set of annotations and are subject to all the limitations of an automated process. Considerable refinement is still necessary to improve the accuracy of these transcript predictions. All the predictions and descriptions of genes and the associated evidence that we present are the product of completely computational processes, not expert curation. We have attempted to enumerate the genes in the human genome in such a way that we have different levels of confidence based on the amount of supporting evidence: known genes, genes with good protein or EST homology evidence, and de novo gene predictions confirmed by modest homology evidence.

3.4 Features of human gene transcripts

The set of transcripts promoted to gene annotations varies in a number of ways. As can be seen from Table 8 and Fig. 9, transcripts predicted by Otto tend to be longer, having on average about 7.8 exons, whereas those promoted from gene-prediction programs average about 3.7 exons. The largest number of exons that we have identified in a transcript is 234 in the titin mRNA. Table 8 compares the amounts of evidence that support the Otto and other predicted transcripts. For example, one can see that a typical Otto transcript has 6.99 of its 7.81 exons supported by protein homology evidence. As would be expected, the Otto transcripts generally have more support than do transcripts predicted by the de novo methods.

4 Genome Structure

4.1 Cytogenetic maps

The remaining ~80% of the genome, the euchromatic component, is divisible into G-, R-, and T-bands (67). These cytogenetic bands have been presumed to differ in their nucleotide composition and gene density, although we have been unable to determine precise band boundaries at the molecular level. T-bands are the most G+C- and gene-rich, and G-bands are G+C-poor (68). Bernardi has also offered a description of the euchromatin at the molecular level as long stretches of DNA of differing base composition, termed isochores (denoted L, H1, H2, and H3), which are >300 kbp in length (69). Bernardi defined the L (light) isochores as G+C-poor (<43%), whereas the H (heavy) isochores fall into three G+C-rich classes representing 24, 8, and 5% of the genome. Gene concentration has been claimed to be very low in the L isochores and 20-fold more enriched in the H2 and H3 isochores (70). By examining contiguous 50-kbp windows of G+C content across the assembly, we found that regions of G+C content >48% (H3 isochores) averaged 273.9 kbp in length, those with G+C content between 43 and 48% (H1+H2 isochores) averaged 202.8 kbp in length, and the average span of regions with <43% (L isochores) was 1078.6 kbp. The correlation between G+C content and gene density was also examined in 50-kbp windows along the assembled sequence (Table 9 and Figs. 10 and 11). We found that the density of genes was greater in regions of high G+C than in regions of low G+C content, as expected. However, the correlation between G+C content and gene density was not as skewed as previously predicted (69). A higher proportion of genes were located in the G+C-poor regions than had been expected.

Table 9. Characteristics of G+C in isochores. Isochore G+C (%) Fraction of genome Fraction of genes Predicted* Observed Predicted* Observed H3 >48 5 9.5 37 24.8 H1/H2 43-48 25 21.2 32 26.6 L <43 67 69.2 31 48.5 * The predictions were based on Bernardi's definitions (70) of the isochore structure of the human genome.

Chromosomes 17, 19, and 22, which have a disproportionate number of H3-containing bands, had the highest gene density (Table 10). Conversely, of the chromosomes that we found to have the lowest gene density, X, 4, 18, 13, and Y, also have the fewest H3 bands. Chromosome 15, which also has few H3 bands, did not have a particularly low gene density in our analysis. In addition, chromosome 8, which we found to have a low gene density, does not appear to be unusual in its H3 banding.

Table 10. Features of the chromosomes. De novo/any refers to the union of de novo predictions that do not overlap Otto predictions and have at least one other type of supporting evidence; de novo/2x refers to the union of de novo predictions that do not overlap Otto predictions and have at least two types of evidence. Deserts are regions of sequence with no annotated genes. Chr. Sequence coverage (CS assembly) Base composition Gene prediction* Gene density (genes/Mbp) Size (Mbp) No. of scaf- folds Largest scaf- fold (Mbp) No. of scaf- folds >500 kbp Se- quence covered by scaf- folds >500 kbp % of total se- quence in scaf-folds >500 kbp % repeat % GC No of CpG islands Otto De novo/any De novo/2× Total (Otto + de novo/any) Total (Otto + de novo/2×) Se- quence in deserts >500/kbp Se- quence in deserts >1 Mbp Otto De novo/ any De novo/2× Otto + de novo/ any Otto + de novo/ 2×  1 220 2,549 11 82 192 88 37 42 2,335 1,743 1,710 710 3,453 2,453 29 6 8 8 3 16 11  2 240 3,263 13 78 217 91 36 40 1,703 1,183 1,771 633 2,954 1,816 55 19 5 7 2 12 7  3 200 3,532 7 78 173 87 37 40 1,271 1,013 1,414 598 2,427 1,611 50 12 5 7 3 12 8  4 186 2,180 10 70 169 91 37 38 1,081 696 1,165 449 1,861 1,145 55 18 4 6 2 10 6  5 182 3,231 11 63 163 89 37 40 1,302 892 1,244 474 2,136 1,366 46 15 5 7 2 11 7  6 172 1,713 13 58 160 93 37 40 1,384 943 1,314 524 2,257 1,467 38 9 6 7 3 13 8  7 146 1,326 14 53 130 89 38 40 1,406 759 1,072 460 1,831 1,219 26 12 5 7 3 12 8  8 146 1,772 11 54 135 92 36 40 948 583 977 357 1,560 940 33 6 4 7 2 11 6  9 113 1,616 8 40 101 89 38 41 1,315 689 848 329 1,537 1,018 22 9 6 7 3 13 8 10 130 2,005 9 55 116 89 36 42 1,087 685 968 342 1,653 1,027 21 8 5 7 2 12 7 11 132 2,814 9 44 116 88 39 42 1,461 1,051 1,134 535 2,185 1,586 27 9 8 8 4 16 12 12 134 2,614 8 51 117 87 38 41 1,131 925 936 417 1,861 1,342 24 9 7 7 3 14 10 13 99 1,038 13 34 91 91 36 38 644 341 691 241 1,032 582 31 16 4 7 2 10 5 14 87 576 11 16 83 95 40 41 913 583 700 290 1,283 873 34 20 7 8 3 14 10 15 80 1,747 8 31 70 87 37 42 722 558 640 246 1,198 804 8 1 7 8 3 15 10 16 75 1,520 8 27 62 82 40 44 1,533 748 673 247 1,421 995 13 3 10 9 3 19 12 17 78 1,683 6 40 61 78 39 45 1,489 897 648 313 1,545 1,210 15 6 12 8 4 19 15 18 79 1,333 13 18 72 92 36 40 510 283 543 189 826 472 21 10 4 7 2 10 6 19 58 2,282 3 31 38 67 57 49 2,804 1,141 534 268 1,675 1,409 3 0 20 9 4 29 23 20 61 580 14 17 58 94 41 44 997 517 469 180 986 697 7 1 8 7 3 16 11 21 33 358 10 6 32 96 38 41 519 184 265 102 449 286 15 9 6 8 3 13 8 22 36 333 11 12 32 88 44 48 1,173 494 341 147 835 641 3 0 14 9 4 23 17 X 128 1,346 4 91 93 73 46 39 726 605 860 387 1,465 992 29 8 5 6 3 11 7 Y 19 638 2 10 12 65 50 39 65 55 155 49 210 104 4 2 3 8 2 11 5 U* 75 11,542 1 479 196 278 132 474 328 Total 2907 53,591 1,059 2,490 28,519 17,764 21,350 8,619 39,114 26,383 606 208 Avg. 116 2,144 9 44 104 87 40 41 1,160 714 812 333 1,526 1,047 25 9 7 7 3 14 9 * Chromosomal assignment unknown.

How valid is Ohno's postulate (71) that mammalian genomes consist of oases of genes in otherwise essentially empty deserts? It appears that the human genome does indeed contain deserts, or large, gene-poor regions. If we define a desert as a region >500 kbp without a gene, then we see that 605 Mbp, or about 20% of the genome, is in deserts. These are not uniformly distributed over the various chromosomes. Gene-rich chromosomes 17, 19, and 22 have only about 12% of their collective 171 Mbp in deserts, whereas gene-poor chromosomes 4, 13, 18, and X have 27.5% of their 492 Mbp in deserts (Table 11). The apparent lack of predicted genes in these regions does not necessarily imply that they are devoid of biological function.

Table 11. Genome overview. Size of the genome (including gaps) 2.91 Gbp Size of the genome (excluding gaps) 2.66 Gbp Longest contig 1.99 Mbp Longest scaffold 14.4 Mbp Percent of A+T in the genome 54 Percent of G+C in the genome 38 Percent of undetermined bases in the genome 9 Most GC-rich 50 kb Chr. 2 (66%) Least GC-rich 50 kb Chr. X (25%) Percent of genome classified as repeats 35 Number of annotated genes 26,383 Percent of annotated genes with unknown function 42 Number of genes (hypothetical and annotated) 39,114 Percent of hypothetical and annotated genes with unknown function 59 Gene with the most exons Titin (234 exons) Average gene size 27 kbp Most gene-rich chromosome Chr. 19 (23 genes/Mb) Least gene-rich chromosomes Chr. 13 (5 genes/Mb), Chr. Y (5 genes/Mb) Total size of gene deserts (>500 kb with no annotated genes) 605 Mbp Percent of base pairs spanned by genes 25.5 to 37.8* Percent of base pairs spanned by exons 1.1 to 1.4* Percent of base pairs spanned by introns 24.4 to 36.4* Percent of base pairs in intergenic DNA 74.5 to 63.6* Chromosome with highest proportion of DNA in annotated exons Chr. 19 (9.33) Chromosome with lowest proportion of DNA in annotated exons Chr. Y (0.36) Longest intergenic region (between annotated + hypothetical genes) Chr. 13 (3,038,416 bp) Rate of SNP variation 1/1250 bp * In these ranges, the percentages correspond to the annotated gene set (26, 383 genes) and the hypothetical + annotated gene set (39,114 genes), respectively.

4.2 Linkage map

We mapped the location of the markers that constitute the Genethon linkage map to the genome. The rate of recombination, expressed as cM per Mbp, was calculated for 3-Mbp windows as shown in Table 12. Higher rates of recombination in the telomeric region of the chromosomes have been previously documented (73). From this mapping result, there is a difference of 4.99 between lowest rates and highest rates and the largest difference of 4.4 between males and females (4.99 to 0.47 on chromosome 16). This indicates that the variability in recombination rates among regions of the genome exceeds the differences in recombination rates between males and females. The human genome has recombination hotspots, where recombination rates vary fivefold or more over a space of 1 kbp, so the picture one gets of the magnitude of variability in recombination rate will depend on the size of the window examined. Unfortunately, too few meiotic crossovers have occurred in Centre d'Étude du Polymorphism Humain (CEPH) and other reference families to provide a resolution any finer than about 3 Mbp. The next challenge will be to determine a sequence basis of recombination at the chromosomal level. An accurate predictor for the rate for variation in recombination rates between any pair of markers would be extremely useful in designing markers to narrow a region of linkage, such as in positional cloning projects.

Table 12. Rate of recombination per physical distance (cM/Mb) across the genome. Genethon markers were placed on CSA-mapped assemblies, and then relative physical distances and rates were calculated in 3-Mb windows for each chromosome. NA, not applicable. Chrom. Male Sex-average Female Max. Avg. Min. Max. Avg. Min. Max. Avg. Min.  1 2.60 1.12 0.23 2.81 1.42 0.52 3.39 1.76 0.68  2 2.23 0.78 0.33 2.65 1.12 0.54 3.17 1.40 0.61  3 2.55 0.86 0.23 2.40 1.07 0.42 2.71 1.30 0.33  4 1.66 0.67 0.15 2.06 1.04 0.60 2.50 1.40 0.77  5 2.00 0.67 0.18 1.87 1.08 0.42 2.26 1.43 0.62  6 1.97 0.71 0.28 2.57 1.12 0.37 3.47 1.67 0.64  7 2.34 1.16 0.48 1.67 1.17 0.47 2.27 1.21 0.34  8 1.83 0.73 0.14 2.40 1.05 0.46 3.44 1.36 0.43  9 2.01 0.99 0.53 1.95 1.32 0.77 2.63 1.66 0.82 10 3.73 1.03 0.22 3.05 1.29 0.66 2.84 1.51 0.76 11 1.43 0.72 0.31 2.13 0.99 0.47 3.10 1.32 0.49 12 4.12 0.76 0.26 3.35 1.16 0.49 2.93 1.55 0.59 13 1.60 0.75 0.01 1.87 0.95 0.17 2.49 1.19 0.32 14 3.15 0.98 0.18 2.65 1.30 0.62 3.14 1.63 0.75 15 2.28 0.94 0.34 2.31 1.22 0.42 2.53 1.56 0.54 16 1.83 1.00 0.47 2.70 1.55 0.63 4.99 2.32 1.12 17 3.87 0.87 0.00 3.54 1.35 0.54 4.19 1.83 0.94 18 3.12 1.37 0.86 3.75 1.66 0.43 4.35 2.24 0.72 19 3.02 0.97 0.10 2.57 1.41 0.49 2.89 1.75 0.87 20 3.64 0.89 0.00 2.79 1.50 0.83 3.31 2.15 1.34 21 3.23 1.26 0.69 2.37 1.62 1.08 2.58 1.90 1.18 22 1.25 1.10 0.84 1.88 1.41 1.08 3.73 2.08 0.93 X NA NA NA NA NA NA 3.12 1.64 0.72 Y NA NA NA NA NA NA NA NA NA Genome 4.12 0.88 0.00 3.75 1.22 0.17 4.99 1.55 0.32

4.3 Correlation between CpG islands and genes

Experimental methods have been used that resulted in an estimate of 30,000 to 45,000 CpG islands in the human genome (74, 80) and an estimate of 499 CpG islands on human chromosome 22 (81). Larsen et al. (76) and Gardiner-Garden and Frommer (75) used a computational method to identify CpG islands and defined them as regions of DNA of >200 bp that have a G+C content of >50% and a ratio of observed versus expected frequency of CG dinucleotide 0.6.

It is difficult to make a direct comparison of experimental definitions of CpG islands with computational definitions because computational methods do not consider the methylation state of cytosine and experimental methods do not directly select regions of high G+C content. However, we can determine the correlation of CpG island with gene starts, given a set of annotated genomic transcripts and the whole genome sequence. We have analyzed the publicly available annotation of chromosome 22, as well as using the entire human genome in our assembly and the computationally annotated genes. A variation of the CpG island computation was compared with Larsen et al. (76). The main differences are that we use a sliding window of 200 bp, consecutive windows are merged only if they overlap, and we recompute the CpG value upon merging, thus rejecting any potential island if it scores less than the threshold.

To compute various CpG statistics, we used two different thresholds of CG dinucleotide likelihood ratio. Besides using the original threshold of 0.6 (method 1), we used a higher threshold of CG dinucleotide likelihood ratio of 0.8 (method 2), which results in the number of CpG islands on chromosome 22 close to the number of annotated genes on this chromosome. The main results are summarized in Table 13. CpG islands computed with method 1 predicted only 2.6% of the CSA sequence as CpG, but 40% of the gene starts (start codons) are contained inside a CpG island. This is comparable to ratios reported by others (82). The last two rows of the table show the observed and expected average distance, respectively, of the closest CpG island from the first exon. The observed average closest CpG islands are smaller than the corresponding expected distances, confirming an association between CpG island and the first exon.

Table 13. Characteristics of CpG islands identified in chromosome 22 (34-Mbp sequence length) and the whole genome (2.9-Gbp sequence length) by means of two different methods. Method 1 uses a CG likelihood ratio of 0.6. Method 2 uses a CG likelihood ratio of 0.8. Chromosome 22 Whole genome (CS assembly) Method 1 Method 2 Method 1 Method 2 Number of CpG islands detected 5,211 522 195,706 26,876 Average length of island (bp) 390 535 395 497 Percent of sequence predicted as CpG 5.9 0.8 2.6 0.4 Percent of first exons that overlap a CpG island 44 25 42 22 Percent of first exons with first position of exon contained inside a CpG island 37 22 40 21 Average distance between first exon and closest CpG island (bp) 1,013 10,486 2,182 17,021 Expected distance between first exon and closest CpG island (bp) 3,262 32,567 7,164 55,811

We also looked at the distribution of CpG island nucleotides among various sequence classes such as intergenic regions, introns, exons, and first exons. We computed the likelihood score for each sequence class as the ratio of the observed fraction of CpG island nucleotides in that sequence class and the expected fraction of CpG island nucleotides in that sequence class. The result of applying method 1 on CSA were scores of 0.89 for intergenic region, 1.2 for intron, 5.86 for exon, and 13.2 for first exon. The same trend was also found for chromosome 22 and after the application of a higher threshold (method 2) on both data sets. In sum, genome-wide analysis has extended earlier analysis and suggests a strong correlation between CpG islands and first coding exons.

4.4 Genome-wide repetitive elements

Table 14. Distribution of repetitive DNA in the compartmentalized shotgun assembly sequence. Repetitive elements Megabases in assembled sequences Percent of assembly Previously predicted (%) (83) Alu 288 9.9 10.0 Mammalian interspersed repeat (MIR) 66 2.3 1.7 Medium reiteration (MER) 50 1.7 1.6 Long terminal repeat (LTR) 155 5.3 5.6 Long interspersed nucleotide element (LINE) 466 16.1 16.7 Total 1025 35.3 35.6

5 Genome Evolution

5.1 Retrotransposition in the human genome

A set of proteins corresponding to all 901 Otto-predicted, single-exon genes were subjected to BLAST analysis against the proteins encoded by the remaining multiexon predicted transcripts. Using homology criteria of 70% sequence identity over 90% of the length, we identified 298 instances of single- to multi-exon correspondence. Of these 298 sequences, 97 were represented in the GenBank data set of experimentally validated full-length genes at the stringency specified and were verified by manual inspection.

We believe that these 97 cases may represent intronless paralogs (see Web table 1 on Science Online at www.sciencemag.org/cgi/content/full/291/5507/1304/DC2) of known genes. Most of these are flanked by direct repeat sequences, although the precise nature of these repeats remains to be determined. All of the cases for which we have high confidence contain polyadenylated [poly(A)] tails characteristic of retrotransposition.

Recent publications describing the phenomenon of functional intronless paralogs speculate that retrotransposition may serve as a mechanism used to escape X-chromosomal inactivation (84, 86). We do not find a bias toward X chromosome origination of these retrotransposed genes; rather, the results show a random chromosome distribution of both the intron-containing and corresponding intronless paralogs. We also have found several cases of retrotransposition from a single source chromosome to multiple target chromosomes. Interesting examples include the retrotransposition of a five exon-containing ribosomal protein L21 gene on chromosome 13 onto chromosomes 1, 3, 4, 7, 10, and 14, respectively. The size of the source genes can also show variability. The largest example is the 31-exon diacylglycerol kinase zeta gene on chromosome 11 that has an intronless paralog on chromosome 13. Regardless of route, retrotransposition with subsequent gene changes in coding or noncoding regions that lead to different functions or expression patterns, represents a key route to providing an enhanced functional repertoire in mammals (87).

Our preliminary set of retrotransposed intronless paralogs contains a clear overrepresentation of genes involved in translational processes (40% ribosomal proteins and 10% translation elongation factors) and nuclear regulation (HMG nonhistone proteins, 4%), as well as metabolic and regulatory enzymes. EST matches specific to a subset of intronless paralogs suggest expression of these intronless paralogs. Differences in the upstream regulatory sequences between the source genes and their intronless paralogs could account for differences in tissue-specific gene expression. Defining which, if any, of these processed genes are functionally expressed and translated will require further elucidation and experimental validation.

5.2 Pseudogenes

We searched the complete set of Otto-predicted transcripts against the genomic sequence by means of BLAST. Genomic regions corresponding to all Otto-predicted transcripts were excluded from this analysis. We identified 2909 regions matching with greater than 70% identity over at least 70% of the length of the transcripts that likely represent processed pseudogenes. This number is probably an underestimate because specific methods to search for pseudogenes were not used.

We looked for correlations between structural elements and the propensity for retrotransposition in the human genome. GC content and transcript length were compared between the genes with processed pseudogenes (1177 source genes) versus the remainder of the predicted gene set. Transcripts that give rise to processed pseudogenes have shorter average transcript length (1027 bp versus 1594 bp for the Otto set) as compared with genes for which no pseudogene was detected. The overall GC content did not show any significant difference, contrary to a recent report (88). There is a clear trend in gene families that are present as processed pseudogenes. These include ribosomal proteins (67%), lamin receptors (10%), translation elongation factor alpha (5%), and HMG-non-histone proteins (2%). The increased occurrence of retrotransposition (both intronless paralogs and processed pseudogenes) among genes involved in translation and nuclear regulation may reflect an increased transcriptional activity of these genes.

5.3 Gene duplication in the human genome

5.4 Large-scale duplications

The first of the methods is based on the idea of searching for blocks of highly conserved homologous proteins that occur in more than one location on the genome. For this comparison, two genes were considered equivalent if their protein products were determined to be in the same family and the same complete Lek cluster (essentially paralogous genes) (89). Initially, each chromosome was represented as a string of genes ordered by the start codons for predicted genes along the chromosome. We considered the two strands as a single string, because local inversions are relatively common events relative to large-scale duplications. Each gene was indexed according to the protein family and Lek complete cluster (89). All pairs of indexed gene strings were then aligned in both the forward and reverse directions with the Smith-Waterman algorithm (90). A match between two proteins of the same Lek complete cluster was given a score of 10 and a mismatch 10, with gap open and extend penalties of 4 and 1. With these parameters, 19 conserved interchromosomal blocks of duplication were observed, all of which were also detected and expanded by the comprehensive method described below. The detection of only a relatively small number of block duplications was a consequence of using an intrinsically conservative method grounded in the conservative constraints of the complete Lek clusters.

In the second, more comprehensive approach, we aligned all chromosomes directly with one another using an algorithm based on the MUMmer system (91). This alignment method uses a suffix tree data structure and a linear-time algorithm to align long sequences very rapidly; for example, two chromosomes of 100 Mbp can be aligned in less than 20 min (on a Compaq Alpha computer) with 4 gigabytes of memory. This procedure was used recently to identify numerous large-scale segmental duplications among the five chromosomes of A. thaliana (92); in that organism, the method revealed that 60% of the genome (66 Mbp) is covered by 24 very large duplicated segments. For Arabidopsis, a DNA-based alignment was sufficient to reveal the segmental duplications between chromosomes; in the human genome, DNA alignments at the whole-chromosome level are insufficiently sensitive. Therefore, a modified procedure was developed and applied, as follows. First, all 26,588 proteins (9,675,713 million amino acids) were concatenated end-to-end in order as they occur along each of the 24 chromosomes, irrespective of strand location. The concatenated protein set was then aligned against each chromosome by the MUMmer algorithm. The resulting matches were clustered to extract all sets of three or more protein matches that occur in close proximity on two different chromosomes (93); these represent the candidate segmental duplications. A series of filters were developed and applied to remove likely false-positives from this set; for example, small blocks that were spread across many proteins were removed. To refine the filtering methods, a shuffled protein set was first created by taking the 26,588 proteins, randomizing their order, and then partitioning them into 24 shuffled chromosomes, each containing the same number of proteins as the true genome. This shuffled protein set has the identical composition to the real genome; in particular, every protein and every domain appears the same number of times. The complete algorithm was then applied to both the real and the shuffled data, with the results on the shuffled data being used to estimate the false-positive rate. The algorithm after filtering yielded 10,310 gene pairs in 1077 duplicated blocks containing 3522 distinct genes; tandemly duplicated expansions in many of the blocks explain the excess of gene pairs to distinct genes. In the shuffled data, by contrast, only 370 gene pairs were found, giving a false-positive estimate of 3.6%. The most likely explanation for the 1077 block duplications is ancient segmental duplications. In many cases, the order of the proteins has been shuffled, although proximity is preserved. Out of the 1077 blocks, 159 contain only three genes, 137 contain four genes, and 781 contain five or more genes.

To illustrate the extent of the detected duplications, Fig. 13 shows all 1077 block duplications indexed to each chromosome in 24 panels in which only duplications mapped to the indexed chromosome are displayed. The figure makes it clear that the duplications are ubiquitous in the genome. One feature that it displays is many relatively small chromosomal stretches, with one-to-many duplication relationships that are graphically striking. One such example captured by the analysis is the well-documented olfactory receptor (OR) family, which is scattered in blocks throughout the genome and which has been analyzed for genome-deployment reconstructions at several evolutionary stages (94). The figure also illustrates that some chromosomes, such as chromosome 2, contain many more detected large-scale duplications than others. Indeed, one of the largest duplicated segments is a large block of 33 proteins on chromosome 2, spread among eight smaller blocks in 2p, that aligns to a paralogous set on chromosome 14, with one rearrangement (see chromosomes 2 and 14 panels in Fig. 13). The proteins are not contiguous but span a region containing 97 proteins on chromosome 2 and 332 proteins on chromosome 14. The likelihood of observing this many duplicated proteins by chance, even over a span of this length, is 2.3 × 10⁶⁸ (93). This duplicated set spans 20 Mbp on chromosome 2 and 63 Mbp on chromosome 14, over 70% of the latter chromosome. Chromosome 2 also contains a block duplication that is nearly as large, which is shared by chromosome arm 2q and chromosome 12. This duplication incorporates two of the four known Hox gene clusters, but considerably expands the extent of the duplications proximally and distally on the pair of chromosome arms. This breadth of duplication is also seen on the two chromosomes carrying the other two Hox clusters.

An additional large duplication, between chromosomes 18 and 20, serves as a good example to illustrate some of the features common to many of the other observed large duplications (Fig. 13, inset). This duplication contains 64 detected ordered intrachromosomal pairs of homologous genes. After discounting a 40-Mb stretch of chromosome 18 free of matches to chromosome 20, which is likely to represent a large insert (between the gene assignments "Krup rel" and "collagen rel" on chromosome 18 in Fig. 13), the full duplication segment covers 36 Mb on chromosome 18 and 28 Mb on chromosome 20. By this measure, the duplication segment spans nearly half of each chromosome's net length. The most likely scenario is that the whole span of this region was duplicated as a single very large block, followed by shuffling owing to smaller scale rearrangements. As such, at least four subsequent rearrangements would need to be invoked to explain the relative insertions and inversions seen in the duplicated segment interval. The 64 protein pairs in this alignment occur among 217 protein assignments on chromosome 18, and among 322 protein assignments on chromosome 20, for a density of involved proteins of 20 to 30%. This is consistent with an ancient large-scale duplication followed by subsequent gene loss on one or both chromosomes. Loss of just one member of a gene pair subsequent to the duplication would result in a failure to score a gene pair in the block; less than 50% gene loss on the chromosomes would lead to the duplication density observed here. As an independent verification of the significance of the alignments detected, it can be seen that a substantial number of the pairs of aligning proteins in this duplication, including some of those annotated (Fig. 13), are those populating small Lek complete clusters (see above). This indicates that they are members of very small families of paralogs; their relative scarcity within the genome validates the uniqueness and robust nature of their alignments.

Two additional qualitative features were observed among many of the large-scale duplications. First, several proteins with disease associations, with OMIM (Online Mendelian Inheritance in Man) assignments, are members of duplicated segments (see web table 2 on Science Online at www.sciencemag.org/cgi/content/full/291/5507/1304/DC1). We have also observed a few instances where paralogs on both duplicated segments are associated with similar disease conditions. Notable among these genes are proteins involved in hemostasis (coagulation factors) that are associated with bleeding disorders, transcriptional regulators like the homeobox proteins associated with developmental disorders, and potassium channels associated with cardiovascular conduction abnormalities. For each of these disease genes, closer study of the paralogous genes in the duplicated segment may reveal new insights into disease causation, with further investigation needed to determine whether they might be involved in the same or similar genetic diseases. Second, although there is a conserved number of proteins and coding exons predicted for specific large duplicated spans within the chromosome 18 to 20 alignment, the genomic DNA of chromosome 18 in these specific spans is in some cases more than 10-fold longer than the corresponding chromosome 20 DNA. This selective accretion of noncoding DNA (or conversely, loss of noncoding DNA) on one of a pair of duplicated chromosome regions was observed in many compared regions. Hypotheses to explain which mechanisms foster these processes must be tested.

Evaluation of the alignment results gives some perspective on dating of the duplications. As noted above, large-scale ancient segmental duplication in fact best explains many of the blocks detected by this genome-wide analysis. The regions of human chromosomes involved in the large-scale duplications expanded upon above (chromosomes 2 to 14, 2 to 12, and 18 to 20) are each syntenic to a distinct mouse chromosomal region. The corresponding mouse chromosomal regions are much more similar in sequence conservation, and even in order, to their human synteny partners than the human duplication regions are to each other. Further, the corresponding mouse chromosomal regions each bear a significant proportion of genes orthologous to the human genes on which the human duplication assignments were made. On the basis of these factors, the corresponding mouse chromosomal spans, at coarse resolution, appear to be products of the same large-scale duplications observed in humans. Although further detailed analysis must be carried out once a more complete genome is assembled for mouse, the underlying large duplications appear to predate the two species' divergence. This dates the duplications, at the latest, before divergence of the primate and rodent lineages. This date can be further refined upon examination of the synteny between human chromosomes and those of chicken, pufferfish (Fugu rubripes), or zebrafish (95). The only substantial syntenic stretches mapped in these species corresponding to both pairs of human duplications are restricted to the Hox cluster regions. When the synteny of these regions (or others) to human chromosomes is extended with further mapping, the ages of the nearly chromosome-length duplications seen in humans are likely to be dated to the root of vertebrate divergence.

The MUMmer-based results demonstrate large block duplications that range in size from a few genes to segments covering most of a chromosome. The extent of segmental duplications raises the question of whether an ancient whole-genome duplication event is the underlying explanation for the numerous duplicated regions (96). The duplications have undergone many deletions and subsequent rearrangements; these events make it difficult to distinguish between a whole-genome duplication and multiple smaller events. Further analysis, focused especially on comparing the estimated ages of all the block duplications, derived partially from interspecies genome comparisons, will be necessary to determine which of these two hypotheses is more likely. Comparisons of genomes of different vertebrates, and even cross-phyla genome comparisons, will allow for the deconvolution of duplications to eventually reveal the stagewise history of our genome, and with it a history of the emergence of many of the key functions that distinguish us from other living things.

6 A Genome-Wide Examination of Sequence Variations

Having a complete genome sequence enables researchers to achieve a dramatic acceleration in the rate of gene discovery, but only through analysis of sequence variation in DNA can we discover the genetic basis for variation in health among human beings. Whole-genome shotgun sequencing is a particularly effective method for detecting sequence variation in tandem with whole-genome assembly. In addition, we compared the distribution and attributes of SNPs ascertained by three other methods: (i) alignment of the Celera consensus sequence to the PFP assembly, (ii) overlap of high-quality reads of genomic sequence (referred to as "Kwok"; 1,120,195 SNPs) (97), and (iii) reduced representation shotgun sequencing (referred to as "TSC"; 632,640 SNPs) (98). These data were consistent in showing an overall nucleotide diversity of ~8 × 10⁴, marked heterogeneity across the genome in SNP density, and an overwhelming preponderance of noncoding variation that produces no change in expressed proteins.

6.1 SNPs found by aligning the Celera consensus to the PFP assembly

6.2 Comparisons to public SNP databases

Table 15. Overlap of SNPs from genome-wide SNP databases. Table entries are SNP counts for each pair of data sets. Numbers in parentheses are the fraction of overlap, calculated as the count of overlapping SNPs divided by the number of SNPs in the smaller of the two databases compared. Total SNP counts for the databases are: Celera-PFP, 2,104,820; TSC, 585,811; and Kwok 438,032. Only unique SNPs in the TSC and Kwok data sets were included. TSC Kwok Celera-PFP 188,694 158,532 (0.322) (0.362) TSC 72,024 (0.164)

One means of assessing whether the three sets of SNPs provide the same picture of human variation is to tally the frequencies of the six possible base changes in each set of SNPs (Table 16). Previous measures of nucleotide diversity were mostly derived from small-scale analysis on candidate genes (101), and our analysis with all three data sets validates the previous observations at the whole-genome scale. There is remarkable homogeneity between the SNPs found in the Kwok set, the TSC set, and in our whole-genome shotgun (46) in this substitution pattern. Compared with the rest of the data sets, Celera-PFP deviates slightly from the 2:1 transition-to-transversion ratio observed in the other SNP sets. This result is not unexpected, because some fraction of the computationally identified SNPs in the Celera-PFP comparison may in fact be sequence errors. A 2:1 transition:transversion ratio for the bona fide SNPs would be obtained if one assumed that 15% of the sequence differences in the Celera-PFP set were a result of (presumably random) sequence errors.

Table 16. Summary of nucleotide changes in different SNP data sets. SNP data set A/G (%) C/T (%) A/C (%) A/T (%) C/G (%) T/G (%) Transition: transversion Celera-PFP 30.7 30.7 10.3 8.6 9.2 10.3 1.59:1 Kwok* 33.7 33.8 8.5 7.0 8.6 8.4 2.07:1 TSC 33.3 33.4 8.8 7.3 8.6 8.6 1.99:1 * November 2000 release of the NCBI database dbSNP (www.nci.nlm.nih.gov/SNP/) with the method defined as Overlap SnpDetectionWithPolyBayes. The submitter of the data is Pui-Yan Kwok from Washington University. November 2000 release of NCBI dbSNP (www.ncbi.nlm.nih.gov/SNP/) with the methods defined as TSC-Sanger, TSC-WICGR, and TSC-WUGSC. The submitter of the data is Lincoln Stein from Cold Spring Harbor Laboratory.

6.3 Estimation of nucleotide diversity from ascertained SNPs

Estimation of nucleotide diversity from a shotgun assembly entails calculating for each column of the multialignment, the probability that two or more distinct alleles are present, and the probability of detecting a SNP if in fact the alleles have different sequence (i.e., the probability of correct sequence calls). The greater the depth of coverage and the higher the sequence quality, the higher is the chance of successfully detecting a SNP (105). Even after correcting for variation in coverage, the nucleotide diversity appeared to vary across autosomes. The significance of this heterogeneity was tested by analysis of variance, with estimates of  for 100-kbp windows to estimate variability within chromosomes (for the Celera-PFP comparison, F = 29.73, P < 0.0001).

Average diversity for the autosomes estimated from the Celera-PFP comparison was 8.94 × 10⁴. Nucleotide diversity on the X chromosome was 6.54 × 10⁴. The X is expected to be less variable than autosomes, because for every four copies of autosomes in the population, there are only three X chromosomes, and this smaller effective population size means that random drift will more rapidly remove variation from the X (106).

Having ascertained nucleotide variation genome-wide, it appears that previous estimates of nucleotide diversity in humans based on samples of genes were reasonably accurate (101, 102, 106, 107). Genome-wide, our estimate of nucleotide diversity was 8.98 × 10^-4 for the Celera-PFP alignment, and a published estimate averaged over 10 densely resequenced human genes was 8.00 × 10⁴ (108).

6.4 Variation in nucleotide diversity across the human genome

somes, and whether this heterogeneity is greater than expected by chance. If SNPs occur by random and independent mutations, then it would seem that there ought to be a Poisson distribution of numbers of SNPs in fragments of arbitrary constant size. The observed dispersion in the distribution of SNPs in 100-kbp fragments was far greater than predicted from a Poisson distribution (Fig. 14). However, this simplistic model ignores the different recombination rates and population histories that exist in different regions of the genome. Population genetics theory holds that we can account for this variation with a mathematical formulation called the neutral coalescent (109). Applying well-tested algorithms for simulating the neutral coalescent with recombination (110), and using an effective population size of 10,000 and a per-base recombination rate equal to the mutation rate (111), we generated a distribution of numbers of SNPs by this model as well (112). The observed distribution of SNPs has a much larger variance than either the Poisson model or the coalescent model, and the difference is highly significant. This implies that there is significant variability across the genome in SNP density, an observation that begs an explanation.

Several attributes of the DNA sequence may affect the local density of SNPs, including the rate at which DNA polymerase makes errors and the efficacy of mismatch repair. One key factor that is likely to be associated with SNP density is the G+C content, in part because methylated cytosines in CpG dinucleotides tend to undergo deamination to form thymine, accounting for a nearly 10-fold increase in the mutation rate of CpGs over other dinucleotides. We tallied the GC content and nucleotide diversities in 100-kbp windows across the entire genome and found that the correlation between them was positive (r = 0.21) and highly significant (P < 0.0001), but G+C content accounted for only a small part of the variation.

6.5 SNPs by genomic class

It is striking how small is the fraction of SNPs that lead to potentially dysfunctional alterations in proteins. In the 10,239 RefSeq genes, missense SNPs were only about 0.12, 0.14, and 0.17% of the total SNP counts in Celera-PFP, TSC, and Kwok SNPs, respectively. Nonconservative protein changes constitute an even smaller fraction of missense SNPs (47, 41, and 40% in Celera-PFP, Kwok, and TSC). Intergenic regions have been virtually unstudied (113), and we note that 75% of the SNPs we identified were intergenic (Table 17). The SNP rate was highest in introns and lowest in exons. The SNP rate was lower in intergenic regions than in introns, providing one of the first discriminators between these two classes of DNA. These SNP rates were confirmed in the Celera SNPs, which also exhibited a lower rate in exons than in introns, and in extragenic regions than in introns (46). Many of these intergenic SNPs will provide valuable information in the form of markers for linkage and association studies, and some fraction is likely to have a regulatory function as well.

Table 17. Distribution of SNPs in classes of genomic regions. Genomic region class Size of region examined (Mb) Celera-PFP SNP density (SNP/Mb) Intergenic 2185 707 Gene (intron + exon) 646 917 Intron 615 921 First intron 164 808 Exon 31 529 First exon 10 592

7 An Overview of the Predicted Protein-Coding Genes in the Human Genome

A preliminary analysis of the predicted human protein-coding genes was conducted. Two methods were used to analyze and classify the molecular functions of 26,588 predicted proteins that represent 26,383 gene predictions with at least two lines of evidence as described above. The first method was based on an analysis at the level of protein families, with both the publicly available Pfam database (114, 115) and Celera's Panther Classification (CPC) (Fig. 15) (116). The second method was based on an analysis at the level of protein domains, with both the Pfam and SMART databases (115, 117).

The results presented here are preliminary and are subject to several limitations. Both the gene predictions and functional assignments have been made by using computational tools, although the statistical models in Panther, Pfam, and SMART have been built, annotated, and reviewed by expert biologists. In the set of computationally predicted genes, we expect both false-positive predictions (some of these may in fact be inactive pseudogenes) and false-negative predictions (some human genes will not be computationally predicted). We also expect errors in delimiting the boundaries of exons and genes. Similarly, in the automatic functional assignments, we also expect both false-positive and false-negative predictions. The functional assignment protocol focuses on protein families that tend to be found across several organisms, or on families of known human genes. Therefore, we do not assign a function to many genes that are not in large families, even if the function is known. Unless otherwise specified, all enumeration of the genes in any given family or functional category was taken from the set of 26,588 predicted proteins, which were assigned functions by using statistical score cutoffs defined for models in Panther, Pfam, and SMART.

For this initial examination of the predicted human protein set, three broad questions were asked: (i) What are the likely molecular functions of the predicted gene products, and how are these proteins categorized with current classification methods? (ii) What are the core functions that appear to be common across the animals? (iii) How does the human protein complement differ from that of other sequenced eukaryotes?

7.1 Molecular functions of predicted human proteins

In our analysis of the 12,731 additional low-confidence predicted genes (those with only one piece of supporting evidence), only 636 (5%) of these additional putative genes were assigned molecular functions by the automated methods. One-third of these 636 predicted genes represented endogenous retroviral proteins, further suggesting that the majority of these unknown-function genes are not real genes. Given that most of these additional 12,095 genes appear to be unique among the genomes sequenced to date, many may simply represent false-positive gene predictions.

The most common molecular functions are the transcription factors and those involved in nucleic acid metabolism (nucleic acid enzyme). Other functions that are highly represented in the human genome are the receptors, kinases, and hydrolases. Not surprisingly, most of the hydrolases are proteases. There are also many proteins that are members of proto-oncogene families, as well as families of "select regulatory molecules": (i) proteins involved in specific steps of signal transduction such as heterotrimeric GTP-binding proteins (G proteins) and cell cycle regulators, and (ii) proteins that modulate the activity of kinases, G proteins, and phosphatases.

7.2 Evolutionary conservation of core processes

We enumerated the "strict orthologs" conserved between human and fly, and between human and worm (Fig. 16) to address the question, What are the core functions that appear to be common across the animals? The concept of orthology is important because if two genes are orthologs, they can be traced by descent to the common ancestor of the two organisms (an "evolutionarily conserved protein set"), and therefore are likely to perform similar conserved functions in the different organisms. It is critical in this analysis to separate orthologs (a gene that appears in two organisms by descent from a common ancestor) from paralogs (a gene that appears in more than one copy in a given organism by a duplication event) because paralogs may subsequently diverge in function. Following the yeast-worm ortholog comparison in (120), we identified two different cases for each pairwise comparison (human-fly and human-worm). The first case was a pair of genes, one from each organism, for which there was no other close homolog in either organism. These are straightforwardly identified as orthologous, because there are no additional members of the families that complicate separating orthologs from paralogs. The second case is a family of genes with more than one member in either or both of the organisms being compared. Chervitz et al. (120) deal with this case by analyzing a phylogenetic tree that described the relationships between all of the sequences in both organisms, and then looked for pairs of genes that were nearest neighbors in the tree. If the nearest-neighbor pairs were from different organisms, those genes were presumed to be orthologs. We note that these nearest neighbors can often be confidently identified from pairwise sequence comparison without having to examine a phylogenetic tree (see legend to Fig. 16). If the nearest neighbors are not from different organisms, there has been a paralogous expansion in one or both organisms after the speciation event (and/or a gene loss by one organism). When this one-to-one correspondence is lost, defining an ortholog becomes ambiguous. For our initial computational overview of the predicted human protein set, we could not answer this question for every predicted protein. Therefore, we consider only "strict orthologs," i.e., the proteins with unambiguous one-to-one relationships (Fig. 16). By these criteria, there are 2758 strict human-fly orthologs, 2031 human-worm (1523 in common between these sets). We define the evolutionarily conserved set as those 1523 human proteins that have strict orthologs in both D. melanogaster and C. elegans.

10

The distribution of the functions of the conserved protein set is shown in Fig. 16. Comparison with Fig. 15 shows that, not surprisingly, the set of conserved proteins is not distributed among molecular functions in the same way as the whole human protein set. Compared with the whole human set (Fig. 15), there are several categories that are overrepresented in the conserved set by a factor of ~2 or more. The first category is nucleic acid enzymes, primarily the transcriptional machinery (notably DNA/RNA methyltransferases, DNA/RNA polymerases, helicases, DNA ligases, DNA- and RNA-processing factors, nucleases, and ribosomal proteins). The basic transcriptional and translational machinery is well known to have been conserved over evolution, from bacteria through to the most complex eukaryotes. Many ribonucleoproteins involved in RNA splicing also appear to be conserved among the animals. Other enzyme types are also overrepresented (transferases, oxidoreductases, ligases, lyases, and isomerases). Many of these enzymes are involved in intermediary metabolism. The only exception is the hydrolase category, which is not significantly overrepresented in the shared protein set. Proteases form the largest part of this category, and several large protease families have expanded in each of these three organisms after their divergence. The category of select regulatory molecules is also overrepresented in the conserved set. The major conserved families are small guanosine triphosphatases (GTPases) (especially the Ras-related superfamily, including ADP ribosylation factor) and cell cycle regulators (particularly the cullin family, cyclin C family, and several cell division protein kinases). The last two significantly overrepresented categories are protein transport and trafficking, and chaperones. The most conserved groups in these categories are proteins involved in coated vesicle-mediated transport, and chaperones involved in protein folding and heat-shock response [particularly the DNAJ family, and heat-shock protein 60 (HSP60), HSP70, and HSP90 families]. These observations provide only a conservative estimate of the protein families in the context of specific cellular processes that were likely derived from the last common ancestor of the human, fly, and worm. As stated before, this analysis does not provide a complete estimate of conservation across the three animal genomes, as paralogous duplication makes the determination of true orthologs difficult within the members of conserved protein families.

7.3 Differences between the human genome and other sequenced eukaryotic genomes

Molecular differences can be correlated with phenotypic differences to begin to reveal the developmental and cellular processes that are unique to the vertebrates. Tables 18 and 19 display a comparison among all sequenced eukaryotic genomes, over selected protein/domain families (defined by sequence similarity, e.g., the serine-threonine protein kinases) and superfamilies (defined by shared molecular function, which may include several sequence-related families, e.g., the cytokines). In these tables we have focused on (super) families that are either very large or that differ significantly in humans compared with the other sequenced eukaryote genomes. We have found that the most prominent human expansions are in proteins involved in (i) acquired immune functions; (ii) neural development, structure, and functions; (iii) intercellular and intracellular signaling pathways in development and homeostasis; (iv) hemostasis; and (v) apoptosis.

Table 18. Domain-based comparative analysis of proteins in H. sapiens (H), D. melanogaster (F), C. elegans (W), S. cerevisiae (Y), and A. thaliana (A). The predicted protein set of each of the above eukaryotic organisms was analyzed with Pfam version 5.5 using E value cutoffs of 0.001. The number of proteins containing the specified Pfam domains as well as the total number of domains (in parentheses) are shown in each column. Domains were categorized into cellular processes for presentation. Some domains (i.e., SH2) are listed in more than one cellular process. Results of the Pfam analysis may differ from results obtained based on human curation of protein families, owing to the limitations of large-scale automatic classifications. Representative examples of domains with reduced counts owing to the stringent E value cutoff used for this analysis are marked with a double asterisk (**). Examples include short divergent and predominantly alpha-helical domains, and certain classes of cysteine-rich zinc finger proteins. Accession number Domain name Domain description H F W Y A Developmental and homeostatic regulators PF02039 Adrenomedullin Adrenomedullin 1 0 0 0 0 PF00212 ANP Atrial natriuretic peptide 2 0 0 0 0 PF00028 Cadherin Cadherin domain 100 (550) 14 (157) 16 (66) 0 0 PF00214 Calc_CGRP_IAPP Calcitonin/CGRP/IAPP family 3 0 0 0 0 PF01110 CNTF Ciliary neurotrophic factor 1 0 0 0 0 PF01093 Clusterin Clusterin 3 0 0 0 0 PF00029 Connexin Connexin 14 (16) 0 0 0 0 PF00976 ACTH_domain Corticotropin ACTH domain 1 0 0 0 0 PF00473 CRF Corticotropin-releasing factor family 2 1 0 0 0 PF00007 Cys_knot Cystine-knot domain 10 (11) 2 0 0 0 PF00778 DIX Dix domain 5 2 4 0 0 PF00322 Endothelin Endothelin family 3 0 0 0 0 PF00812 Ephrin Ephrin 7 (8) 2 4 0 0 PF01404 EPh_Ibd Ephrin receptor ligand binding domain 12 2 1 0 0 PF00167 FGF Fibroblast growth factor 23 1 1 0 0 PF01534 Frizzled Frizzled/Smoothened family membrane region 9 7 3 0 0 PF00236 Hormone6 Glycoprotein hormones 1 0 0 0 0 PF01153 Glypican Glypican 14 2 1 0 0 PF01271 Granin Grainin (chromogranin or secretogranin) 3 0 0 0 0 PF02058 Guanylin Guanylin precursor 1 0 0 0 0 PF00049 Insulin Insulin/IGF/Relaxin family 7 4 0 0 0 PF00219 IGFBP Insulin-like growth factor binding proteins 10 0 0 0 0 PF02024 Leptin Leptin 1 0 0 0 0 PF00193 Xlink LINK (hyaluron binding) 13 (23) 0 1 0 0 PF00243 NGF Nerve growth factor family 3 0 0 0 0 PF02158 Neuregulin Neuregulin family 4 0 0 0 0 PF00184 Hormone5 Neurohypophysial hormones 1 0 0 0 0 PF02070 NMU Neuromedin U 1 0 0 0 0 PF00066 Notch Notch (DSL) domain 3 (5) 2 (4) 2 (6) 0 0 PF00865 Osteopontin Osteopontin 1 0 0 0 0 PF00159 Hormone3 Pancreatic hormone peptides 3 0 0 0 0 PF01279 Parathyroid Parathyroid hormone family 2 0 0 0 0 PF00123 Hormone2 Peptide hormone 5 (9) 0 0 0 0 PF00341 PDGF Platelet-derived growth factor (PDGF) 5 1 0 0 0 PF01403 Sema Sema domain 27 (29) 8 (10) 3 (4) 0 0 PF01033 Somatomedin_B Somatomedin B domain 5 (8) 3 0 0 0 PF00103 Hormone Somatotropin 1 0 0 0 0 PF02208 Sorb Sorbin homologous domain 2 0 0 0 0 PF02404 SCF Stem cell factor 2 0 0 0 0 PF01034 Syndecan Syndecan domain 3 1 1 0 0 PF00020 TNFR_c6 TNFR/NGFR cysteine-rich region 17 (31) 1 0 0 0 PF00019 TGF- Transforming growth factor -like domain 27 (28) 6 4 0 0 PF01099 Uteroglobin Uteroglobin family 3 0 0 0 0 PF01160 Opiods_neuropep Vertebrate endogenous opioids neuropeptide 3 0 0 0 0 PF00110 Wnt Wnt family of developmental signaling proteins 18 7 (10) 5 0 0 Hemostasis PF01821 ANATO Anaphylotoxin-like domain 6 (14) 0 0 0 0 PF00386 C1q C1q domain 24 0 0 0 0 PF00200 Disintegrin Disintegrin 18 2 3 0 0 PF00754 F5_F8_type_C F5/8 type C domain 15 (20) 5 (6) 2 0 0 PF01410 COLFI Fibrillar collagen C-terminal domain 10 0 0 0 0 PF00039 Fn1 Fibronectin type I domain 5 (18) 0 0 0 0 PF00040 Fn2 Fibronectin type II domain 11 (16) 0 0 0 0 PF00051 Kringle Kringle domain 15 (24) 2 2 0 0 PF01823 MACPF MAC/Perforin domain 6 0 0 0 0 PF00354 Pentaxin Pentaxin family 9 0 0 0 0 PF00277 SAA_proteins Serum amyloid A protein 4 0 0 0 0 PF00084 Sushi Sushi domain (SCR repeat) 53 (191) 11 (42) 8 (45) 0 0 PF02210 TSPN Thrombospondin N-terminal-like domains 14 1 0 0 0 PF01108 Tissue_fac Tissue factor 1 0 0 0 0 PF00868 Transglutamin_N Transglutaminase family 6 1 0 0 0 PF00927 Transglutamin_C Transglutaminase family 8 1 0 0 0 PF00594 Gla Vitamin K-dependent carboxylation/gamma- carboxyglutamic (GLA) domain 11 0 0 0 0 Immune response PF00711 Defensin_beta Beta defensin 1 0 0 0 0 PF00748 Calpain_inhib Calpain inhibitor repeat 3 (9) 0 0 0 0 PF00666 Cathelicidins Cathelicidins 2 0 0 0 0 PF00129 MHC_I Class I histocompatibility antigen, domains alpha 1 and 2 18 (20) 0 0 0 0 PF00993 MHC_II_alpha** Class II histocompatibility antigen, alpha domain 5 (6) 0 0 0 0 PF00969 MHC_II_beta** Class II histocompatibility antigen, beta domain 7 0 0 0 0 PF00879 Defensin_propep Defensin propeptide 3 0 0 0 0 PF01109 GM_CSF Granulocyte-macrophage colony-stimulating factor 1 0 0 0 0 PF00047 Ig Immunoglobulin domain 381 (930) 125 (291) 67 (323) 0 0 PF00143 Interferon Interferon alpha/beta domain 7 (9) 0 0 0 0 PF00714 IFN-gamma Interferon gamma 1 0 0 0 0 PF00726 IL10 Interleukin-10 1 0 0 0 0 PF02372 IL15 Interleukin-15 1 0 0 0 0 PF00715 IL2 Interleukin-2 1 0 0 0 0 PF00727 IL4 Interleukin-4 1 0 0 0 0 PF02025 IL5 Interleukin-5 1 0 0 0 0 PF01415 IL7 Interleukin-7/9 family 1 0 0 0 0 PF00340 IL1 Interleukin-1 7 0 0 0 0 PF02394 IL1_propep Interleukin-1 propeptide 1 0 0 0 0 PF02059 IL3 Interleukin-3 1 0 0 0 0 PF00489 IL6 Interleukin-6/G-CSF/MGF family 2 0 0 0 0 PF01291 LIF_OSM Leukemia inhibitory factor (LIF)/oncostatin (OSM) family 2 0 0 0 0 PF00323 Defensins Mammalian defensin 2 0 0 0 0 PF01091 PTN_MK PTN/MK heparin-binding protein 2 0 0 0 0 PF00277 SAA_proteins Serum amyloid A protein 4 0 0 0 0 PF00048 IL8 Small cytokines (intecrine/chemokine), interleukin-8 like 32 0 0 0 0 PF01582 TIR TIR domain 18 8 2 0 131 (143) PF00229 TNF TNF (tumor necrosis factor) family 12 0 0 0 0 PF00088 Trefoil Trefoil (P-type) domain 5 (6) 0 2 0 0 PI-PY-rho GTPase signaling PF00779 BTK BTK motif 5 1 0 0 0 PF00168 C2 C2 domain 73 (101) 32 (44) 24 (35) 6 (9) 66 (90) PF00609 DAGKa Diacylglycerol kinase accessory domain (presumed) 9 4 7 0 6 PF00781 DAGKc Diacylglycerol kinase catalytic domain (presumed) 10 8 8 2 11 (12) PF00610 DEP Domain found in Dishevelled, Egl-10, and Pleckstrin (DEP) 12 (13) 4 10 5 2 PF01363 FYVE FYVE zinc finger 28 (30) 14 15 5 15 PF00996 GDI GDP dissociation inhibitor 6 2 1 1 3 PF00503 G-alpha G-protein alpha subunit 27 (30) 10 20 (23) 2 5 PF00631 G-gamma G-protein gamma like domains 16 5 5 1 0 PF00616 RasGAP GTPase-activator protein for Ras-like GTPase 11 5 8 3 0 PF00618 RasGEFN Guanine nucleotide exchange factor for Ras-like GTPases; N-terminal motif 9 2 3 5 0 PF00625 Guanylate_kin Guanylate kinase 12 8 7 1 4 PF02189 ITAM Immunoreceptor tyrosine-based activation motif 3 0 0 0 0 PF00169 PH PH domain 193 (212) 72 (78) 65 (68) 24 23 PF00130 DAG_PE-bind Phorbol esters/diacylglycerol binding domain (C1 domain) 45 (56) 25 (31) 26 (40) 1 (2) 4 PF00388 PI-PLC-X Phosphatidylinositol-specific phospholipase C, X domain 12 3 7 1 8 PF00387 PI-PLC-Y Phosphatidylinositol-specific phospholipase C, Y domain 11 2 7 1 8 PF00640 PID Phosphotyrosine interaction domain (PTB/PID) 24 (27) 13 11 (12) 0 0 PF02192 PI3K_p85B PI3-kinase family, p85-binding domain 2 1 1 0 0 PF00794 PI3K_rbd PI3-kinase family, ras-binding domain 6 3 1 0 0 PF01412 ArfGAP Putative GTP-ase activating protein for Arf 16 9 8 6 15 PF02196 RBD Raf-like Ras-binding domain 6 (7) 4 1 0 0 PF02145 Rap_GAP Rap/ran-GAP 5 4 2 0 0 PF00788 RA Ras association (RalGDS/AF-6) domain 18 (19) 7 (9) 6 1 0 PF00071 Ras Ras family 126 56 (57) 51 23 78 PF00617 RasGEF RasGEF domain 21 8 7 5 0 PF00615 RGS Regulator of G protein signaling domain 27 6 (7) 12 (13) 1 0 PF02197 RIIa Regulatory subunit of type II PKA R-subunit 4 1 2 1 0 PF00620 RhoGAP RhoGAP domain 59 19 20 9 8 PF00621 RhoGEF RhoGEF domain 46 23 (24) 18 (19) 3 0 PF00536 SAM SAM domain (Sterile alpha motif) 29 (31) 15 8 3 6 PF01369 Sec7 Sec7 domain 13 5 5 5 9 PF00017 SH2 Src homology 2 (SH2) domain 87 (95) 33 (39) 44 (48) 1 3 PF00018 SH3 Src homology 3 (SH3) domain 143 (182) 55 (75) 46 (61) 23 (27) 4 PF01017 STAT STAT protein 7 1 1 (2) 0 0 PF00790 VHS VHS domain 4 2 4 4 8 PF00568 WH1 WH1 domain 7 2 2 (3) 1 0 Domains involved in apoptosis PF00452 Bcl-2 Bcl-2 9 2 1 0 0 PF02180 BH4 Bcl-2 homology region 4 3 0 1 0 0 PF00619 CARD Caspase recruitment domain 16 0 2 0 0 PF00531 Death Death domain 16 5 7 0 0 PF01335 DED Death effector domain 4 (5) 0 0 0 0 PF02179 BAG Domain present in Hsp70 regulators 5 (8) 3 2 1 5 PF00656 ICE_p20 ICE-like protease (caspase) p20 domain 11 7 3 0 0 PF00653 BIR Inhibitor of Apoptosis domain 8 (14) 5 (9) 2 (3) 1 (2) 0 Cytoskeletal PF00022 Actin Actin 61 (64) 15 (16) 12 9 (11) 24 PF00191 Annexin Annexin 16 (55) 4 (16) 4 (11) 0 6 (16) PF00402 Calponin Calponin family 13 (22) 3 7 (19) 0 0 PF00373 Band_41 FERM domain (Band 4.1 family) 29 (30) 17 (19) 11 (14) 0 0 PF00880 Nebulin_repeat Nebulin repeat 4 (148) 1 (2) 1 0 0 PF00681 Plectin_repeat Plectin repeat 2 (11) 0 0 0 0 PF00435 Spectrin Spectrin repeat 31 (195) 13 (171) 10 (93) 0 0 PF00418 Tubulin-binding Tau and MAP proteins, tubulin-binding 4 (12) 1 (4) 2 (8) 0 0 PF00992 Troponin Troponin 4 6 8 0 0 PF02209 VHP Villin headpiece domain 5 2 2 0 5 PF01044 Vinculin Vinculin family 4 2 1 0 0 ECM adhesion PF01391 Collagen Collagen triple helix repeat (20 copies) 65 (279) 10 (46) 174 (384) 0 0 PF01413 C4 C-terminal tandem repeated domain in type 4 procollagen 6 (11) 2 (4) 3 (6) 0 0 PF00431 CUB CUB domain 47 (69) 9 (47) 43 (67) 0 0 PF00008 EGF EGF-like domain 108 (420) 45 (186) 54 (157) 0 1 PF00147 Fibrinogen_C Fibrinogen beta and gamma chains, C-terminal globular domain 26 10 (11) 6 0 0 PF00041 Fn3 Fibronectin type III domain 106 (545) 42 (168) 34 (156) 0 1 PF00757 Furin-like Furin-like cysteine rich region 5 2 1 0 0 PF00357 Integrin_A Integrin alpha cytoplasmic region 3 1 2 0 0 PF00362 Integrin_B Integrins, beta chain 8 2 2 0 0 PF00052 Laminin_B Laminin B (Domain IV) 8 (12) 4 (7) 6 (10) 0 0 PF00053 Laminin_EGF Laminin EGF-like (Domains III and V) 24 (126) 9 (62) 11 (65) 0 0 PF00054 Laminin_G Laminin G domain 30 (57) 18 (42) 14 (26) 0 0 PF00055 Laminin_Nterm Laminin N-terminal (Domain VI) 10 6 4 0 0 PF00059 Lectin_c Lectin C-type domain 47 (76) 23 (24) 91 (132) 0 0 PF01463 LRRCT Leucine rich repeat C-terminal domain 69 (81) 23 (30) 7 (9) 0 0 PF01462 LRRNT Leucine rich repeat N-terminal domain 40 (44) 7 (13) 3 (6) 0 0 PF00057 Ldl_recept_a Low-density lipoprotein receptor domain class A 35 (127) 33 (152) 27 (113) 0 0 PF00058 Ldl_recept_b Low-density lipoprotein receptor repeat class B 15 (96) 9 (56) 7 (22) 0 0 PF00530 SRCR Scavenger receptor cysteine-rich domain 11 (46) 4 (8) 1 (2) 0 0 PF00084 Sushi Sushi domain (SCR repeat) 53 (191) 11 (42) 8 (45) 0 0 PF00090 Tsp_1 Thrombospondin type 1 domain 41 (66) 11 (23) 18 (47) 0 0 PF00092 Vwa von Willebrand factor type A domain 34 (58) 0 17 (19) 0 1 PF00093 Vwc von Willebrand factor type C domain 19 (28) 6 (11) 2 (5) 0 0 PF00094 Vwd von Willebrand factor type D domain 15 (35) 3 (7) 9 0 0 Protein interaction domains PF00244 14-3-3 14-3-3 proteins 20 3 3 2 15 PF00023 Ank Ank repeat 145 (404) 72 (269) 75 (223) 12 (20) 66 (111) PF00514 Armadillo_seg Armadillo/beta-catenin-like repeats 22 (56) 11 (38) 3 (11) 2 (10) 25 (67) PF00168 C2 C2 domain 73 (101) 32 (44) 24 (35) 6 (9) 66 (90) PF00027 cNMP_binding Cyclic nucleotide-binding domain 26 (31) 21 (33) 15 (20) 2 (3) 22 PF01556 DnaJ_C DnaJ C terminal region 12 9 5 3 19 PF00226 DnaJ DnaJ domain 44 34 33 20 93 PF00036 Efhand** EF hand 83 (151) 64 (117) 41 (86) 4 (11) 120 (328) PF00611 FCH Fes/CIP4 homology domain 9 3 2 4 0 PF01846 FF FF domain 4 (11) 4 (10) 3 (16) 2 (5) 4 (8) PF00498 FHA FHA domain 13 15 7 13 (14) 17 PF00254 FKBP FKBP-type peptidyl-prolyl cis-trans isomerases 15 (20) 7 (8) 7 (13) 4 24 (29) PF01590 GAF GAF domain 7 (8) 2 (4) 1 0 10 PF01344 Kelch Kelch motif 54 (157) 12 (48) 13 (41) 3 102 (178) PF00560 LRR** Leucine Rich Repeat 25 (30) 24 (30) 7 (11) 1 15 (16) PF00917 MATH MATH domain 11 5 88 (161) 1 61 (74) PF00989 PAS PAS domain 18 (19) 9 (10) 6 1 13 (18) PF00595 PDZ PDZ domain (Also known as DHR or GLGF) 96 (154) 60 (87) 46 (66) 2 5 PF00169 PH PH domain 193 (212) 72 (78) 65 (68) 24 23 PF01535 PPR** PPR repeat 5 3 (4) 0 1 474 (2485) PF00536 SAM SAM domain (Sterile alpha motif) 29 (31) 15 8 3 6 PF01369 Sec7 Sec7 domain 13 5 5 5 9 PF00017 SH2 Src homology 2 (SH2) domain 87 (95) 33 (39) 44 (48) 1 3 PF00018 SH3 Src homology 3 (SH3) domain 143 (182) 55 (75) 46 (61) 23 (27) 4 PF01740 STAS STAS domain 5 1 6 2 13 PF00515 TPR** TPR domain 72 (131) 39 (101) 28 (54) 16 (31) 65 (124) PF00400 WD40** WD40 domain 136 (305) 98 (226) 72 (153) 56 (121) 167 (344) PF00397 WW WW domain 32 (53) 24 (39) 16 (24) 5 (8) 11 (15) PF00569 ZZ ZZ-Zinc finger present in dystrophin, CBP/p300 10 (11) 13 10 2 10 Nuclear interaction domains PF01754 Zf-A20 A20-like zinc finger 2 (8) 2 2 0 8 PF01388 ARID ARID DNA binding domain 11 6 4 2 7 PF01426 BAH BAH domain 8 (10) 7 (8) 4 (5) 5 21 (25) PF00643 Zf-B_box** B-box zinc finger 32 (35) 1 2 0 0 PF00533 BRCT BRCA1 C Terminus (BRCT) domain 17 (28) 10 (18) 23 (35) 10 (16) 12 (16) PF00439 Bromodomain Bromodomain 37 (48) 16 (22) 18 (26) 10 (15) 28 PF00651 BTB BTB/POZ domain 97 (98) 62 (64) 86 (91) 1 (2) 30 (31) PF00145 DNA_methylase C-5 cytosine-specific DNA methylase 3 (4) 1 0 0 13 (15) PF00385 Chromo chromo' (CHRromatin Organization MOdifier) domain 24 (27) 14 (15) 17 (18) 1 (2) 12 PF00125 Histone Core histone H2A/H2B/H3/H4 75 (81) 5 71 (73) 8 48 PF00134 Cyclin Cyclin 19 10 10 11 35 PF00270 DEAD DEAD/DEAH box helicase 63 (66) 48 (50) 55 (57) 50 (52) 84 (87) PF01529 Zf-DHHC DHHC zinc finger domain 15 20 16 7 22 PF00646 F-box** F-box domain 16 15 309 (324) 9 165 (167) PF00250 Fork_head Fork head domain 35 (36) 20 (21) 15 4 0 PF00320 GATA GATA zinc finger 11 (17) 5(6) 8 (10) 9 26 PF01585 G-patch G-patch domain 18 16 13 4 14 (15) PF00010 HLH** Helix-loop-helix DNA-binding domain 60 (61) 44 24 4 39 PF00850 Hist_deacetyl Histone deacetylase family 12 5 (6) 8 (10) 5 10 PF00046 Homeobox Homeobox domain 160 (178) 100 (103) 82 (84) 6 66 PF01833 TIG IPT/TIG domain 29 (53) 11 (13) 5 (7) 2 1 PF02373 JmjC JmjC domain 10 4 6 4 7 PF02375 JmjN JmjN domain 7 4 2 3 7 PF00013 KH-domain KH domain 28 (67) 14 (32) 17 (46) 4 (14) 27 (61) PF01352 KRAB KRAB box 204 (243) 0 0 0 0 PF00104 Hormone_rec Ligand-binding domain of nuclear hormone receptor 47 17 142 (147) 0 0 PF00412 LIM LIM domain containing proteins 62 (129) 33 (83) 33 (79) 4 (7) 10 (16) PF00917 MATH MATH domain 11 5 88 (161) 1 61 (74) PF00249 Myb_DNA-binding Myb-like DNA-binding domain 32 (43) 18 (24) 17 (24) 15 (20) 243 (401) PF02344 Myc-LZ Myc leucine zipper domain 1 0 0 0 0 PF01753 Zf-MYND MYND finger 14 14 9 1 7 PF00628 PHD PHD-finger 68 (86) 40 (53) 32 (44) 14 (15) 96 (105) PF00157 Pou Pou domain--N-terminal to homeobox domain 15 5 4 0 0 PF02257 RFX_DNA_binding RFX DNA-binding domain 7 2 1 1 0 PF00076 Rrm RNA recognition motif (a.k.a. RRM, RBD, or RNP domain) 224 (324) 127 (199) 94 (145) 43 (73) 232 (369) PF02037 SAP SAP domain 15 8 5 5 6 (7) PF00622 SPRY SPRY domain 44 (51) 10 (12) 5 (7) 3 6 PF01852 START START domain 10 2 6 0 23 PF00907 T-box T-box 17 (19) 8 22 0 0 PF02135 Zf-TAZ TAZ finger 2 (3) 1 (2) 6 (7) 0 10 (15) PF01285 TEA TEA domain 4 1 1 1 0 PF02176 Zf-TRAF TRAF-type zinc finger 6 (9) 1 (3) 1 0 2 PF00352 TBP Transcription factor TFIID (or TATA-binding protein, TBP) 2 (4) 4 (8) 2 (4) 1 (2) 2 (4) PF00567 TUDOR TUDOR domain 9 (24) 9 (19) 4 (5) 0 2 PF00642 Zf-CCCH Zinc finger C-x8-C-x5-C-x3-H type (and similar) 17 (22) 6 (8) 22 (42) 3 (5) 31 (46) PF00096 Zf-C2H2** ZInc finger, C2H2 type 564 (4500) 234 (771) 68 (155) 34 (56) 21 (24) PF00097 Zf-C3HC4 Zinc finger, C3HC4 type (RING finger) 135 (137) 57 88 (89) 18 298 (304) PF00098 Zf-CCHC Zinc knuckle 9 (17) 6 (10) 17 (33) 7 (13) 68 (91)