De novo assembly of human genomes with massively parallel short read sequencing

Overall strategy for large genome assembly

We sequenced 200 Gb of Illumina GA reads for the Asian individual, including 72-Gb single-end and 128-Gb paired-end reads. The read lengths ranged from 35 bp to 75 bp, and the insert sizes of the paired-end libraries were 140 bp, 440 bp, 2.6 kb, 6 kb, and 9.6 kb (Supplemental Table 1). To manage the huge number of short reads effectively and handle them in a standard supercomputer with 512 Gb memory installed, we modularized the assembly method and organized it as a pipeline by loading only the necessary data at each step. Since some (∼5%) chimeric reads in long paired-end (≥2 kb) sequencing are generated in the circularizing and fragmentation process (Bentley et al. 2008), we only used single-end and paired-end reads with insert sizes of 140 bp and 440 bp for contig assembly; all paired-end data were used for scaffold construction.

We used genomic DNA to construct sequencing libraries and generated short reads from both ends of the clones (Fig. 2A). The read sequences were loaded into the computer and de Bruijn graph data structure was used to represent the overlap among the reads (Fig. 2B). Next, erroneous connections in the graph were removed, tiny repeats were resolved by read path, and the graph was simplified by merging unambiguously connected nodes into one (Fig. 2C). There are three major types of erroneous connections that need to be addressed: (1) tips—short and low-coverage dead ends, which are likely to be caused by sequencing errors at read ends; (2) low-coverage links—nodes connected by only one or a few reads, which are likely to be chimeric connections; and (3) bubbles—redundant paths with minor differences, which may represent polymorphisms between either homogenous chromosomes or repeat copies.

Open in a separate window

Figure 2.

Schematic overview of the assembly algorithm. (A) Genomic DNA was fragmented randomly and sequenced using paired-end technology. Short clones with sizes between 150 and 500 bp were amplified and sequenced directly; while long range (2–10 kb) paired-end libraries were constructed by circularizing DNA, fragmentation, and then purifying fragments with sizes in the range of 400–600 bp for cluster formation. (B) The raw or precorrected reads were then loaded into computer memory and de Bruijn graph data structure was used to represent the overlap among the reads. (C) The graph was simplified by removing erroneous connections (in red color on the graph) and solving tiny repeats by read path: (i) Clipping the short tips, (ii) removing low-coverage links, (iii) solving tiny repeats by read path, and (iv) merging the bubbles that were caused by repeats or heterozygotes of diploid chromosomes. (D) On the simplified graph, we broke the connections at repeat boundaries and output the unambiguous sequence fragments as contigs. (E) We realigned the reads onto the contigs and used the paired-end information to join the unique contigs into scaffolds. (F) Finally, we filled in the intrascaffold gaps, which were most likely comprised by repeats, using the paired-end extracted reads.

Once these erroneous edges were corrected, repeat connections on the graph were broken and linear sequences were output as contigs (Fig. 2D). By realigning the reads onto the contigs and transferring read pairing onto contig pairing relationships, we ordered the unique contigs and constructed them into scaffolds (Fig. 2E). Finally, the intrascaffold gaps were filled in through local assembly of the extracted reads inside the gap regions using paired-end information (Fig. 2F).

Detailed steps for genome assembly

Preassembly sequencing error correction

The rate of sequencing errors in Illumina GA reads is about 1%–2%. Even though errors primarily accumulate at the 3′-end of reads, many of the 25-mers will contain errors, which will make the total number of 25-mers much greater than expected. Error correction before assembly for small data sets is less important (and therefore optional) since the erroneous connections can easily be removed in the graph during assembly. This step, however, is essential for large data sets, as doing so tremendously reduces memory usage, making it feasible to load the complete number of read sequences and construct the de Bruijn graph.

For the Asian genome data, the total number of distinct 25-mers was reduced from 14.6 billion to 5.0 billion (2.9 timers smaller) through this correction (Table 1). With the majority of the errors corrected, the ratio of error-free reads increased from 60.1% to 74.0%. A very small fraction (0.29%) of the reads might have been incorrectly revised in regions where sequence coverage was not deep enough, but these are unlikely to cause misassembly since paired-end information will be used in a later step to confirm the sequence overlap.

Table 1.

Summary of preassembly error correction in the Asian genome sequencing

Open in a separate window

Scaffolding

