A sequencer coming of age: De novo genome assembly using MinION reads

2.1 Segmentation and event detection

Before basecalling takes place, current basecalling tools require the continuous signal to be divided into discrete sections, called events, each supposedly representing one combination of bases. Furthermore, the poly-T lead preceding the sequence needs to be removed, and, if 2D chemistry is used, the signal derived from the template strand needs to be separated from that of the hairpin and the complement strand. The latter process is commonly referred to as segmentation.

For locally basecalled reads, both segmentation and subsequent event detection are performed by the MinKNOW software provided by ONT. The software identifies the hairpin by a characteristic double current spike, while a long poly-T signal denotes the lead. For event detection, MinKNOW was reported to calculate a simple t-statistic between sliding adjacent windows of set size. Peaks in the t-statistic above a certain threshold are then assumed to signify the border between two segments. Optionally, the raw current signal can also be reported for further analysis.

2.2 Basecallers

Several dedicated basecalling tools are already available to MinION users. In this section, the underlying principles and implementation of these tools are explored, along with their reported subsequent strengths and weaknesses. It should be noted that an all-inclusive fair comparison, in which all basecallers are optimized for the latest MinION chemistry and then run on the same dataset, is currently lacking. Any comparison in this chapter is therefore based on accuracy reports made by the authors of the open-source basecallers in their publications (i.e. between the Metrichor basecallers and their own).

Metrichor basecallers Metrichor, a spin-off company of ONT and its main developer of proprietary analysis software, maintains a range of basecallers that have remained the go-to option for most MinION users. Initially, the Metrichor basecallers relied on hidden Markov models (HMM) to find the biological sequence underlying the segmented signal. As of early 2016, this model was replaced by a more accurate recurrent neural network (RNN)-implementation. Currently several RNN-based versions exist under different names; Albacore, Nanonet and the MinKNOW integrated basecaller. Albacore and the MinKNOW version are stable versions intended for regular MinION users. As of Albacore v1.0.1 and MinKNOW v1.6.11, both were reported to improve homopolymer calling through the implementation of a transducer, a feature that was transferred from the Scrappie basecaller (see below). Nanonet is an unsupported version under continuous development and makes use of the CURRENNT library to implement its RNN³³. A cloud-based version was previously integrated in the EPI2ME platform, but this service has been discontinued. Due to the proprietary nature of the software, the source code of Albacore and the MinKNOW basecaller is currently only open to developer users, while that of Nanonet is openly available. The community manual describes the general procedure followed by all Metrichor basecallers as follows. Long reads are broken up into partially overlapping stretches, referred to as chunks. This keeps the Viterbi-decoding process, required after the RNN has determined k-mer probabilities, tractable. The RNN then assigns a likelihood for every possible five-mer to each event in the chunks.

Metrichor basecallers were reported to use a bidirectional long-short term memory (BLSTM)-RNN. This architecture allows the RNN to overcome the vanishing gradient problem, encountered often during the training of deep neural networks, and the usage of multiple previous events in the sequenceⁱⁱ.

The RNN assigns probabilities for each k-mer to each event, which are used to estimate basecalling quality scores (i.e. q-scores). Five-mers that would entail a forbidden move, considering the likelihood of surrounding k-mers, are removed at this point. The maximum likelihood sequence is then extracted from the k-mer likelihoods through Viterbi decoding. Finally, basecalled chunks are merged again.

The accuracy of the Metrichor basecallers is currently considered to be the highest of all available basecallers. Before the switch to RNNs, an international performance analysis estimated the accuracy at 88% for 2D basecalling⁵. David et al. reported an accuracy of approximately 68%⁶ for 1D basecalling on human data, while Boza et al. reported 77% for 1D and 87% for 2D on Escherichia coli and Klebsiella pneumoniae genomes. After moving to an RNN-implementation (as well as numerous improvements in used chemistry and hardware) this percentage has risen to over 95% (see here for a publicly available dataset supporting this), even before the mentioned correction of homopolymer regions in Albacore and MinKNOW was applied.

