De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds

See allHide authors and affiliations

Science  23 Mar 2017:
DOI: 10.1126/science.aal3327


The Zika outbreak, spread by the Aedes aegypti mosquito, highlights the need to create high-quality assemblies of large genomes in a rapid and cost-effective fashion. Here, we combine Hi-C data with existing draft assemblies to generate chromosome-length scaffolds. We validate this method by assembling a human genome, de novo, from short reads alone (67X coverage). We then combine our method with draft sequences to create genome assemblies of the mosquito disease vectors Aedes aegypti and Culex quinquefasciatus, each consisting of three scaffolds corresponding to the three chromosomes in each species. These assemblies indicate that virtually all genomic rearrangements among these species occur within, rather than between, chromosome arms. The genome assembly procedure we describe is fast, inexpensive, accurate, and can be applied to many species.

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 sequence, 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. This data is 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, due to repetitive sequence 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-15 megabases (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 distance, in base pairs, between a pair of loci. For instance, loci separated by 10kb in the human genome form contacts 8 times more often than those at a distance of 100kb. In absolute terms, a typical distribution of Hi-C contacts from a given locus is 15% to loci within 10kb; 15% to loci 10kb – 100kb away; 18% to loci 100kb – 1Mb away; 13% to loci 1Mb – 10Mb away; 16% to loci 10Mb – 100Mb away; 2% to loci on the same chromosome, but more than 100Mb 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 datasets. For the maps reported in (4, 5), summing the span of each individual contact reveals that 1X of sequence coverage of the target genome translates, on average, into 23,000X 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 (68). 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 misassemblies 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 one-dimensional genome. Finally, we merge contigs and scaffolds that correspond to overlapping regions of the genome by identifying pairs of scaffolds exhibiting both strong sequence homology as well as strong similarity in long-range contact pattern (Fig. 1 and fig. S1).

Fig. 1 Starting with a draft assembly, we use Hi-C data to correct misjoins, scaffold, and merge overlaps, thereby generating an assembly of the Aedes aegypti mosquito genome with chromosome-length scaffolds.

Here, we show contact matrices generated by aligning a Hi-C dataset to both the AaegL2 assembly (17) 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 co-locate in the nucleus. The loci corresponding to each row and column are illustrated using chromograms. On the left, the chromogram depicts the 3 linkage groups (Lnk1, Lnk2, Lnk3, or Unassigned) reported in AaegL2; on the right, it depicts the 3 chromosome-length scaffolds in AaegL4 (chr1, chr2, chr3). To create the chromogram, each AaegL4 arm is assigned 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 (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 (middle row). Finally, segments exhibiting a similar 3D signal are examined for evidence of overlapping sequence (green rectangle) and merged (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 (67X coverage). We created a draft assembly from 250bp paired-end reads (60X coverage, generated by Illumina sequencing with a PCR-free protocol, downloaded from Sequence Read Archive (SRX297987); assembled with DISCOVAR de novo (13)). This assembly, dubbed Hs1, comprises 2.82 Gb of sequence (contig N50 length: 103kb) partitioned among 73,770 scaffolds (scaffold N50: 126kb; Table 1).

Table 1 Assembly statistics for the Hs2-HiC, AaegL2, and CpipJ3 assemblies.

We did not attempt to further assemble tiny scaffolds contained in each draft assembly. The other scaffolds in each draft were assembled using Hi-C to create huge, chromosome-length scaffolds, and additional small scaffolds.

View this table:

We then used in situ Hi-C data (6.7X sequence coverage) to improve Hs1. We set aside the tiny scaffolds (43,231 scaffolds shorter than 15kb, whose N50 length is 6.1kb). Together, these contain 5.4% of sequenced bases in Hs1. Due to 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.8Mb to 225.2Mb) containing 99.5% of the total sequence, together with an additional 811 small scaffolds (N50 length of 30kb; maximum length of 231kb) 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 provides less information to resolve the fine-structure order of short scaffolds. However, the agreement was 99% for scaffolds of length at least 120kb. Similarly, the orientation was correct for 93% of scaffolds, with errors mostly due to 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 using 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 (8X coverage) (17). This assembly, ‘AaegL2’, contains 1.3Gb of sequence (contig N50: 83kb) partitioned among 4756 scaffolds (scaffold N50: 1.5Mb).

To improve AaegL2, we generated in situ Hi-C data (40X sequence coverage). After setting aside 2222 scaffolds shorter than 10kb (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 (307Mb, 472Mb, and 404Mb in length) comprising 93.6% of the input sequence, together with an additional 3981 small scaffolds (N50 of 65kb, maximum of 474kb) comprising the remainder. The three huge scaffolds correspond to chromosomes 1, 2, and 3 of the Ae. aegypti genome (18) (Table 1 and tables S1 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. Strikingly, 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).

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 (100X sequence coverage) and used it 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).

We note that the mosquito Hi-C data and the draft assemblies were generated from different strains, although 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). Strikingly, 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 278Mb long, and Ae. Aegypti, which is ~1.3Gb 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 (1820) (figs. S17 and S18).

Fig. 3 The content of chromosome arms is strongly conserved across mosquitos.

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

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 approximately 150-200 million years ago. The preference for within-arm rearrangement in mosquitos is stronger than has been observed in mammals (25).

Interestingly, 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 sequencing costs 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, while the Hi-C data provides 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.

Supplementary Materials

Materials and Methods

Figs. S1 to S21

Tables S1 to S20

References (26-41)

References and Notes

  1. Acknowledgments: This work was supported by a Center for Theoretical Biological Physics Postdoctoral Fellowship to O.D., an NIH New Innovator Award (1DP2OD008540-01), an NSF Physics Frontier Center (PHY-1427654, Center for Theoretical Biological Physics), the NHGRI Center for Excellence for Genomic Sciences (HG006193), the Welch Foundation (Q-1866), an NVIDIA Research Center Award, an IBM University Challenge Award, a Google Research Award, a Cancer Prevention Research Institute of Texas Scholar Award (R1304), a McNair Medical Institute Scholar Award, and the President’s Early Career Award in Science and Engineering to E.L.A. We thank Chad Nusbaum, Leslie Vosshall, Ben Matthews, and Raul Andino for their thoughtful comments on this manuscript. All datasets reported in this paper are available at the Gene Expression Omnibus (GEO), series accession number GSE95797. Code and details to reproduce the assemblies are available at O.D., S.S.B., A.D.O., S.K.N., N.C.D., E.S.L., A.P.A., and E.L.A. are inventors on U.S. provisional patent application 62/347,605, filed June 8, 2016, by the Baylor College of Medicine and the Broad Institute, relating to the assembly method in this manuscript.
View Abstract

Navigate This Article