Generating a high-quality genome sequence is a critical foundation for the analysis of any organism. Yet it remains a challenge, especially for genomes containing substantial repetitive sequences, such as Aedes aegypti, the principal vector of the Zika virus. Recently, an international consortium was organized to better understand Zika’s principal vector by improving the quality of the Ae. aegypti genome (1).

Currently, most genomes are assembled from a deep collection of short DNA sequence reads. These data are combined with linking information, which makes it possible to estimate the distance between individual sequences; such linking information is typically obtained by sequencing paired ends from a DNA clone library with a characteristic insert size. On the basis of sequence overlap, the reads are assembled into contiguous sequences (contigs); by means of the linking information, the contigs are ordered and oriented with respect to one another into larger scaffolds (2). Within scaffolds, adjacent contigs are often separated by a gap, which corresponds to a region that is hard to assemble from the available sequence reads (for example, because of repetitive sequences or low coverage) but that can nevertheless be spanned by using the linking information to determine the contigs at either end of the gap. Long links, from large-insert clones such as Fosmids, have been especially valuable (3). Such clone libraries provide physical coverage (defined as the average number of clones spanning a point in the genome), often in the range of 1000-fold. With this strategy, it has been possible to produce mammalian genome assemblies with scaffolds ranging from 1 to 15 Mb (2, 3). However, it has generally not been feasible to achieve scaffolds that span entire chromosomes, because some repetitive regions are too large and difficult to be spanned by the available clone libraries.

Hi-C is a sequencing-based approach for determining how a genome is folded by measuring the frequency of contact between pairs of loci (4, 5). Contact frequency depends strongly on the one-dimensional (1D) distance, in base pairs, between a pair of loci. For instance, loci separated by 10 kb in the human genome form contacts eight times more often than those at a distance of 100 kb. In absolute terms, a typical distribution of Hi-C contacts from a given locus is 15% to loci within 10 kb; 15% to loci 10 kb to 100 kb away; 18% to loci 100 kb to 1 Mb away; 13% to loci 1 to 10 Mb away; 16% to loci 10 to 100 Mb away; 2% to loci on the same chromosome but more than 100 Mb away; and 21% to loci on a different chromosome.

Hi-C data can provide links across a variety of length scales, spanning even whole chromosomes. However, unlike paired-end reads from clone libraries, any given Hi-C contact spans an unknown length and may connect loci on different chromosomes. This challenge may be mitigated, in part, by the physical coverage achieved by Hi-C data sets. For the maps reported in (4, 5), summing the span of each individual contact reveals that 1× of sequence coverage of the target genome translates, on average, into 23,000× of physical coverage. This suggests that a statistical approach analyzing the pattern of Hi-C contacts as a whole could generate extremely long scaffolds.

Computational experiments with simulated input data have suggested that Hi-C should be able to produce chromosome-length scaffolds (6–8). Indeed, Hi-C has been used to improve draft genome assemblies (7, 9) and to create chromosome-length scaffolds for large genomes (10). In this process, Hi-C data are used to assign draft scaffolds to chromosomes and then to order and orient the draft scaffolds within each chromosome. Unfortunately, the resulting predictions contain large errors, including chromosome-scale inversions and misjoins that fuse chromosomes (10). Such mis-assemblies may be caused by errors in the original draft assembly (10). One approach to avoiding such errors might be additional types of information, such as longer reads or optical mapping data [see, e.g., (11, 12)].

We therefore sought to develop a robust procedure for using Hi-C linking information to generate accurate genome assemblies with chromosome-length scaffolds. A key aspect of the approach is to first use Hi-C data to identify and correct errors in the scaffolds of the initial assembly. Briefly, we correct misjoins by identifying positions where a scaffold’s long-range contact pattern changes abruptly, which is unlikely for a correctly assembled scaffold. Next, we use a novel algorithm to anchor, order, and orient the resulting sequences, employing the contact frequency between a pair of sequences as an indicator of their proximity in the 1D genome. Finally, we merge contigs and scaffolds that correspond to overlapping regions of the genome by identifying pairs of scaffolds that exhibit both strong sequence homology as well as strong similarity in long-range contact pattern (Fig. 1 and fig. S1).

Open in a separate window

Fig. 1

Starting with a draft assembly, we used Hi-C data to correct mis-joins, scaffold, and merge overlaps, thereby generating an assembly of the Ae. aegypti mosquito genome with chromosome-length scaffolds