Scrappie A basecaller by Metrichor and platform for ongoing development, Scrappie was the first basecaller reported to specifically address homopolymer basecalling, through the implementation of what is referred to as a transducer. It was written in C and its source code is publicly available. In previous MinION sequencing efforts, correctly determining the length of homopolymers proved challenging, as the separation between states in a homopolymer stretch was not clear^20,21,35. Scrappie can process the raw current signal, rather than the event-detected data.

As an initial performance assessment, an early-development version of Scrappie was run in parallel to other Metrichor basecallers on reads of the human chromosome 20. It was found that Scrappie indeed called homopolymer stretches more accurately than the EPI2ME basecaller and Nanonet; homopolymeric stretches of up to 16 bases were called accurately (although a slight over estimation of repeat units was seen between 5 and 15 bases), after which the accuracy fell. In contrast, the EPI2ME basecaller and Nanonet failed to call stretches longer than five bases. The effect was also shown in an analysis of 5-mer counts; of the tested basecallers, Scrappie stayed closest to the reference genome and did not under-represent homopolymeric 5-mers³. Although these results are promising, a thorough assessment of Scrappie in a later stage of development is required to evaluate the true potential of the transducer-based homopolymer calling.

Nanocall Nanocall⁶ was the first open-source basecaller for the MinION offered as an alternative to the proprietary Metrichor software. It was written in C++. Nanocall accepts the segmented signal from MinKNOW and assigns k-mers to the events using a HMM.

An HMM assumes that, in a sequence of events, the probability of the next event being of a certain nature - its state - is dependent only on the current state. This makes sense for the basecalling of subsequent k-mers, as the prefix of the next k-mer that is called should be the suffix of the previous. The "hidden" property of HMMs refers to the fact that the sequence of states - the path - is not known. Instead, it has to be inferred with some probability from a derived signal. In the case of a basecaller, this signal consists of electrical current levels which are derived of a sequence of k-mers. The goal of the HMM is to find the path that is most likely to have emitted the given signal (for a more thorough introduction to HMMs, see 36).

As described by David et al.⁶, the HMM underlying Nanocall takes stay, step and skip scenarios into account at each event transition (Figure 6). Although skips of more than one base may occur in a sequence, Nanocall only considers one-base skips to cut down on computational requirements, thus reducing the number of neighboring states to 21 (i.e. stay, four different steps and sixteen different one-base skips). During an initial training phase, state transition probabilities for stays and skips of any number of bases, p_stay and p_skip respectively, are estimated using the Baum-Welch algorithm (p_step is derived as 1 − p_stay − p_skip)³⁷. The transition probability for exactly one base p_skip1 is then derived as p_skip1 = p_skip /(1 + p_skip). Note that all sixteen one-base skips have the same transition probability (i.e. $\frac{1}{16}$ · p_skip1), as do the four step scenarios ($\frac{1}{4}$ · p_step). To calculate emission probabilities, Nanocall relies on ONT’s pore models, which provide parameters for current level distributions per k-mer, which can be fitted to individual reads. The most probable sequence is obtained through Viterbi decoding.

Figure 6. Example of the most likely path of k-mers, derived by Nanocall from event data by a hidden Markov model (HMM).

The path is chosen by maximizing the product of transition probabilities p_stay, p_step and p_skip1 (modified by fractions in the case of skips and steps, as there are multiple states to which these transitions could lead) and emission probabilities of k-mers. The data used in the construction of this figure is mock data.

Performance of Nanocall at the moment of its publication, in March of 2016, was reportedly on par with that of 1D analysis by the Metrichor basecallers, or around 68% accuracy⁶. However, unlike the Metrichor basecallers, Nanocall does not support 2D analysis, which allowed for a much higher accuracy in the comparison by David et al. (around 85%)⁶. Since the Metrichor basecallers switched to a recurrent neural network-approach, its 1D accuracy proved significantly higher than that of Nanocall as well. Moreover, Nanocall is incapable of calling homopolymer stretches longer than its k-mer length, while current Metrichor basecallers are able to do so. Therefore, as the authors currently state on the Nanocall Github, it should be seen as a platform for testing out novel ideas, rather than a prime choice for basecalling MinION data.