After obtaining the contig sequences, we realigned the short reads onto the contigs. Since the repeat copies had been merged into consensus sequences in the graph and in the output contigs, each short read always mapped unambiguously to one contig. We used a minimum of three read pairs as the criteria to define the order and distance between two contigs, so the small fraction of chimeric reads will not create misassembly. Then the relationship among all the contigs was displayed as a graph. We masked the repeat contigs that had multiple and conflicting connections to the unique contigs. The remaining contigs with compatible connections to each other were linearized and constructed into scaffolds.

Starting from small and moving to larger insert sizes, we used the read mate pairs to join contigs into scaffolds step by step. By adopting 140- and 440-bp insert size paired-ends, the N50 of constructed scaffolds was 17.3 kb (Table 2). Adding 2.6 kb insert size paired-ends, the N50 size was improved to 103.5 kb. Adding 6- and 9.6-kb insert size libraries, the N50 of final scaffolds reached 446.3 kb. As shown in theoretical estimation, further improvement of N50 scaffold size depends on even larger insert size libraries.

Table 2.

Summary of the African and Asian genome assembly

Open in a separate window

All read sequences were used in contig assembly, while paired-end libraries with different insert sizes were used step-by-step additively on scaffold construction. N50 of contig or scaffold was calculated by ordering all sequences, then adding the lengths from longest to shortest until the summed length exceeded 50% of the total length of all sequences. N90 is similarly defined. NCBI build 36.1 was used as the reference genome and RefSeq was used as the gene set to evaluate genome and gene region coverage. Since both genomes were sequenced of male individuals, chromosomes X and Y only have half-sequencing depths of the autosomes, and hence were excluded in calculation genome and gene coverage. For calculating scaffold N50 and total length, the intrascaffold gaps were included.

Gap closure

The majority of the gaps inside the scaffolds were composed of repeats that were masked during scaffold construction. To disassemble the repeat copies and fill in the gaps, we used the paired-end information to retrieve the read pairs that had one read well-aligned on the contigs and another read located in the gap region, then, a local assembly for the collected reads was done. We closed 83.5% of the 6.3 million intrascaffold gaps, or 45.0% of the 717-Mb sum gap length (Table 3). The contig N50 size grew from 1050 bp to 7.4 kb (if we ignore gaps <50 bp in length, the N50 size was 11.5 kb), and the genome coverage improved from 80.3% to 87.4% (Table 2).

Table 3.

Percentage of the intrascaffold gaps that were closed

Open in a separate window

Summary and comparison of the two assembled genomes

We applied the same assembly method to the African genome, and the final assembly of both genomes is summarized in Table 2. The assembly of the Asian genome had a longer N50 size than the African genome for contigs (7.4 kb vs. 5.9 kb) and for scaffolds (446.3 kb vs. 61.9 kb), which is likely due to the longer average read length (55 bp vs. 35 bp) and longer paired-end insert sizes (9.6 kb vs. 2 kb) of the Asian genome sequencing data.

Although there is structural variation between human genomes (Bentley et al. 2008; Kidd et al. 2008; Wang et al. 2008), comparison of the Asian and African genome assemblies against the NCBI human reference genome gave us a general assessment of genome coverage and assembly accuracy for both. The Asian sequence assembly had higher coverage of the NCBI reference genome than did the African genome (87.4% vs. 85.4%). This may be due to the Asian genome assembly having a longer total length, and may also be because the Asian genome is more similar to the NCBI reference genome than is the African genome (International HapMap Consortium 2007; Bentley et al. 2008; Wang et al. 2008). The lower coverage of the NCBI reference genome by the assembly method than by the mapping-based method (92% coverage) (Wang et al. 2008) can be explained by the fact that regions with insufficient sequence depth cannot be assembled and that small contigs (<100 bp) were filtered in the final assembly. According to the location of the RefSeq genes on the NCBI reference genome, the Asian and African genome assemblies covered 95.5% and 89.2% of the gene region, respectively.

Sequence accuracy of the assemblies

By mapping all the reads onto the assembled genomes, we calculated the allele frequency in the reads at each genomic location to measure assembly quality at a single-base level. The peak read depth of the Asian and African genome was 55 and 40, and over 99% of both assembled genomes had more than 20-read coverage (Supplemental Fig. 2). This indicated that deep short read sequencing has a higher base-level sequence accuracy than traditional Sanger sequencing, which normally has 4–10× coverage.

