Assembly and diploid architecture of an individual human genome via single-molecule technologies

Contiguity and accuracy of scaffolds relative to hg19

We compared our assembly against the published Allpaths-LG NA12878 assembly, which used short-read sequencing of insert and fosmid libraries⁵. A high-level comparison of the two assemblies, using metrics from refs. 22 and 23, is shown in Supplementary Table 1. Normalizing to hg19, our assembly has a higher contig N50 (886 kb versus 19 kb), scaffold N50 (26 Mb versus 10 Mb) and ‘scaffold accuracy’ (98.7% versus 94.9%), which represents the odds of being correctly connected at a distance of 100 kb. Additionally, fewer hg19 reference bases are missing (14.9% versus 7.6%), and more new assembly sequence was potentially added (58 Mb versus 9 Mb).

However, our sequence identity compared to hg19 is lower (99.7% versus 99.8%); though some of this deviation may be due to detection of true variants or alternative alleles, much of it represents miscalled small indels that result from the higher, indel-based error rate of SMRT sequencing. Such errors can be resolved by mapping short-read data¹ to contigs and using variant calls to correct contigs, leading to sequence identity consistent with the Allpaths-LG assembly (Supplementary Table 2). Using heterozygous SNVs, over 2 Gb of sequence was resolvable into haplotype blocks, with a haplotype N50 of 145 kb (Supplementary Table 3). Last, we measured the structural fidelity of both scaffolds by performing nick-site mapping (or in silico nick-site mapping) relative to each other and hg19 (Supplementary Fig. 3). Fewer chimeras were observed in the V2 scaffolds; moreover, when discrepancies existed between V2 and Allpaths scaffold mapping, the V2 scaffold was 15 times more likely to be consistent with hg19 (Supplementary Results).

Phasing

Phasing was performed using a combination of short-and long-read approaches, enabling long haplotype blocks with low switch error rates and resolving unphased variants from trio-based approaches. SNVs and indels previously identified by deep Illumina sequencing of the NA12878 trio²⁴ represented 2,367,085 heterozygous events (1,925,040 phased by trio analysis), far more than those detected by PacBio sequencing alone (Supplementary Results). The consistency of SMRT sequencing–based phasing with trio results markedly improved when SNV filtering was performed with a modified reference (Supplementary Table 1 and Online Methods). The switch error rate was estimated by measuring concordance of the predicted haplotype blocks with the SNVs labeled by trio-based phasing. After filtering, the estimated switch error rate was reduced to 0.1, which was lower than the estimated switch error rate of 0.9 in HuRef²⁵. There is a trade-off between reducing the switch error rate and eliminating SNVs from analysis (using multiple parameters for both long- and short-read data sets; Supplementary Fig. 4). To increase accuracy, 369,785 SNVs were eliminated from the analysis; however, a similar number of additional variants not amenable to trio-based phasing were resolved (314,630) via the long reads. For all heterozygous TR and SV events, we used local phasing in an attempt to assign events to either the maternal or paternal haplotype. Both alleles were assigned to distinct haplotypes for 9,196 TRs and 3,562 SVs. The assignment approach also serves to assess the accuracy of heterozygous calls. True events should have distinct haplotypes assigned to the two read clusters; here, 97% of predicted heterozygous SVs form consistent haplotypes (Supplementary Fig. 5 and Supplementary Results).

Structural variation

SV calls from PacBio data were generated using the de novo sequence assembly and read-mapping approaches (Online Methods). Some SVs evaluated by locally comparing our de novo sequence assembly to syntenic intervals of hg19 (using tools developed in ref.26) were shared with those previously detected in the CHM1 haploid cell line²⁶ for both insertions (39%) and deletions (12%). Although each genome clearly had many unique variants, they were largely comparable in the magnitude of calls (Supplementary Result and Supplementary Table 4).