DeepNano Before Metrichor made its own switch from HMM- to RNN-basecalling, the open-source basecaller DeepNano⁷ already implemented a form of RNN-basecalling, booking a significant improvement in accuracy with respect to the then-current Metrichor version (late 2015). DeepNano was written in Python, using the Theano library³⁸.

Rather than LSTM-nodes, as currently used in Metrichor basecallers, DeepNano implements gated recurrent units (GRUs)³⁹. Although the design of GRUs is simpler than that of LSTM units, performance assessments show that GRUs are able to achieve similar or slightly better accuracy in different learning tasks^40,41.

DeepNano’s RNN was trained separately for (hairpin-connected) 1D basecalling and true 2D basecalling. As the exact alignment to the reference may be uncertain, the algorithm alternately trains the RNN weights for 100 iterations and then re-aligns the reads to the reference. At the moment of its release in March of 2016, DeepNano was reported by its authors to make more accurate calls than the Metrichor agent for both 1D calling (approximately 70% versus 77%) and 2D reads (approximately 88% versus 87%). However, as the Metrichor agent has since then switched to an RNN-approach as well, these comparisons are outdated.

2.3 Assemblers

In the assembly stage, basecalled sequences are combined to reconstruct continuous parts of the genome as accurately and completely as possible. Assembly of nanopore reads requires a different approach than that of SGS reads; as nanopore reads are longer, finding a correct overlap should be easier, yet they are more error-prone, which increases the uncertainty of overlaps. Because of these differences, a return of interest in overlap layout consensus (OLC) algorithms - which were at the peak of their popularity in the era of Sanger sequencing - is seen. De-Bruijn graph (DBG) assemblers, the more popular choice for SGS reads, were reported to return lower quality assemblies of nanopore reads than OLC-based methods, but proved faster in some cases⁴². A more detailed comparison between OLC and DBG assemblers, along with implementations of both approaches is discussed in this section. Software using traditional greedy extension algorithms (e.g. SSAKE) was found to perform decidedly less well in a de novo assembly setting, both in terms of assembly quality and required computational resources⁴², and is therefore omitted here.

PBcR & Canu Originally developed for the first human genome draft, the Celera assembly pipeline⁴³ and its extensions [Koren et al., 2017²⁹, Miller et al., 2008⁴⁴, Zimin, 2013⁴⁵, i.a.] have remained a popular choice in a growing landscape of OLC assemblers. Briefly, the Celera assembler uses read overlaps to find contigs of which the structure can unambiguously be derived from overlap information, referred to as unitigs. It then separates repetitive unitigs from unique ones and attempts to orient the unique unitigs with respect to each-other. Where possible, gaps between unique unitigs are filled with repetitive elements. As a high read error rate is detrimental to the quality of the assembly⁴⁶, two different modifications to the pipeline are available. The PacBio corrected Reads (PBcR) algorithm, originally developed for the correction of PacBio reads suffering from similar error rates, uses accurate short reads mapped with high confidence to the long reads to correct errors. The assembly then proceeds as is usually done by Celera⁴⁷. PBcR’s successor, Canu²⁹, provides a solution that is more integrated with Celera and does not require short accurate reads. The pipeline includes three stages; correction, trimming and assembly. Overlaps are found using the efficient MHAP²⁸, which hashes k-mers using different hash functions and for each hash function stores the smallest integer to which a k-mer of the sequence is hashed. Comparing the hashed k-mers per read results in initial overlap hits, which are then used to perform error correction by consensus seeking. By selecting overlaps for correction on quality, but limiting the number of overlaps a read can contribute to, Canu attempts to prevent masking of true repeat variants. Shorter reads are used at this stage to improve accuracy of longer reads. In the trimming step, overlaps are recalculated to locate and filter out regions of low coverage and high error. Reads are overlapped two more times to correct specific types of errors (i.e. missed hairpin sections/adapters, chimeric reads) and adjust error rate per overlap, before the actual assembly phase starts. With adjustments to account for erroneous alignments and residual errors, assembly essentially follows the same procedure as CABOG, another Celera-based pipeline⁴⁴.