Here we show contact matrices generated by aligning a Hi-C data set to both the AaegL2 assembly (18) that we used as input (left) and the final AaegL4 assembly generated by our algorithm (right). Pixel intensity in the contact matrix indicates how often a pair of loci colocate in the nucleus. The loci corresponding to each row and column are illustrated using chromograms. The chromograms on the left depict the three linkage groups [Lnk1, Lnk2, Lnk3, or unassigned (U)] reported in AaegL2; the chromograms on the right depict the three chromosome-length scaffolds in AaegL4 (chr1, chr2, and chr3). To create the chromogram, we assigned each AaegL4 arm a linear color gradient, thereby specifying a color for each AaegL4 locus. The same colors are then used for the corresponding loci in AaegL2 (left) and in the illustration of our procedure (center, though with increased contrast). Chromogram discontinuities indicate differences with AaegL4. In the center, we illustrate our assembly algorithm using an input scaffold from Lnk1 of AaegL2 (“supercontig 1.12,” see bracket). First, the scaffold is examined for misjoins and split such that the resulting segments each exhibit a continuous Hi-C signal (center, top row). Next, the segments are used as input for iterative scaffolding. Ultimately, only one of the segments is assigned to chromosome 1 of AaegL4. The rest of supercontig 1.12 is assigned to 2q, in the vicinity of several scaffolds that were not anchored in AaegL2 (center, middle row). Finally, segments exhibiting a similar 3D signal are examined for evidence of overlapping sequence (green rectangle) and merged (center, bottom row). The final contact map is consistent with the Rabl configuration, i.e., the spatial clustering of centromeres and telomeres.

We validated our approach by creating a de novo assembly of a human genome (the GM12878 cell line), comprising 23 chromosome-length scaffolds, using only short Illumina reads (67× coverage). We created a draft assembly from 250–base pair (bp) paired-end reads [60× coverage, generated by Illumina sequencing with a polymerase chain reaction (PCR)–free protocol, downloaded from Sequence Read Archive (SRX297987); assembled with DISCOVAR de novo (13)]. This assembly, termed Hs1, comprises 2.82 Gb of sequence (contig N50 length: 103 kb) partitioned among 73,770 scaffolds (scaffold N50: 126 kb, Table 1).

We then used in situ Hi-C data (6.7× sequence coverage) to improve Hs1. We set aside the tiny scaffolds (43,231 scaffolds shorter than 15 kb, whose N50 length is 6.1 kb). Together, these contain 5.4% of sequenced bases in Hs1. Because of their small size, they have relatively few Hi-C contacts and are more difficult to analyze. We then used Hi-C data to split, anchor, order, and orient the remaining 30,539 scaffolds.

The resulting assembly (Hs2-HiC) consisted of 23 huge scaffolds (lengths from 28.8 to 225.2 Mb) containing 99.5% of the total sequence, together with an additional 811 small scaffolds (N50 length of 30 kb; maximum length of 231 kb) making up the remaining 0.5% of the genome (Table 1 and tables S1 to S6). Crucially, the assembly was generated entirely de novo.

We assessed the quality of Hs2-HiC by comparing it to the human genome reference, hg38 (fig. S9). The 23 scaffolds correspond to the 23 human chromosomes, spanning 99% of the length and containing 91% of the sequence in the chromosome-length scaffolds (table S1). These scaffolds are comparable in length to those reported by the International Human Genome Sequencing Consortium (14) and longer than those reported by (15).

Of the 29,344 scaffolds that were incorporated into chromosome-length scaffolds in Hs2-HiC and that could be uniquely placed in hg38, 99.70% (comprising 99.88% of the sequenced bases) were assigned to the correct chromosome. For randomly selected pairs of scaffolds assigned to the same chromosome-length scaffold in Hs2-HiC, the order in Hs2-HiC agreed with the order in hg38 in 99% of cases. The agreement was 96% for pairs of scaffolds that were adjacent in Hs2-HiC, reflecting the fact that the Hi-C data provide less information to resolve the fine-structure order of short scaffolds. However, the agreement was 99% for scaffolds of at least 120 kb in length. Similarly, the orientation was correct for 93% of scaffolds, with errors mostly resulting from short scaffolds.

Taken together, the chromosome-length, small, and tiny scaffolds accounted for 97.3% of the chromosome-length scaffolds of hg38; the remainder was mostly due to repetitive sequences that could not be adequately assembled from short reads. Our method was further validated by obtaining similar results with a draft assembly generated with Pacific Biosciences long reads, which contained longer contigs (16).

Next, we applied our approach to Ae. aegypti, which was previously assembled from Sanger reads (8× coverage) (17). This assembly, AaegL2, contains 1.3 Gb of sequence (contig N50: 83 kb) partitioned among 4756 scaffolds (scaffold N50: 1.5 Mb).

To improve AaegL2, we generated in situ Hi-C data (40× sequence coverage). After setting aside 2222 scaffolds shorter than 10 kb (spanning 1% of the bases in the initial assembly), we used Hi-C data to split, anchor, order, orient, and merge the remaining 2534 scaffolds. Notably, our pipeline identified apparent misjoins in 1422 of these input scaffolds (56%).

The resulting assembly, AaegL4, contained three huge scaffolds (307, 472, and 404 Mb in length) comprising 93.6% of the input sequence, together with an additional 3981 small scaffolds (N50 of 65 kb, maximum of 474 kb) comprising the remainder. The three huge scaffolds correspond to chromosomes 1, 2, and 3 of the Ae. aegypti genome (18) (Table 1 and tables S2 to S7).