SVs were further evaluated by aligning individual raw and error-corrected reads to hg19 (Supplementary Tables 5 and 6 and Supplementary Fig. 6). Short-read calls from tandem duplications, deletions and inversions were evaluated on both short and long insert data using Delly²⁷. These were compared to the PacBio call set and further evaluated by manually inspecting break-point-spanning reads (Supplementary Result and Online Methods). Of the callable short-read predictions, 95% agreed in approximate variant type with the PacBio data (Table 2). Substantially more SVs were predicted in the PacBio data set (Supplementary Table 5), even when considering assembly- mapping or read-mapping approaches separately. For insertions, as expected, the most frequent mobile element insertions corresponded to Alu elements, L1s and SVAs (SINE-VNTR-Alu) (Supplementary Table 4 and Supplementary Fig. 7), elements known to be active in the human genome and expanded in the human lineage^28,29. To estimate the false discovery rate of the mapping-based calls without short-read support, we interrogated events using long-range PCR. For PCR validation, SV predictions were divided into bins on the basis of predicted event size (200–500 bp, 500–1,000 bp, 1,000–1,500 bp, 1,500–5,000 bp). Of the 59 successful PCR reactions, 58 positively supported the predicted event (Supplementary Results and Supplementary Table 7). Additionally, whole-genome data from Tru-seq, an orthogonal long-read platform, were largely consistent with PacBio calls (Supplementary Result).

Tandem repeats

TRs represent an important source of variation that are associated with a broad range of diseases³⁰ but are not easily addressed by NGS technologies. The combination of single-molecule and long-read approaches not only identified large TRs outside the range of short-read approaches but also suggested that many TRs are substantially larger than indicated in hg19 (Fig. 2). Some of these differences may be explained by allelic variation at a given locus (Fig. 2b), but there is clearly a systematic underrepresentation of repeat signals in the reference (Fig. 2a), with certain regions showing increased variability (Supplementary Results and Supplementary Fig. 8). Small events were compared to short-read predictions by RepeatSeq³¹ and showed >90% concordance (Supplementary Fig. 9). As events become larger, they seem to be more consistent with the reference. However, this is most likely due to the exponential read-length distribution observed in SMRT sequencing. As fewer reads reliably span larger TR intervals, high-quality alignments are more likely to be observed when consistent with the reference. Thus, very large TRs (Fig. 2c) can only be directly examined using genome mapping data. An example of this is a TR within the LPA gene on chromosome 6q26 (the ~5.6-kb ‘kringle-IV’ type 2 (KIV2)-like domain; Fig. 2c); long molecules spanning over 100 kb are needed to reliably span the TR. Correctly identifying its multiplicity is particularly relevant as LPA size has been associated with plasma lipoprotein level and risk of cerebrovascular and cardiovascular diseases in the human population³².

Large variation via assembly and scaffolding

The scaffolding results largely validate the layout predicted by hg19 (Fig. 1); however, large structural variation events are observed in the hybrid scaffold (Supplementary Table 8), and a number of large SVs were directly identified by genome mapping (Supplementary Table 9). To validate these events, we compared the genome maps, and the raw molecules used to construct them, to hg19. For example, a 206.6-kb insertion seems to be a large TR expansion (Fig. 3a). A number of raw molecules spanning the event support the BioNano assembly, whereas no spanning molecules confirm the smaller reference allele. The observed difference could be due to variability in the population or artificial compression from traditional assembly approaches. In another example, a 577.3-kb inversion spans previously unresolved regions in hg19 (Fig. 3b), suggesting potential misassembly during BAC tiling layout. This is supported by higher concordance with the updated GRCh38 (hg38) assembly. Other large events (Supplementary Fig. 10) show our assembly to be more consistent with hg38, suggesting that the hg38 assembly has fixed errors in hg19 or is more representative of dominant haplotypes. Yet, despite these improvements, gaps still persist in hg38. Our sequence assembly resolves 28 previously defined ‘interstitial gap’ intervals²⁶, yielding 34 kb of assembled sequence that spans 621 kb in hg38 (Supplementary Table 10). The resulting gap sequence is enriched for simple repeats (Supplementary Table 10), consistent with previous long-read gap closure results in hg19 (ref.26).

Complexity of variant sequences

As mentioned earlier, spanning long reads enable direct observation of breakpoints and inserted sequences. This allows one to distinguish mobile element insertions that only contain the repeat element from those that contain other inserted sequences (Supplementary Fig. 7a). Although many of these events may be duplications not derived from mobile element insertion, some of these intervals are the result of SVA^33,34 or L1 (ref.35) DNA transduction: for example, a 5′ truncated SVA element mediating the 3′ transduction of a proximal Alu sequence (Supplementary Fig. 11). SMRT sequencing long reads can also be used to distinguish subtle SV insertions within TR intervals³⁶ (Supplementary Table 11). A common feature of these internal SVs appears to be distinct repeat substructures within the putative inserted SV. In one example, the canonical reference repeat of AGG is interrupted by three distinct (but related) repeat subpatterns (Fig. 2d).