Due to its thorough yet relatively efficient correction steps, Canu is significantly more accurate than both its predecessor Celera/PBcR and Minimap/Miniasm²⁵. In a benchmarking effort on Enterobacter kobei reads, it produced an assembly of higher contiguity, with fewer indels and mismatches. This is in line with the accuracy assessment by its authors²⁹.

Minimap & Miniasm In terms of speed and computational efficiency, the OLC-based pipeline consisting of Minimap and Miniasm³² has a definite advantage over other existing tools²⁵. This efficiency was reached through the omission of the consensus step and the use of minimizers. Much like the output of the MinHash algorithm used by Canu²⁹, a minimizer is a memory-efficient hashed representation of a sequence. Minimap computes the set of minimizers of a sequence, the “sketch”, by finding the k-mers represented by the smallest hash value within a certain window size of each position of the sequence (Figure 7). The inverse of each k-mer is also considered. Decreasing the window size will increase the returned number of minimizers and allow for more accurate alignment at the cost of increased computational requirements. Minimap then performs all-versus-all mapping by identifying hits between minimizers of different sequences. The found overlaps are passed on to Miniasm, which constructs an assembly graph. First, artefacts are trimmed from each read by removing all sequences outside the longest stretch with more than 3 times coverage. Then reads contained within other reads are removed and small bubbles, less than 50 kb in length, are popped (i.e. a consensus is taken in cases where paths split and later join up again). Finally, sequences can be extracted from stretches of the graph without multi-edges to form unitigs. The error rate at this point is practically the same as that in the raw reads, emphasizing that correct basecalling is essential for the eventual quality of the assembly. The graphical fragment assembly (GFA) output format of Miniasm conveniently allows both graphing of the uncorrected assembly and addition of consensus error correction tools, such as Nanopolish or Racon, to the pipeline, though the latter increases walltime and computational cost severely²⁵.

Figure 7. The minimizer sketch procedure, demonstrated on a 15-base read, using a window size of 5 and a k-mer value of 3.

(A) To extract one minimizer m, all k-mers and their complements are listed and hashed to a distinct integer using hash function φ. The smallest hash value is then stored with its starting position and a binary value denoting the strand on which it was found (0 for forward and 1 for reverse). If a multiple of the smallest hash value is found within a window, all are stored. (B) The minimizers of all windows in the read are collected and stored as the sketch, M. As M is a set, double minimizers are removed.

In March of 2016, the authors of Minimap and Miniasm reported assembly of MinION reads of an E. coli genome in a single contig. In May of the same year, Judge et al. assembled an Enterobacter kobei genome in 16 contigs with an N50 of 662 kbase in two minutes, while the next fastest assembler (Canu) took two hours. However their benchmark showed that the omission of an error correction step caused the eventual assembly quality of E. kobei to be too low to properly assess by the QUAST analysis tool²⁵. The inclusion of the Nanopolish consensus correction tool after Miniasm improved the quality of the assembly enough to perform a reliable quality analysis. However, the pipeline was still under-performing compared to other tools.

TULIP As more reads are required to cover larger genomes, and as the time required for all-vs-all overlapping increases quadratically with an increasing number of reads, it follows that the overlap step of OLC assemblers may take unfeasibly long for very large genomes. To tackle this issue, The Uncorrected Long read Integration Process (TULIP) takes a different approach to read overlapping²⁰. Instead of all-vs-all alignment, short seed sequences are selected, which the assembler then attempts to align with long reads. This drastically cuts down the overlapping complexity and makes efficient use of long reads to cover long stretches of the genome between the seed regions. The resulting graph represents seeds as vertices and the connecting reads as edges. In a graph cleaning step, vertices with multiple in- or outgoing edges are revisited. Spurious and superfluous edges are removed aggressively, thus producing a linear graph. Note that, as the name implies, TULIP does not perform basecalling error correction.