We compared our assembly to a genetic map of Ae. aegypti (19). Of the 2006 markers in the genetic map, 1826 markers could be unambiguously mapped in AaegL4. Notably, our assembly agreed with the genetic map for 1822 of these 1826 markers (Fig. 2). All exceptions were due to misjoins in AaegL2 that had not been detected in AaegL4. We also observed close correspondence with a physical map of the Ae. aegypti genome (fig. S12).

Open in a separate window

Fig. 2

Comparison of AaegL4 and CpipJ3 with genetic maps

(A) We compared AaegL4 with a genetic map of Ae. aegypti (19). Our assembly agreed with the genetic map on 1822 out of 1826 markers. The exceptions are due to misjoins in AaegL2 that were not corrected in AaegL4. (B) Similarly, CpipJ3 is in agreement with a genetic map of Cx. quinquefasciatus (21).

Next, we used our approach to create a genome assembly of the mosquito Culex quinquefasciatus, which, like Ae. aegypti, is a disease vector—in this case for West Nile virus, St. Louis encephalitis, and lymphatic filariasis. We generated in situ Hi-C data (~100× sequence coverage) and used them to improve the previous assembly, CpipJ2 (20), obtaining a new assembly, CpipJ3, with three chromosome-length scaffolds that together contain 94% of the sequence in the initial assembly (Table 1 and tables S2 to S7). We validated CpipJ3 by comparing it to existing genetic and physical maps of the Cx. quinquefasciatus genome (20, 21) (Fig. 2 and figs. S13 and S14).

The mosquito Hi-C data and the draft assemblies were generated from different strains. Ideally, the same strain would be used in both cases.

The creation of chromosome-length scaffolds for Ae. aegypti and Cx. quinquefasciatus allowed us to use our Hi-C data to create a Hi-C heatmap showing proximity relationships between chromosomal loci throughout both genomes (22, 23) (Fig. 1 and fig. S15). Notably, the distal ends of the three chromosomes show spatial clustering in both species. Both species also exhibited a second spatial cluster, comprising three loci: one locus from each chromosome, positioned roughly in the middle. This clustering is consistent with the spatial clustering of centromeres, which is known to be present in many organisms. Taken together, the 3D maps are consistent with a spatial arrangement known as the Rabl configuration (24). Our findings also suggest the position of each chromosome’s centromere and thereby partition each mosquito chromosome into two arms.

The assemblies of the Ae. aegypti and Cx. quinquefasciatus genomes allowed us to study genome evolution. We began by examining a whole-genome alignment between the published Anopheles gambiae genome, which is 278 Mb long, and the Ae. aegypti genome, which is ~1.2 Gb long. This analysis identified 1389 large blocks of conserved synteny (fig. S16). Similar results were observed for Cx. quinquefasciatus. Despite extensive rearrangements, we observed correspondence of sequence content among chromosome arms in An. gambiae, Cx. quinquefasciatus, and Ae. aegypti. Specifically, for the vast majority of DNA sequences on a particular chromosome arm in one of the three species, the homologous sequences were all found on a single chromosome arm in the other two species. The only exception is the observation that a single arm in An. gambiae (2R) corresponds to two arms in both Ae. aegypti (1q and 3p) and Cx. quinquefasciatus (1q and 3q). This is consistent with the breakage of this arm in the lineage leading to the shared ancestor of Ae. aegypti and Cx. quinquefasciatus (Fig. 3 and tables S8 to S11). These observations are consistent with cytogenetic analyses (18–20) (figs. S17 and S18).

Open in a separate window

Fig. 3

The content of chromosome arms is strongly conserved across mosquitoes

Here each 100-kb locus in Ae. aegypti is assigned a color. For the other species, each 100-kb locus is assigned a combination of the colors of the corresponding DNA sequences in Ae. aegypti, weighted by length. MYA, million years ago.

Taken together, these results suggest that—with the exception of the breakage event noted above—each chromosome arm in the Aedes, Culex, and Anopheles species descends from a single arm present in their common ancestor about 150 to 200 million years ago. The preference for within-arm rearrangement in mosquitoes is stronger than has been observed in mammals (25).

Notably, the left arm of chromosome 2 in Drosophila melanogaster has a clear counterpart in all three mosquito species. Thus, all four arms derive from a single chromosome arm present in their dipteran ancestor a quarter of a billion years ago (Fig. 3 and fig. S19).

Overall, our results show that incorporating Hi-C data into genome assembly provides a rapid, inexpensive methodology for generating highly accurate de novo assemblies with chromosome-length scaffolds. At present, the total sequencing costs for a 3D de novo assembly are below $10,000 for mammalian genomes and less for smaller genomes (table S12).

It is important to bear in mind that these assemblies still contain errors. For example, although the Hi-C data provide extensive links covering large distances, the current approach is not perfect for local ordering of small adjacent contigs. This might be circumvented by more sophisticated analysis of Hi-C data. Additional data (such as long or paired-end reads) could also improve the results. The ability to rapidly and reliably generate genome assemblies with chromosome-length scaffolds should accelerate genomic analysis of many organisms.