Spanning long reads also elucidate complex rearrangements typically missed by conventional NGS in which multiple events are located together. The assembly-based approach identified 4.2% of events as complex, and whereas Delly short-read predictions were largely confirmed, substantially greater complexity was observed in the variants (Table 2), with 3.4% showing added complexity. Inversions seem to be particularly enriched for complexity (55%; Supplementary Table 6), a feature we are exploring further in the context of a large population cohort (T.R., M.H.-Y.F., A.M.S., A.B. and J.O.K., unpublished data). We find predicted inversions located together with insertions (Fig. 4a), deletions (Fig. 4a,b) and duplications (Fig. 4c). A number of inversions also showed overlapping boundaries (for example, inverted repeat structure at the inversion boundaries), making it challenging to resolve precise breakpoints (Supplementary Table 6). Another arrangement frequently observed in the data is the insertion or deletion of proximally located sequences, which we refer to as proximal duplicated or deleted substrings. These appear in both forward (Fig. 4d) and inverted (Fig. 4e) orientation, and highly similar substrings can excise and insert multiple times within the same genomic interval (Fig. 4e). These events are particularly challenging to detect using short reads and are often mischaracterized as tandem duplications or inversions (Table 2).

High-molecular-weight DNA extraction, DNA labeling and data collection

NA12878 cells were washed with PBS, and the final cell pellet was resuspended in cell suspension buffer (CHEF Genomic DNA Plug Kit). Cells were embedded in a thin LMP agarose layer and lysed, protease treated and washed. Purified DNA embedded in a thin agarose layer was labeled following the IrysPrep Reagent Kit protocol (BioNano Genomics). Briefly, DNA was digested with Nt.BspQI nicking endonuclease (New England BioLabs) for 2 h at 37 °C. Nicked DNA was then incubated for 1 h at 50 °C with fluorescently labeled dUTP and Taq Polymerase (New England BioLabs). Taq ligase (New England BioLabs) was used in the presence of dNTPs for ligation of nicks. Recovered DNA was counterstained with YOYO-1 (Life Technologies).

Labeled and counterstained DNA samples were loaded into IrysChips (BioNano Genomics) and run on the Irys (BioNano Genomics) imaging instrument. Data were collected for each sample until >50-fold coverage of long molecules (>180 kb) was achieved. The IrysView (BioNano Genomics) software package was used to detect individual linearized DNA molecules using the YOYO-1 counterstain and to determine the localization of labeled nick sites along each DNA molecule. Sets of single-molecule maps for each sample were then used to build a full genome assembly.

De novo assembly of genome maps

De novo assembly of single molecules was accomplished using a custom BioNano assembler software program based on an Overlap-Layout-Consensus paradigm^51–53. First, we started with pairwise comparison of all molecules longer than 180 kb and nine labels to find all overlaps with P < 1 × 10⁻¹⁰, then we constructed a draft consensus map on the basis of these overlaps. The draft map was further refined by mapping single molecules to it and recalculating the label positions. Next the consensus maps were extended by aligning overhanging molecules to the consensus maps and calculating a consensus in the extended regions. Finally, the consensus maps were compared and merged where patterns matched with P < 10⁻¹⁵. The process of extension and merge was repeated five times before a final refinement was applied to ‘finish’ all genome maps. The result of this assembly is a genome map set entirely independent of any known reference or external data.

Error correction and assembly

Error correction of all reads was performed using Falcon, following the general principles proposed in ref.15 (Supplementary Note 1). In short, all long reads greater than 3 kb were first aligned to one another using Blasr⁵⁴. These reads were then grouped together by selecting the top alignments (using a coverage cutoff of 40). A consensus was formed for each read; the resulting read was trimmed at the ends to eliminate potential chimeras and low-quality sequence (here we require at least 5× coverage of a given base). Error-corrected reads were passed to the Celera Assembler to form contigs (Supplementary Note 1). The resulting raw contigs were passed back to the Quiver pipeline (SMRTAnalysis v.2.2.0) on the subset of raw reads corresponding to the newer chemistry to provide the final, high-quality sequences. These sequences were then merged with genome mapping scaffolds to produce our initial V1 assemblies. For scaffolding purposes, an alternative long-read assembly was generated using a modified form of the Falcon pipeline, which yielded a more aggressive assembly with a 2.1 Mb N50 (Supplementary Fig. 12 and Supplementary Table 1). Final assemblies were cleaned to remove contaminants (Supplementary Note 1).