Since we previously annotated the SNPs between the Asian individual and the NCBI reference sequences (Wang et al. 2008), we aligned the Asian genome assembly with the NCBI reference to detect differing alleles and checked the overlap of these alleles to the previously identified SNPs in the Asian genome. There were 1.87 million mismatched alleles that comprised 0.09% of the aligned region. Only 78 alleles (0.004%) were inconsistent with the annotated SNPs.

Structural difference to the NCBI reference genome

From the comparison between the assembled genomes and the NCBI reference sequence, we observed structural differences that could be structural variations or misassemblies. To distinguish these two categories of differences and evaluate the rate of misassembly, we checked the number of supportive paired-ends and conflicting paired-ends to the assembly at the discrepant regions (Supplemental Fig. 3).

There were 2195 and 2406 contigs in the Asian and African genome that showed greater than 100-bp insertion or deletion against the NCBI reference sequence (Supplemental Table 3). The insert sizes of paired-ends were consistent with the span of their alignment on the assembly at these regions, so these insertions or deletions are more likely to be true. There were 117 and 3339 small contigs in the Asian and African assembly that failed to be placed into the gaps of scaffolds due to insufficient paired-end information at the regions. We found 3516 and 3339 contig insertions in Asian and African genomes compared to the NCBI reference. Only eight (0.2%) and three (0.1%) cases were potential misassemblies that have clearly more (over two times) conflicting than supportive paired-ends, while the others were true insertions. During the scaffolding process, some flanking pairs of small contigs may have been placed in an incorrect order because the contig sizes were in the same range of paired-end insert size deviation and were thus difficult to order. We observed 4715 and 3094 such cases in the Asian and African genomes, of which 1681 (35.7%) and 593 (19.2%) were likely misassembled. Including inversions and long-range translocations, the total length of discrepant regions were about 6 Mb and 5 Mb in the Asian and African genome, comprising 0.3% and 0.25% of the assembly, respectively.

Carrying out de novo assembly allowed us to identify small deletions and insertions; whereas, this is not possible by mapping-based methods when the length of deletion or insertion is comparable or smaller than the standard deviation of paired-end insert sizes. De novo assembly also has the advantage of resolving structural variations to a single-base level and obtaining the inserted sequences. Figure 3A shows an example of a detected 17-bp deletion. Figure 3B shows a detected 7926-bp insertion. Case-by-case analysis and further experiment validation will be required to fully characterize all of the structural variations.

Open in a separate window

Figure 3.

Examples of deletion and insertion identified in the comparison of the assembled individual human genomes and the NCBI reference genome. (A) A 17-bp deletion in scaffold27122121 of the African genome located on chromosome 7. (B) A 7926-bp insertion in scaffold4928 of the Asian genome located on chromosome 7. The inserted sequence fragment was validated by a human BAC clone AC153461.2 in GenBank, and also exists in the chimpanzee genome.

Sequence depth effect on genome assembly

To determine the minimal sequence depth required for achieving a proper assembly of the human genome, we randomly sampled subsets of reads with different average depths from the error corrected reads of the Asian genome. This showed that contig size increased nearly linearly by improving sequence depth from 10× to 30×. Contig N50 and N90 size at 30× sequence depth were 1035 bp and 201 bp, respectively. We found that further improvement of sequence depth (to 40× and 50×) made only slight changes in contig size (N50 size of 1045 bp and 1050 bp) (Fig. 4). Accordingly, the total length of assembled contigs at 10×, 20×, and 30× sequence depth was 1.70, 2.09, and 2.15 Gb, respectively.

Open in a separate window

Figure 4.

N50 and N90 size of assembled contigs by different sequence depths. We sampled subsets of randomly selected reads from the Asian genome data for de novo assembly of contigs. The same K-mer (K = 25) size was used for all the assemblies.

We also simulated reads with different lengths to investigate the optimal sequence depth for each read length. With a 35-bp read length, 30–50× provided the best results; with a 50-bp read length, 30× was best; and with a 75-bp read length, 20× sequence depth was sufficient (Supplemental Table 4). Considering that in real sequencing experiments DNA fragmentation is not completely random and there are unavoidable sequencing errors, we suggest that 5–10× more reads than the theoretical estimate would be best for achieving optimal assembly.