The success of assembly using TULIP highly depends on proper seed selection. To avoid spurious connections between reads, the seeds need to be sufficiently unique in the genome and contain few sequencing errors. If available, SGS reads may be used to construct seeds, although with the increasing accuracy of MinION reads, the ends of long reads may be used as well. Apart from cutting out the need for SGS methods, the latter approach has the added advantage that pairs of seeds are connected by at least one long read. Furthermore, as TULIP is not able to assemble regions in which the gap between seeds is larger than the read length, a proper seed density over the entire genome is required. If a marker map is available for the genome, this information can be used to control the distribution of seeds in the selection process.

As a first demonstration of TULIP’s efficiency, Jansen et al. assembled the genome of the European eel Anguilla anguilla (approximately 850Mbp) with 18x coverage in three hours (excluding sequence polishing), requiring only 4.4GB of RAM and four threads²⁰. The resulting assembly was more continuous than the SGS-based reference genome. As was the case with Minimap/Miniasm however, the current quality of MinION reads combined with the lack of an error correction step necessitates post-assembly correction. The authors further showed that missed seed alignments were the most commonly encountered issue during graph simplification, followed by tangled alignments due to repetitive seeds and spurious alignments. The seeds, constructed from short SGS reads, only underwent selection by uniqueness, which did not lead to an equal distribution over the entire genome; however, density remained high enough for successful assembly. The authors noted that assembly using the tips of MinION reads as seeds proved successful for Escherichia coli genomes, but this has not been attempted for larger genomes yet (personal communication, May 1, 2017).

HINGE Although long reads provide a definite edge when attempting to resolve repeat regions, issues may still occur if not all individual repeats are spanned by at least one whole read. In such cases, HINGE may provide a solution. Rather than attempting to resolve frayed rope structures in the assembly graph afterwards, HINGE preprocesses the reads to separate repeat regions that are entirely spanned by a read (and are thus more easily resolvable) from those that are not, and collapses the latter beforehand³⁰.

First, HINGE attempts to identify reads that wholly or partly overlap a repeat region. It does so by performing all-vs-all alignment and then selecting those reads of which a stretch aligns to a proportionally larger number of other reads than the rest of the read. The intuition behind this is that reads from all copies of a repeat region existing in the genome align to each other, thus causing a characteristic abrupt increase in alignments for reads that overlap these repeat regions. Repeat regions covered entirely by at least one read can be easily resolved and are omitted from the following procedure. Of the reads lining the same repeat region, the reads that extend furthest into the repeat region (regardless of the location of the actual copy), are designated “hinges”. In the subsequent greedy extension of the hinges, the path will split at the hinge regions. Like Miniasm, HINGE outputs its assembly in the form of a graph. As its authors show, this is particularly useful for circular genomes.

HINGE provides an elegant solution to long repeat resolution, by separating resolvable regions from unresolvable ones on forehand. Its authors compared HINGE to Miniasm on PacBio reads of 997 circular bacterial genomes and found that overall, HINGE produced a completed genome in more cases than Miniasm could³⁰. Whether the precaution taken by HINGE is necessary is dependent on the genome under consideration and the used reads; if the genome is known to contain repeats longer than most of the reads, the described approach would be justified.

ABruijn While more traditional DBG assemblers performed relatively less well than OLC assemblers on assembling long error-prone reads⁴², the adapted approach taken by the ABruijn assembler has shown more promise³¹. To account for the high error rate, ABruijn filters all k-mers occurring in the reads by their frequency; if a k-mer occurs relatively few times, it is assumed that it contains basecalling errors and it is removed. Then, k-mers are fused into so-called “solid strings”, sequences that contain no other occurring sequences as substring. The ABruijn graph is then drawn by representing solid strings as vertices and connecting them where connections exist in the reads. The edges are weighted by the number of positions between the first bases of the connected solid strings. The assembler consults the weights in this graph to quickly identify overlaps between reads, allowing to select on a minimum overlap length and maximum overhang length. The assembly graph is constructed by starting with the graph for an arbitrary read and continuously extending it by overlapping it with other reads. ABruijn also includes an error correction routine, during which a best consensus between reads is found by identifying low-error stretches and, in between those stretches, choosing the consensus sequence that maximizes the likelihood of the read sequences.