Short-read mapping and variant calling and correction

Short reads from ref.1 were mapped using BWA-MEM⁵⁵ (version 0.7.12-r1039) with default parameters. Variants were called using Freebayes⁵⁶ (version 0.9.18-3-gb72a21b) with default parameters but were filtered for variants with Q50 or higher.

PacBio structural variation

Events were called using the methodology from ref.26, PBHoney⁵⁷ and a custom pipeline. For comparison to the haploid CHM1 assembled breakpoint data set more directly, only events spanned by the de novo assembly were considered (though some alternative haplotypes persisted in the assembly). A brief description of the custom pipeline follows (see Supplementary Note 2 for detailed overview). Reads were first aligned to the reference using Blasr via an iterative process. First, the full read was mapped, maintaining the top ten highest-scoring alignments relative to the reference. Next, unmapped portions of each read were extracted from the input read set and remapped to the reference to identify potential highly divergent rearrangements that were missed in the initial mapping step. The top alignment for each query was then passed into a three-state HMM to identify potential insertions and deletions contained within a single long-read alignment. The HMM is needed because of a lack of affine gap alignment⁵⁸ in the initial Blasr results (1.3.1) and to reduce false positive calls in high-degeneracy raw PacBio reads, both of which lead to sporadic alignment of query bases within a true deletion region (or reference bases within a true insertion region). Note that using an updated version of Blasr with affine gap alignment (Supplementary Note 2) resulted in improved specificity at improved, or similar sensitivities, for most read types and event categories (Supplementary Table 5).

For more complex or larger rearrangements, a secondary step was performed in which a directed alignment graph is created. Alignments (nodes) were ordered relative to their position on the query; a directed edge was drawn between alignments a and b if the end of alignment a preceded the end of alignment b. A source, s, node and a sink, k, node are created for each query, and two edges, (s, u) and (u, k), are added to each alignment node (Supplementary Fig. 13). A simple dynamic programming algorithm determines the highest-scoring path from source to sink, where overlapping alignments are rescaled to eliminate double counting of overlapping intervals. This highest-scoring path is returned if it indicates a nonreference alignment path for the query. Although this step was not rate-limiting, it is shown in ref.59 that sparse dynamic programming approaches can yield O(n log n) runtimes, and these approaches have been applied in the context of detecting rearrangements^60,61. The resulting set of individual read calls is then clustered by event type across the entire data set to yield the final set of predicted SVs. Both error-corrected and raw PacBio reads were processed via the same protocol (Supplementary Note 2).

Detection of structural variants with CLRs differs from paired reads in that the detection of SVs often implies complete resolution of the spanned event. However, given the potential for chimeras, and specifically the known issue of inverted tandem repeat–like chimeras due to missed adaptor sequence, singleton events are not sufficient to accurately resolve a sequence¹³. Therefore, to provide and validate predictions, we performed consensus calling across all putative events using partial order alignments of all ‘event’-containing reads⁶². For insertions, these consensus sequences were then scanned through Dfam HMMs to identify putative repeats using “nhmmscan” with default parameters⁶³.

Inversions were broken up into several distinct categories for custom analysis: (i) spanned inversions with a single contiguous alignment, (ii) spanned inversions with multiple subalignments and (iii) inversions in which only a single breakpoint is observed (Supplementary Fig. 14). Additionally, all inversion calls were passed through an additional step that used custom dot-plotting script followed by manual analysis from spanning reads to reduce false positives.

Genome map structural variation

The structural variation algorithm begins by aligning genome maps to the reference (hg19). The alignment algorithm uses a Smith-Waterman style dynamic programming algorithm where the units of comparison are distance intervals between detected label sites (as opposed to base pairs). Intervals are compared using maximum likelihood and a noise model designed for BioNano data. Both orientations are aligned separately for each map, and the best scoring alignment for each genome map, over the entire reference, is recorded. The algorithm scores each interval: positive scores are given when the interval in the reference and genome maps agree according to the model, and negative scores are given when they do not agree. If there are outliers in the alignment, or if one side of the genome map does not align (i.e., all intervals have a negative score), the map is split at the outlier positions. Here, we set the outlier cutoff to be 10⁻⁶, which represents the probability of the interval being similar by random chance. The split subintervals are then realigned to the entire reference (with each piece again permitted to align in either orientation). In the case of a large insertion or deletion (>3 kb), the split portions of the map on either side of the SV will align next to one other. In the case of inversions, the split portions will again align next to each other, but in opposite orientations.