As was shown in Figure 1B and Table 1, the length of paired-end insert sizes determined the scaffold size of de novo genome assembly. SOAPdenovo required a minimum of 50× physical clone coverage for each of the libraries of ∼200 bp, ∼500 bp, ∼2 kb, ∼5 kb, ∼10 kb, etc.

Error correction

For deep sequencing, the correct K-mers appear multiple times in the reads set, while random sequencing error-containing K-mers have low frequency. Our error correction method used K-mer frequency information. The K-mer size was determined by the genome size, read length, and supercomputer memory. Since we only need to correct the low-frequency K-mer, to save memory, we used one byte for each K-mer to store the frequency and assign all counts over 255 as 255. Here, we chose K = 17 bp, because 4¹⁷ = 16 G is larger than the genome size; thus, error-containing 17-mers were unlikely to exist in real genomes. The peak frequency of correct 17-mers would be about 20 in the read sets. We built a hash table to store the frequency of all 17-mers, which occupied 16 Gb of memory. Then, for each read, we started from high-frequency regions and extended both sides to infer potential erroneous sites of low-frequency (<3) 17-mers. For each inferred erroneous site, we tested the impact of changing it to the other three allele types, and these changes were picked up as candidates if all 17-mers containing the allele had a frequency equal to or over 3. If we obtained no candidates that satisfied these criteria, we did not change it; otherwise, the allele was revised to that with the highest 17-mer frequency. A dynamic programming algorithm was used to find the optimal solution with minimal changes. To increase speed, we used threaded parallelization to split read sets and handled them in parallel by sharing the same 17-mer hash table.

De Bruijn graph construction

For the de Bruijn graph, each node is a K-mer. Two nodes that overlap K − 1 bp and appeared in a neighboring read sequence were connected as an edge. Small K-mers make the graph very complex with a lot of edges created by repeat sequences; while large K-mers can have poor overlap in regions with low sequencing depth. After assessing different K-mer sizes, we found that the 25-mer provided the best tradeoff.

Tip removal

A single base-pair sequencing error on a read will create K consecutive incorrect K-mers. If the error occurs in the middle of a read and the K-mers at both ends are correct, the path created by the error would appear as a bubble (the bubble is discussed below) in the graph; otherwise, it would cause a “dead end,” or a tip, in the graph. We removed the tips that were shorter than 2K (50 bp if K = 25) and had a lower frequency than other alternative paths that connected at a common destination node in the graph.

Solving tiny repeats

Tiny repetitive sequences in the graph that are shorter than the read length may be able to be resolved by read paths. To avoid misassembly, we only tried to solve repeat nodes with equal N incoming and outgoing edges. If each of the N incoming edges had read path support from one of the N outgoing edges and had no conflicts, we then removed the repeat node and split the connections into N parallel paths.

Merging bubbles

We used Dijkstra's algorithm to detect bubbles, which is similar to the “Tour-bus” method in Velvet. We merged the detected bubbles into a single path if the sequences of the parallel paths were very similar; that is, only had a single base pair difference or had fewer than four base pairs difference with >90% identity.

Contig linkage graph

The first step of scaffolding is to realign the reads onto the contig sequences. Then the paired-end relationship between the reads was transferred to linkage between contigs. The linkages among all contigs formed a graph. We used the number of read pairs between two contigs to weight the linkage, and used the read paired-end insert sizes to estimate the gap size between the two contigs. Theoretically, if the insert size of a paired-end clone library obeys Normal distribution with a mean value μ and variance σ², the gap size estimated from N paired-ends will also have mean μ, but variance σ²/N. For our method, we required as least three read pairs to form a linkage.

This project is supported by the Chinese Academy of Science (GJHZ0701-6; KSCX2-YWN-023), the National Natural Science Foundation of China (30725008; 30890032; 90608010, 30811130531), the Chinese 973 program (2007CB815701; 2007CB815703; 2007CB815705), the Chinese 863 program (2006AA02Z177; 2006AA10A121), the International Science and Technology Cooperation Project (0806), the Shenzhen hundred-million project, the Danish Platform for Integrative Biology, the Ole Rømer grant from the Danish Natural Science Research Council, the pig bioinformatics grant from Danish Research Council, the Solexa project (272-07-0196), the Danish Strategic Research Council grant no. 2106-07-0021 (Seqnet), and the Lundbeck Foundation Centre of Applied Medical Genomics for Personalized Disease Prediction, Prevention and Care (LUCAMP). We thank Laurie Goodman for editing the manuscript.