Unfortunately, the authors of ABruijn evaluated their assembler only on E. coli MinION reads of an older and more error-prone chemistry (R7.3), and compared their assembly to that of a Nanocorrect and Celera pipeline, which is currently considered obsolete. In this comparison ABruijn performs better, producing an error rate of 1.1% versus 1.5% for Nanocorrect/Celera, with the majority of errors being deletions. How ABruijn stacks up against currently popular assemblers optimized for MinION data on reads produced with up-to-date chemistries is unclear. The more extensive assessment on PacBio reads provided by the authors shows a dramatic decrease in error rate in the ABruijn assembly in comparison to Canu’s, but requiring comparable running time.

2.4 Post-assembly correction tools

A number of tools attempt to improve assemblies by remapping reads to the assembly and adapting the assembly to increase local resemblance to the reads. These tools may be essential to use after assemblers that do not include a consensus step themselves, such as Minimap/Miniasm and SMARTdeNovo, but have also frequently been used to polish assemblies produced by the other assemblers.

Nanopolish Nanopolish attempts to find an optimal consensus between an assembly and the raw current signal output by the MinION, by iteratively proposing and evaluating small adaptations to the assembly based on the original reads³⁵. The proposal mechanism for adaptations works in two steps. First, reads are aligned to the assembly and the resulting multiple alignment is divided in 50 bp subsequences of the assembly. For each read aligning partly or fully to a subsequence, sections in which 5-mers perfectly align to the assembly are detected. The consensus sequence between each pair of aligning sections is replaced by the aligned read subsequence, creating an initial set of alternative candidate sequences. In the second step, this set is further extended by proposing every possible one-base deletion, insertion and substitution in the previously generated candidate sequences. Of this set, the sequence maximizing the likelihood of observing the raw signal is picked. This process allows Nanopolish to explore a decent number of likely modifications, while remaining computationally tractable.

Polishing an assembly with Nanopolish was found to improve assembly quality, regardless of the assembly tool used. One study on E. coli sequencing data reported that identity to the reference genome rose from 89% to 99% when Nanopolish was applied after Minimap/Miniasm, while improvement after Canu was more modest (98.2% to 99.6%)⁴⁸. Although it addresses all types of errors, a large part of the increase in accuracy is reached due to the correction of homopolymeric stretches. Despite its efficient searching heuristic of block replacement and mutation, running Nanopolish remains a time-consuming step; on the error-prone assembly produced by the otherwise quick Minimap/Miniasm pipeline, it may take up to a few extra days of CPU time to finish polishing^25,48.

Racon Racon corrects MinION assemblies by finding a consensus sequence between reads and the assembly through the construction of partial order alignment (POA) graphs. After alignment of the reads by a mapper of choice (e.g. Minimap or Graphmap), Racon segments the sequence and finds the best alignment between a POA graph of the reads and the assembly. By default, the alignment is performed using the Needleman-Wunsch algorithm, which can align sequence and POA graph with little adaptation (Figure 8). The alignment process is sped up by parallelization. Racon was reported by its authors to be two orders of magnitude faster than the popular (yet currently deprecated) Nanocorrect after assembly of an E. coli genome by Miniasm, albeit not quite as good at diminishing the error rate (to 1.31% versus 0.62% for Nanocorrect). Compared to consensus steps in Falcon⁴⁹ and Canu²⁹ on that same assembly, Racon remains an order of magnitude faster while producing similar error rates. A closer look at the error rate reveals that the majority of errors are indels, validating the need for homopolymer and tandem repeat corrections even in a pipeline containing a polishing step with Racon. Finally, the total genome size estimate following application of Racon was closer to the reference genome size than the estimates of Canu, Falcon and Nanocorrect.

Figure 8. Needleman-Wunsch alignment of two sequences (left), compared to alginment of a partial order alignment (POA) graph and a sequence.

In the latter case, the graph dictates which steps in the matrix are allowed. Taken from Jonathan Dursi’s blog on POA multiple alignment, with permission.