Tandem repeats

The tandem repeat detection pipeline uses PacBio long reads, alignment of reads to hg19 (via Blasr) and a reference TR table (available from University of California–Santa Cruz) as input. Only the top scoring alignment was used for each read, and only reads which had at least 100 bp of sequence anchoring upstream and downstream of the TR were considered for analysis. The method is based on the work in ref.36. A summary proceeds as follows. First, a three-stage dynamic programming algorithm step more robustly identifies the TR region in the long read and the putative boundaries. The TR region is then passed into the pairHMM-based method to provide a more robust and probabilistic estimate of TR multiplicity using the appropriate error profile. The processed pairHMM output is used for clustering. The objective of the clustering routine is to call the allele based on the estimated number of TRs for all the reads that span a particular TR event. We keep track of several key features on the clustering, the binomial probability of the clustering split, the minimum number of reads in each TR allele, the total number of reads given as input (to distinguish potential copy number abnormalities that could lead to spurious calls) and the c-separation, which specifies the separation between the means of two clusters. We return this information along with the cluster means and s.d. for each cluster. The sequences identified within each cluster of a TR event are used as an input for the partial order alignment (POA) consensus generation step. SVs internal to TRs are obtained by traversing the raw pairHMM output to find intervals in the query that are of low quality relative to the consensus TR element. TRs were filtered to exclude those that are segmental duplications as well as those that are contained within another repeat (if two TRs overlap but are not fully contained, we excluded the region). In a situation where two alternative TRs were present and both were completely contained within one another, we selected the TR with larger period.

Phasing

SNVs and indels identified by high-depth Illumina sequencing of the NA12878 trio²⁴ were used as a starting point for phasing (Broad Institute, GATK 2.5; Haplotype Caller version 2.7). Long reads were phased relative to hg19 using HapCut⁶⁴ version 0.7, which implements a graph-based optimization heuristic and has been previously applied to phase HuRef. PacBio reads have a well-known reference bias in which SNVs are more likely to be called as the reference owing to alignment artifacts created by the high insertion and deletion error rates of raw reads⁶⁵. To mitigate this process, we first created a ‘variant-free’ version of hg19 in which all known variant bases in the reference were converted to ‘Ns’. Reads were then mapped to this and assessed for consistency with short-read trio predicted variants in NA12878. Short reads were also mapped to assembled PacBio contigs, and heterozygous positions were phased using the same approach (after correction of high-quality homozygous variants). Variants were included in phasing if they were covered by at least ten PacBio reads with at least 25% of reads supporting both alleles and 20 Illumina reads with at least 25% of reads supporting both alleles.

To phase tandem repeats and structural variants, we extracted reads spanning the event of interest. Reads were separated into two sets on the basis of which allele they supported (in the case of ambiguous alignment, the read was discarded from analysis). SNVs immediately upstream or downstream of a putative event were used to assess phase by performing variant free alignment, as described above. For each allele, a consensus SNV was called at each position, and a label (maternal or paternal) was placed on the allele on the basis of aforementioned trio calls. In many cases, insufficient SNVs were available in the flanking region, and these regions were listed as ambiguous and eliminated from haplotype consistency analysis.

We attempted to phase each of the predicted tandem repeats and structural variants by placing them in the context of the previous trio-based phasing. Each tandem repeat and structural variant was evaluated to produce ‘high-confidence’ heterozygous calls. The set of reads corresponding to each allele of a high-quality call was retrieved and used to locally determine the haplotype. In short, each read set for an allele was evaluated at all known SNV or indel locations proximal to the SV, where at least two reads covered the event using a POA alignment. As before, the reference SNV or indel position was eliminated from the reference sequence to eliminate reference bias. A consensus haplotype was then established for each allele; if the consensus haplotypes were consistent with trio-based phasing, we then assigned each allele to its corresponding haplotype.