A survey on de novo assembly methods for single-molecular sequencing

Abstract

Background: The single-molecular sequencing (SMS) is under rapid development and generating increasingly long and accurate sequences. De novo assembly of genomes from SMS sequences is a critical step for many genomic studies. To scale well with the developing trends of SMS, many de novo assemblers for SMS have been released. These assembly workflows can be categorized into two different kinds: the correction-and-assembly strategy and the assembly-and-correction strategy, both of which are gaining more and more attentions.

Results: In this article we make a discussion on the characteristics of errors in SMS sequences. We then review the currently widely applied de novo assemblers for SMS sequences. We also describe computational methods relevant to de novo assembly, including the alignment methods and the error correction methods. Benchmarks are provided to analyze their performance on different datasets and to provide use guides on applying the computation methods.

Conclusion: We make a detailed review on the latest development of de novo assembly and some relevant algorithms for SMS, including their rationales, solutions and results. Besides, we provide use guides on the algorithms based on their benchmark results. Finally we conclude the review by giving some developing trends of third generation sequencing (TGS).

Graphical abstract

Keywords

third generation sequencing / single-molecular real-time sequencing / sequence alignment / sequence error correction / de novo assembly

INTRODUCTION

Third generation sequencing (TGS), also called single-molecule, real-time (SMRT) sequencing, including PacBio RS II platform developed by Pacific Biosciences (PacBio) and MinION platform developed by Oxford Nanopore, captures sequence information directory in the process of DNA molecule replication. TGS generates average read length in the order of 10 kb, which is much longer than next generation sequencing (NGS) read size (100–500 bp) and significantly improves the integrity and continuity of de novo assembly of genomes. Furthermore, the weak effect of classical sequencing bias, such as GC content, allows for resolving complex genome regions and the exploration of structural variants with an unprecedented accuracy and resolution. Many research results based on TGS have been published in different fields, such as de novo assembly, haplotype phasing, N6-methyladenine DNA modification, transcriptome and structural variation detection. TGS has promoted the studies of genomics in an impressive way. However, TGS technology has a serious drawback: its sequencing error rate (11%–15% for PacBio raw reads, 5%–35% for Nanopore raw reads) is much higher than NGS (1%). High error rate and long read length bring great challenges to de novo assembly of TGS sequences. Due to the dramatic difference of data characteristics between TGS sequences and NGS sequences, computation methods designed for NGS sequences always do not work for TGS sequences. Furthermore, sequencing errors between PacBio sequences and Nanopore sequences are also remarkably different, which requires us to address them separately. Currently the widely used de novo assembly workflows developed for TGS sequences can be categorized into two strategies: “correction-and-assembly” and “assembly-and-correction”. The “correction-and-assembly” strategy corrects the sequencing errors in TGS raw reads first, and then assembles genome using corrected reads. On the other hand, the “assembly-and-correction” strategy assembles genome using error-prone TGS raw reads, and then corrects the genome assembly with raw reads. Due to the computational intensive error correction process, the “correction-and-assembly” approach is usually slower. However, directly assembling genome using raw reads with high sequencing error rates may increase assembly errors in genome sequence, especially in genome containing repetitive regions, which then affects the quality of reference genome and causes results bias in downstream analysis. Practices have demonstrated that the “correction-and-assembly” strategy can reconstruct highly continuous and accurate genome assemblies.

In this review we first take a look at the sequencing errors of PacBio sequences and Oxford Nanopore sequences respectively, and then introduce the currently widely used “correction-and-assembly” and “assembly-and-correction” de novo assembly workflows and several relevant methods, such as alignment methods and error correction methods designed for TGS sequences. We summarize all these computation methods in Table 1. To analyze the performance of the workflows, we give two benchmarks and talk about their assembly results on each dataset. Based on their performance on different datasets, we then list some use guides on applying these computation methods. Finally, we talk about the development prospects of sequencing technologies and conclude this paper.

CHARACTERISTICS OF SMS SEQUENCES

Initially the sequence (also called continuous long read, CLR) generated by the PacBio RS system with the first generation of chemistry (C1 chemistry) produces average lengths only around 1500 bp [¹]. Up to now, the PacBio RS II system with the sixth generation of polymerase and the fourth generation of chemistry (P6-C4 chemistry) has improved average read lengths over 10 kb. The maximum read lengths can exceed 60 kb and the N50 can reach up to 20 kb [²]. Error rates of PacBio raw reads are 11%–15%. The large majority of sequencing errors are insertions and deletions (indels).

Due to the random distribution of errors in CLRs, a sufficient sequencing pass yields a circular consensus sequence (CCS) having significantly reduced error rate [³]. For example, an accuracy>99% can be achieved with a coverage of 15 passes. However, considering that the total length of a CLR is limited by the life time of the polymerase [⁴], there exists a tradeoff between the number of sequencing passes and the CCS read length. Recently, the CCS has been optimized to improve the accuracy and length of single-molecule real-time sequencing [⁵]. The average sequence length has been increased to 13.5 kb and the accuracy is as high as 99.8%. This improvement has been successfully used for variant detection and de novo assembly of human genomes.

The lack of bias and the random nature in sequencing errors [⁶] make it possible and relatively easy to correct PacBio reads with high accuracy. Heuristic gapped alignment methods are able to effectively find overlaps between the template (the raw read to be corrected) and its supporting reads (reads sampled from the same regions of the genome with the template) having difference less than 30%. Based on these overlaps, the template can be corrected using some consensus algorithm. The accuracy of corrected reads can be as high as 97%–99%.

The MinION sequencer developed by Oxford Nanopore Technologies (ONT) generates long reads base on the transit of a DNA molecule through a nanoscopic pore which has been thought of as one of the most promising approaches to detect polymeric molecules since the 90s [⁷]. Over the past years, ONT released several different MinION chemistry versions (R6.0, R7.0, R7.3, R9 and R9.4) and several updates of the base-calling tools that strongly improved the performance of the MinON device. While the earliest R6.0 chemistry version allowed to generate a total yield in the order of tens of Mb and 2D reads with 4 kb of median length and 40% of the average error rate, the more recent R7.0 and R7.3 improved the throughput to hundreds of Mb and 2D read size and accuracy to 6 kb and 80%–85%, respectively [⁷]. In October 2016, new flow cells containing R9.4 1D chemistry were released by ONT, which was used to generate the first nanopore high-coverage (30×) human genomes with about 40 flow cells [⁸]. At the same time, a new protocol was also developed to generate ultra-long reads (5×) whose N50 was longer than 100 kb and the maximum read length can reach up to 882 kb. De novo assembly of this 35× raw reads data yielded a contiguous assembly with NG50 as long as 6.4 Mb [⁸].

Compared to PacBio sequences, errors in Nanopore reads are much more complex [⁹,¹⁰]. First, error rates of Nanopore reads are much higher and broadly distributed among reads and among subsequences of a read, ranging from 5%–35%. Second, there commonly exist subsequences with extreme high error rates (50%) in Nanopore reads [⁹]. This means that not only the error rates of two different reads, but also two different subsequences sampled from the same read, can differ greatly. These characteristics of errors make it much more difficult to align, correct and de novo assemble Nanopore reads. Furthermore, due to the higher error rates, average identity of reads corrected from Nanopore reads are generally lower than those corrected from PacBio raw reads. For example, the corrected reads of the 35× human raw sequences [⁸] mentioned above average only 92% identity to the reference; while reads corrected from PacBio raw reads are typically>99% identical to the reference genome. The decreased accuracy of Nanopore corrected sequences introduces troubles for accurately assembling genomes, especially for complex genomes.

CORRECTION-AND-ASSEMBLY METHODS

The “correction-and-assembly” methods achieve high consensus accurate reads by correcting the noisy raw reads and then assemble genomes with the corrected reads. The assembly process generally consists of four stages. (i) Perform a pairwise mapping on the raw reads to find supporting reads for each template. (ii) Align the supporting reads to the template and correct the template based on these overlaps. (iii) Perform a pairwise mapping on the corrected reads and assemble the genome based on the mapping results. (iv) Polish the genome with various data sources, such as the raw reads and NGS sequences using polishing algorithms such as Nanopolish [¹¹], Pilon [¹²], smrtlink [¹³] and Racon [¹⁴] to further improve the accuracy of the genome. It is noted that these polishing methods are also used for correcting errors in genomes assembled by the “assembly-and-correction” approaches.

On the early stage, intuitive and automated hybrid workflows are developed for constructing finished and high-quality small bacterial genomes (1–10 Mb) [¹⁵–¹⁷]. In these hybrid methods, the accurate short reads (from SMRT circular consensus sequencing reads [⁵] or NGS reads) are used to correct errors in the long SMRT sequencing reads. These methods have been applied successfully to a variety of microbes and eukaryotic organisms. They generally involve at least two different sequencing libraries and at least two types of sequencing runs (and sometimes different sequencing platforms). For more efficient genome assembly, especially for large and complex genomes, a homogeneous workflow that requires only one library and sequencing platform is required to simplify the assembly process and to reduce the sequencing costs. Another shortcoming of the hybrid approach is that the constructed genomes always break in regions with insufficient NGS data coverage (owing to the GC or sequence context biases) [¹⁸].

Currently the monohybrid methods that rely solely on TGS data have become the mainstream assembly approaches. These methods avoid the drawbacks in hybrid approaches and have been successfully applied for constructing large genomes (such as the human genome).

Both the hybrid and monohybrid assembly methods involve sequence alignment and error correction as intermediate steps. The sequence alignment for TGS reads is computationally intensive and still remains the most time-consuming step in genome assemblies. In this section, we will first talk about the existing alignment algorithms developed for TGS reads (see the first category of Table 1). We then proceed to introduce the consensus methods that are widely used in the mainstream assembly workflows (see the second category of Table 1). Finally we describe the assembly workflows that have been extensively use for both PacBio sequences and Nanopore sequences.

Alignment methods for TGS sequences

Strictly speaking, the existing alignment methods designed for TGS sequences are all variants of the standard Smith-Waterman algorithm [¹⁹]. They employ various heuristics based on one or more characteristics of TGS reads to reduce the searching space and decrease memory usage of the original Smith-Waterman algorithm. Many of them follow BLAST’s “seed-and-extension” searching paradigm [²⁰].

BLASR

BLASR (basic local alignment with successive refinement) [²¹] uses a successive refinement approach to map SMRT sequences to the reference genome in which case the divergence between them are dominated by insertion and deletion sequencing errors. Its core idea is to narrow down the search space from the whole genome positions to a relatively small number of candidate genome intervals where the read are likely to map. To further validate which candidate interval is mostly likely to be mapped, a detailed alignment is performed on each candidate interval and the best one is chosen as its searching results.

BLASR takes a read r={{\displaystyle r}}_{\mathrm{1}}{{\displaystyle r}}_{\mathrm{2}}\mathrm{...}{{\displaystyle r}}_R and a genome g={{\displaystyle g}}_{\mathrm{1}}{{\displaystyle g}}_{\mathrm{2}}\mathrm{...}{{\displaystyle g}}_G as inputs. The searching process consists of three phases. Step 1: Detecting candidate intervals by clustering short and exact matches. BLASR builds either a suffix array (SA) or BWT-FM index [²²] on the genome according to the time and space requirements, and queries for exact matches of subsequences (of length at least K) between r and g. The BWT-FM index has been used substantially for aligning NGS short reads [²³,²⁴] and accelerating the traditional BLAST search [²⁵,²⁶]. It is a very flexible data structure that allows for searching for fixed length matches, maximal exact matches, and matches in which mismatch bases are allowed efficiently. The set of matches (also called anchors) is denoted as A. The anchors in A are then clustered with global chaining [²⁷]. To do so, the anchors are sorted by genome and query positions and the clusters of anchors are found in intervals that roughly the length of the read. For each anchor {{\displaystyle a}}_i\in A, a set {{\displaystyle A}}_i:=\{{{\displaystyle a}}_j\in A:\mathrm{0}\le G({{\displaystyle a}}_j)+l({{\displaystyle a}}_j)-G({{\displaystyle a}}_j)\le R\} is created, where G({{\displaystyle a}}_j) and l({{\displaystyle a}}_j) are the genome position and length of the anchor {{\displaystyle a}}_j respectively. In each {{\displaystyle A}}_i, the global chaining algorithm [²⁷] is used for finding a maximal subset of anchors {{\displaystyle C}}_i\subseteq {{\displaystyle A}}_i that are not overlapping and listed in increasing read and genome positions. For each cluster {{\displaystyle C}}_i BLASR assigns it with a frequency weighted score {{\displaystyle \Sigma }}_{{{\displaystyle a}}_{j\in }{{\displaystyle C}}_i}\log \left({\displaystyle \frac{\mathrm{1}}{\mathrm{F}\mathrm{r}\mathrm{e}\mathrm{q}({{\displaystyle a}}_j)}}\right), where Freq({{\displaystyle a}}_j) is the frequency of the frequency of subsequence {{\displaystyle a}}_j in the genome. After scoring only the top MAXCANDIDATES (10 by default) clusters are retained. Step 2: Mapping r to the candidate intervals using a sparse dynamic programming (SDP) [²⁸]. For candidate interval {{\displaystyle C}}_i, let a^f and a^l be the anchors in {{\displaystyle C}}_i that have the minimum and maximum genome positions respectively. Suppose the maximum insertion rate of the sequencing device is \delta, then the minimal starting genome position of the interval is s=G(a^f)-(\mathrm{1}+\delta )R(a^f), and the maximal ending genome position is f=G(a^l)+(\mathrm{1}+\delta )(R-(R(a^l)+l(a^l))), here R(a) is the read position of anchor a. The length of this interval is {{\displaystyle l}}_c=f-s. Similar to step 1, a set of matches are found between r and the candidate interval. The matches used here are fixed short length K^{\mathrm{S}\mathrm{D}\mathrm{P}} and denoted as set A^{\mathrm{S}\mathrm{D}\mathrm{P}}. The sparse dynamic programming finds the largest subset of anchors C^{\mathrm{S}\mathrm{D}\mathrm{P}}\subseteq A^{\mathrm{S}\mathrm{D}\mathrm{P}} that are of increasing read and genome positions. Step 3: The sparse dynamic programming alignment is used as a guide for performing a detailed banded alignment. The SDP in Step 2 only aligns some subsequences in r that from C^{\mathrm{S}\mathrm{D}\mathrm{P}} to the candidate interval. To align every base of r to the interval, BLASR triggers a banded dynamic programming. Due to the large indel rate of SMRT sequences, the size of the band used to contain the entire alignment can be prohibitively large. To reduce the searching space, BLASR uses C^{\mathrm{S}\mathrm{D}\mathrm{P}} as a guide where the band follows the layout of the anchors in C^{\mathrm{S}\mathrm{D}\mathrm{P}}. After this step, the read base positions are assigned to reference positions. In BLASR, the authors also provide detailed theory analysis to support their clustering and to compute the mapping quality values.

DALIGN

DALIGN [²⁹] proposes a sensitive and threaded filter for detecting seed points that are likely to have a significant local alignment passing through them and a linear expected time heuristic local alignment algorithm. It is developed for detecting read-to-read overlaps amongst SMRT noisy long reads. Specifically, given two SMRT long noisy read blocks denoted as A and B respectively, DALIGN aims to find local alignments between reads that are sufficiently long and stringent.

The first step of DALIGN is rapid seed detection with a threaded filter. Unlike BLASR, which employs Suffix Array or BWT-FM index to search for variable-length anchors, DALIGN uses a much simpler series of highly optimized sorts to discover k-mer matches in a given diagonal band between two reads. The filter proceeds as follows:

1. Build a k-kmer list {{\displaystyle \mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}}}_A:={{\displaystyle \{(kmer(A^a,i),a,i)\}}}_{a,i} for block A, where kmer(R,i) is the subsequence R[i-k+\mathrm{1},i].

2. Similarly build a k-mer list {{\displaystyle \mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}}}_B for block B.

3. Sort both {{\displaystyle List}}_A and {{\displaystyle List}}_B in order of hash values of their k-mers.

4. Build a k-mer match list {{\displaystyle \mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}}}_M:=\{(a,b,i,j):kmer(A^a,i)=kmer(B^b,j)\} by a merge sweep of the two sorted k-mer lists.

5. Sort {{\displaystyle \mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}}}_M lexicographically on a, b, and i.

Steps 1, 2 and 4 are easily seen to take linear time and space. The most time-consuming step is the sorting process in Steps 3 and 5, which in theory takes O(L\log L) time, where L is the list size. DALIGN optimizes the encoding of the lists in Steps 3 and 5 by squeezing their elements into 64-bit integers. The problem of sorting lists in Steps 3 and 5 now becomes realizing a threadable, memory coherent sort of an array src[\mathrm{0..}N-\mathrm{1}] of N 64-bit integers. DALIGN implements a radix sort [³⁰] where each 64-bit integer is considered as a vector of P=\mathrm{64}/B, B-bit digits ({{\displaystyle x}}_P,{{\displaystyle x}}_{P-\mathrm{1}},\mathrm{...},{{\displaystyle x}}_{\mathrm{1}}) and B is a free parameter whose typical value is 8. The radix sort then sorts the numbers by stably sorting the array on the first B-bit digit {{\displaystyle x}}_{\mathrm{1}}, then on the second {{\displaystyle x}}_{\mathrm{2}}, and so on to {{\displaystyle x}}_p in \mathrm{P} sorting passes. The sorting process asymptotically takes \mathrm{O}(\mathrm{P}(\mathrm{N}+{\mathrm{2}}^{\mathrm{B}})) time. As B and \mathrm{P} are fixed small constants the algorithm is effectively \mathrm{O}(\mathrm{N}) in time complexity.

The second step of DALIGN is rapid local alignment from a seed \rho =(i,j) reported by the filter in Step 1, where i and j are the positions of the k-mer in two reads. The original O(nd) algorithm [³¹] centers on the idea of computing progressive “waves” of furthest reaching (f.r.) points. Starting from \rho =(i,j) on diagonal k=i-j, the goal is to find the longest possible paths starting at \rho, first with 0-differences, then with 1-differences, 2-differences, and so on. The problem is that the waves become wider and wider as we compute f.r. waves away from \rho in each direction. To narrow down the waves, DALIGN uses several strategies to trim the span of a wave by removing the f.r. points that are extremely unlikely to be in the desired local alignment. The key idea is that a desired local alignment should not over any reasonable segment have an exceedingly low correlation. This is also the key idea behind the X-drop criteria in many heuristic alignment tools such as BLAST. In DALIGN, the first principle for trimming a wave is to remove f.r. points for which the last C columns of the alignment have less than M matches. The second trimming principle is keep only f.r. points which are within L antidiagonals of the maximal anti-diagonal reached by its wave. The computation of successive waves eventually ends because either (a) the boundary of either sequence is reached, or (b) all the f.r. points fail the regional alignment quality criterion in which case it is assumed that the two reads no longer correlate with each other.

Experiments show that DALIGN is much more sensitive and can be 22 to 39 times faster than BLASR on the PacBio human dataset. Furthermore, it takes only 15,600 core hours for overlapping the 54× PacBio human dataset. Compared to the 404,000 core hours using BLASR when running on the Google “Exacycle” platform, this represents a substantial 25× reduction in computation time.

Minimap and minimap2

Minimap is an alignment free mapping approach for TGS sequences. In addition to pairwise overlapping, minimap also supports read-to-genome and genome-to-genome mapping. The most distinctive features in minimap is that it takes minimizer matches as seeds and skip the detailed alignment step so that it can be very fast compared to other mapping tools.

Minimizers [³²] are proposed as a replacement of \rho-mers for reducing storage requirements and not reducing the searching sensitivity at the same time. minimap maps a \rho-mer into a 64-bit integer hash values in three \phi (T)=\mathrm{3}steps. First, choose a hash function \phi such that \phi (A)=\mathrm{0}, \phi (C)=\mathrm{1}, \phi (G)=\mathrm{2} and \phi (T)=\mathrm{3}. Next, for a \rho-mer s={{\displaystyle a}}_{\mathrm{1}}{{\displaystyle a}}_{\mathrm{2}}\mathrm{...}{{\displaystyle a}}_k, define

\phi (s):=\phi ({{\displaystyle a}}_{\mathrm{1}})\times {\mathrm{4}}^{k-\mathrm{1}}+\phi ({{\displaystyle a}}_{\mathrm{2}})\times {\mathrm{4}}^{k-\mathrm{2}}+\cdots +\phi ({{\displaystyle a}}_k).

Finally, the function \phi \text{'}:=b\circ \phi is used as the hash function that transfers each \rho-mer to a 64-bit integer hash value. Here b (see algorithm 2 of [³³]) is an invertible integer hash function such that the complex \rho-mer would have the small hash values while the repetitive \rho-mers such as a long run of A’s have large hash values. A (w,k)-minimizer of a string is the \rho-mer having smallest hash values in a surrounding window of w consecutive \rho-mers. Note that there can be more than one minimizer in a window of w consecutive \rho-mers.

Minimap works in two steps: indexing and mapping. In the indexing stage, it collects minimizers in all target sequences, sorts the minimizers by their hash values using radix sort introduced in DALIGN [²⁹], and builds a hash table for them. Minimap only samples roughly \mathrm{1}/w of the total \rho-mers (note that a minimizer is the \rho-kmer having smallest hash value among w consecutive \rho-mers) so that its storage requirement is much less than that in DALIGN. Minimap samples the most informative \rho-mers so that the sensitivity will not reduce significantly. In the mapping stage, minimap collects all the minimizer hits between the query and all the target sequences. A single-linkage clustering is performed to detect the maximal collinear subset of hits by solving a longest increasing sequencing problem (Algorithm 4 of [³³]). This subset is output as the final mapping result.

Minimap has been optimized and improved into minimap2 [³⁴] by adopting a more accurate chaining method and supporting recovering base-level alignments by performing DP-based global alignments between adjacent anchors in a chain. Minimap2 supports aligning many more types of sequences, such as accurate short reads, SMRT long noisy reads, full-length noisy Direct RNA or cDNA reads, as well as assembly contigs or closely related full chromosomes of hundreds of megabases in length.

To chain the minimizer hits (anchors), minimap2 sorts them by ending reference position and let f(i) be the maximal chaining score up to the i-th anchor in the list, which is calculated with a dynamic programming:

f(i):=\max \{{{\displaystyle \max }}_{i\rangle j\ge \mathrm{1}}\{f(i)+\alpha (j,i)-\beta (j,i),{{\displaystyle w}}_i\}\},

where \alpha (j,i):=\min \{\min \{{{\displaystyle y}}_i-{{\displaystyle y}}_j{{\displaystyle x}}_i-{{\displaystyle x}}_i\},{{\displaystyle w}}_i\} is the number of matching bases between the two anchors and \beta (j,i) is the gap cost that defined by

\beta (j,i):=\{\begin{array}{l}\infty ,if\hspace{0.17em}{{\displaystyle y}}_j>{{\displaystyle y}}_i\hspace{0.17em}\mathrm{o}\mathrm{r}\hspace{0.17em}\max \{{{\displaystyle y}}_i-{{\displaystyle y}}_j,{{\displaystyle x}}_i-{{\displaystyle x}}_j\}>G\\ {{\displaystyle r}}_c(({{\displaystyle y}}_i-{{\displaystyle y}}_j)-({{\displaystyle x}}_i-{{\displaystyle x}}_j)),\hspace{0.17em}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}.

In implementation, a gap of length l>\mathrm{0} costs {{\displaystyle r}}_c:=\mathrm{0.01}\cdot \bar{w}\cdot l+\mathrm{0.5}\cdot {{\displaystyle \log }}_{\mathrm{2}}l, where \bar{w} is the average anchor length. A direct computation of all the f(i)’s obviously takes O(N^{\mathrm{2}}) time. Minimap2 reduces this time to O(N) by introducing a heuristic.

To recover detailed alignments, minimap2 performs DP-based global alignments between adjacent anchors in each chain with a 2-piece affine gap cost [³⁵] {{\displaystyle r}}_a(l):=\min \left\{q+l\cdot e,\tilde{q}+l\cdot \tilde{e}\right\}, which helps to recover longer indels. To further accelerate the global alignment, a SSE implementation based on a difference-based formulation developed in [³⁶] is used here.

Edlib

Edlib [³⁷] is developed for efficiently performing an important category of sequence alignments: the calculation of Levenshtein distance (also called edit distance) which is the minimum number of single-character edits (insertion, deletion or substitution). The calculation of Levenshtein distance is equivalent to the calculation of maximum alignment score using standard dynamic programming under the unit cost scoring system. Edlib [³⁷] in principle is an extended implementation of the original bit-vector dynamic programming algorithm [³⁸]. The original bit-vector algorithm accelerates the dynamic programming alignment algorithm which computes the minimal edit distance of two sequences by packing multiple cells as a bit-vector into one CPU register and thus enables the computation the values of multiple cells at a time. It only supports semi-global alignment where gaps at the start and end of the query sequence are not penalized. Edlib extends the bit-vector algorithm to support both global alignments and semi-global alignments in which gaps at the end of the query are not penalized by the extended banded algorithm which supports all three kinds of alignments.

Another deficiency of the original bit-vector algorithm is that it returns no traceback information (a sequence of operations including substitution, insertion and deletion that performed on one sequence to transform it into the other sequence). The traceback is critical to many sequence analysis such as sequencing error correction and variant detection. Edlib remedies this defect by finding the optimal alignment path for all three supported alignments in linear space with Hirschberg’s divide-and-conquer strategy [³⁹].

The consensus algorithms for TGS sequences

The correction of noisy TGS long reads is critical to many genomic studies with high resolution, such as structural variation detection. The currently widely used monohybrid correction methods work by aligning supporting reads to the template sequence. The overlaps are then used to build a multiple sequence alignment, from which a consensus sequence is resolved with some scoring scheme.

DAGCon

DAGCon is developed as part of the hierarchical genome assembly workflow HGAP [⁶] to correct seeding sequences. It draws the idea from [⁴⁰] in which multiple sequence alignment is represented as a directed acyclic graph.

DAGCon constructs a directed acyclic graph which is denoted as the sequence alignment graph (SAG) with alignments between supporting reads and the template in five steps. Step 1: Normalize mismatches and degenerated alignment positions in each alignment. It is known that mismatches in one alignment are indistinguishable from an adjacent insertion-deletion pair to the correction algorithm. Furthermore, given a specific gap cost function, there are cases where there is more than one way to place gaps in an alignment. To reduce such potential inconsistency introduced by the two cases mentioned above and to improve the final consensus accuracy, each mismatch is replaced by a insertion-deletion pair and every gap is moved to the right most equivalent position (see Supplementary Figure S1 in [⁶] for an example). Step 2: Setup an initial directed acyclic graph from the template. Specifically, an initial backbone graph consists of the nodes {{\displaystyle v}}_i with \mathrm{l}\mathrm{a}\mathrm{b}\mathrm{e}\mathrm{l}({{\displaystyle v}}_i)={{\displaystyle R}}_i, where {{\displaystyle R}}_i is the residue of the template at position i, and the directed edges {{\displaystyle e}}_{i,j+\mathrm{1}}:={{\displaystyle v}}_i\to {{\displaystyle v}}_{i+\mathrm{1}} is built in this step. Step 3: Construct a multigraph from the alignments and merge the edges that are connecting the same nodes. Each alignment is represented as a list A=\{({{\displaystyle q}}_i,{{\displaystyle r}}_i):{{\displaystyle q}}_i,{{\displaystyle r}}_i\in A,C,G,T,-\} and added to the backbone graph (Algorithm 1 of [⁶]). Step 4: Merge nodes. To reduce the complexity of the backbone graph, DAGCon merges those nodes that are similar to each other (for example, the nodes having the same label and are connected to the same out-node or are connected from the same in-node) while preserves the directed acyclic property of the graph at the same time (Algorithm 4 of [⁶]). Step 5: Generate consensus sequence from the simplified backbone graph. Each node in the graph is given a score by checking the weight of its out-edges. Each path in the graph is assigned a score by summing up the scores of its nodes. The consensus is constructed by finding the path having the highest path score among all the possible paths with a simple dynamic programming (Algorithms 5 and 6 of [⁶]).

FALCON-sense

FALCON-sense is the consensus algorithm that released together with the FALCON assembler [⁴¹]. FALCON-sense starts with performing a pairwise mapping with DALIGN [²⁹] on the noisy raw reads. Each template is then processed as follows. First, a variation of the O(nd) algorithm [³¹] is applied to align the supporting reads to the template where all the differences in the alignments are encoded as insertions and deletions only. For each position in every alignment, an “alignment tag” that encodes both the positions and the bases in the supporting read and the template is created. These tags are then used for constructing an alignment graph in the following way. For each tag, if it does not exist in the graph, a corresponding node is added to the graph. For consecutive tags from an alignment, an edge is created with edge\_count equaling to 1 is added to the graph. If the edge already exists in the graph, its corresponding edge\_count is increased by 1. After processing all the alignment tags we obtain a directed acyclic graph. Every node in the graph has a tag. The edge\_count in each edge indicates how many reads support the connection between the two nodes. Finally a standard dynamic programming algorithm is applied on the alignment graph to find the highest confidence path. By concatenating the bases of the tags on the best path the consensus sequence is generated. Supplementary Figure 12 of [⁴¹] illustrates this consensus process.

In implementation, FALCON-sense is much simpler, more robust and more efficient than the DAGCon algorithm. It is also integrated in another widely used de novo assembly workflow Canu as the error correction method.

The correction-and-assembly workflows

There are three widely used “assembly-and-correction” assemblers, which are described below.

FALCON

The FALCON assembler [⁴¹] follows the same workflow of the hierarchical genome assembly pipeline HGAP [¹⁸] with many computationally optimized implementations. The most distinctive feature of FALCON is that it also supports phased diploid genome assembly with the FALCON-Unzip algorithm.

The FALCON workflow starts with error correction using the FALCON-sense algorithm (Section FALCON-sense). After that the DALIGN [²⁹] algorithm is used for identifying overlaps between all pairs of corrected reads. These overlaps are used to construct a directed string graph [⁴²], which is further simplified with assembly string graph reduction (Supplementary Figure 14 of [⁴¹]) into the final “haplotype-fused assembly graph”. After the reduction process, the primary and associated contigs are constructed based on the corrected reads and their overlap relationships. For each primary contig, FALCON collects all the raw reads that associated with it and with its associated contigs according to the overlapping information. The raw reads are then aligned to the contigs with BLASR [²¹] for phasing heterozygous SNPs and identifying the haplotype of each raw read (Supplementary Figure 15 of [⁴¹]).

To construct haplotype-specific contigs (haplotigs), the FALCON-Unzip algorithm constructs a haplotype-specific assembly graph for each contig from its corresponding raw reads and then combines this graph to the fused assembly subgraph that contains the paths of this contig to construct a complete subgraph, which has complete read representation from both haplotypes. See Figure 1 and Supplementary Figure 13 of [⁴¹] for details. As a result, FALCON-Unzip build a new primary contig and several haplotigs from the original assembly graph of that contig. It assigns each phased read to the primary contig or one of the haploptigs according to the phasing information.

Canu

The PBcR assembly workflow [⁴³] is developed for effectively de novo assembling large genomes such as the human genome. It proposes MHAP for overlapping noisy long reads using probabilistic and locality-sensitive hashing. The DAGCon [⁶] and FalconSense are used for correcting the errors in the long reads. And the Celera Assembler is used for assembling the corrected long reads.

To reduce the running time of the overlapping step in the assembly pipeline, MHAP (MinHash Alignment Process) uses a dimensionality reduction technique named MinHash sketches for efficient alignment filtering. The MinHash sketch of a read is created by extracting all the \rho-mers and converting them to fingerprints by H hash functions. As a result, we obtain H hash sets {\mathrm{\Gamma }}_{\mathrm{1}},\cdots ,{\mathrm{\Gamma }}_H. In MHAP, after the first hash set {\mathrm{\Gamma }}_{\mathrm{1}}, subsequent fingerprints are generated using XOR shift pseudo-random number generators {\mathrm{\Gamma }}_{\mathrm{2}},\cdots ,{\mathrm{\Gamma }}_H. The \rho-mer that generates the minimum value for each hash is referred to as the min-mer for that hash. The sketch of that sequence is then composed of the ordered set of its H min-mer finerprints that is much smaller than the set of all of the read’s \rho-mers. To find overlaps between two reads, MHAP builds two sketches for them and estimates their Jaccard similarity by computing the fraction of min-mers common to both sketches. If these two sketches are sufficiently similar, the shared min-mers are located in two reads and the median difference in positions is computed to determine the overlap offsets. See Figure 1 of [⁴³] for an example. The concept of sketch presented here and the concept of minimizer that used in minimap and minimap2 share the same idea. They are both reduced representations of the original sequence that used for decreasing memory usage and accelerating the overlapping process.

The PBcR is later optimized into a more efficient assembler named Canu [⁴⁴]. Canu implements an updated version of MHAP which uses a two-stage overlap filter, where the first stage identifies read pairs that are likely to share an overlap and the second stage estimates the extent and quality of the overlap. For the first step, a tf-idf (term frequency, inverse document frequency) weighting is implemented in MHAP to increase sensitivity of true overlaps and reduce the number of false and repetitive overlaps. For a raw read set Q, let {{\displaystyle f}}_{\max } and {{\displaystyle f}}_{\min } be the maximum and minimum observed frequency of all the \rho-mers in Q respectively, and {{\displaystyle f}}_q be the frequency of a specific \rho-mer q in Q. The inverse document frequency function for q, denoted as {{\displaystyle idf}}_q, is defined as

{{\displaystyle idf}}_q:=T\left(\log \left({\displaystyle \frac{{{\displaystyle f}}_{\max }}{{{\displaystyle f}}_q}}-a\right)\right).

The parameter a\in [\mathrm{0},\mathrm{1}] controls how strongly less common \rho-mers are preferred in relation to the more common ones, and T is a constant used for linearly rescaling the idf value to fall in [\mathrm{1},{{\displaystyle f}}_{\max }]. A \rho-mer does not exist in Q will be given {{\displaystyle idf}}_{\max } as its idf value. The tf-idf value of q is then computed by the formulation {{\displaystyle w}}_q:={{\displaystyle tf}}_q\cdot {{\displaystyle idf}}_q, where {{\displaystyle tf}}_q is the number of occurrences of q in Q. After computing the tf-idf weight {{\displaystyle w}}_q for each k-mer q, the original MinHash computation is now replaced by applying {{\displaystyle w}}_q hash functions {\mathrm{\Gamma }}_{i,\mathrm{1}},\cdots ,{\mathrm{\Gamma }}_{i,{{\displaystyle w}}_q} per entry rather than the single {{\displaystyle \Gamma }}_i. For each sketch entry {{\displaystyle S}}_i, the min-mer is chosen as the minimum hash value computed across all the hash functions. Since highly weighted k-mers are hashed more times, this increases the chance that they will be chosen as a min-mer.

In the correction stage, Canu uses two filtering steps to avoid the introduction of overlaps from copy-specific repeat variants. The first is a global filter where each read chooses where it will supply correction evidence, and the second is a local filter where each read accepts or rejects the evidence supplied by other reads. The FALCON-sense consensus algorithm [⁴¹] is used for generating the corrected reads. Canu also reduces the runtime by an order of magnitude of the Celera Assembler. In the assembly stage, Canu uses a variant of the greedy “best overlap graph” algorithm [⁴⁵] for constructing a sparse overlap graph. The greedy algorithm loads only the “best” (longest) overlaps for each read end into memory and thus consumes much less memory compared the string graph formulation [⁴²], which requires loading the full overlap graph into memory. Furthermore, to improve repeat and haplotype separation, Canu introduces a new “Bogart” algorithm to statistically filter repeat-induced overlaps and retrospectively inspect the graph for potential errors.

MECAT

MECAT [⁴⁶] is developed to effectively distinguish true overlaps and random read pairs that share similar subsequences. Traditional aligners such as BLAST [²⁰] employ the \rho-mer match methods to filer out random read pairs and quickly find seeding positions for local gapped alignments. However, the repetitive subsequences in reads always find a large number of \rho-mer matches and results in excessive candidate seeding positions. On the other hand, simply masking these low-complexity regions, or ignoring the \rho-mers with high frequencies can cause loss of correct overlaps. To address this issue, MECAT proposes a pseudolinear alignment distance difference factors scoring (DDFs) algorithm to quickly filter excessive alignments.

For two k-mer matches whose positions are {{\displaystyle s}}_{\mathrm{1}}=({{\displaystyle q}}_{\mathrm{1}},{{\displaystyle r}}_{\mathrm{1}}) and {{\displaystyle s}}_{\mathrm{2}}=({{\displaystyle q}}_{\mathrm{2}},{{\displaystyle r}}_{\mathrm{2}}) respectively, where {{\displaystyle q}}_i is the read position and {{\displaystyle r}}_i is the reference position, their distance difference factor is defined as

\mathrm{D}\mathrm{D}\mathrm{F}({{\displaystyle s}}_{\mathrm{1}},{{\displaystyle s}}_{\mathrm{2}}):=\left|\mathrm{1}-{\displaystyle \frac{{{\displaystyle q}}_{\mathrm{1}}-{{\displaystyle q}}_{\mathrm{2}}}{{{\displaystyle r}}_{\mathrm{1}}-{{\displaystyle r}}_{\mathrm{2}}}}\right|.

If , then this two \rho-mer matches are said to support each other in positions. Now suppose the \rho-mer matches between a read and a target sequence are {{\displaystyle s}}_{\mathrm{1}}, {{\displaystyle s}}_{\mathrm{2}}, …, {{\displaystyle s}}_n. The DDF score of each match is initialized with zero. For every match pair {{\displaystyle s}}_i and {{\displaystyle s}}_j, if , we increase the DDF scores of both matches by one. The k-mer having highest DDF score provides a position to which a significant local alignment is likely to pass through. It is evidently that the DDF scoring process takes O(n^{\mathrm{2}}) time, which is very slow for long SMRT reads. To speed up the scoring process, MECAT breaks it into two steps so that the quadratic mutual scoring only takes place in local regions. In the first step the database is decomposed into consecutive blocks. Each block has the same length of B whose default value is 2000. The scoring starts with randomly selecting a block that having at least b\rho-mer matches. The mutual DDF scoring is performed on the matches into that block. If the match with the highest score is significant (greater than the threshold), we set it as the seed position. In the second step, the scoring process is extended to its neighbor blocks. For each neighbor block, we calculate the DDF between its \rho-mer matches and the seed k-mer match. If , we increase the score of the seed \rho-mer match by 1. The mutual scoring in step 1 is conducted in O(N^{\mathrm{2}}) time, and the extension scoring in step 2 is conducted in O(N) time, where N is the number of \rho-mer matches. Since the number of \rho-mer matches in mutual scoring is small, the overall scoring process can be finished in pseudolinear time.

Experiments show that the scores of the seeding positions grow linearly with their alignment lengths. After DDF scoring filtering, the candidate seeding positions are reduced by 50%–70% before extending the seeding positions to local alignments. In the correction stage, MECAT employs DDF scoring to effectively choose supporting reads for each template and then uses an adaptive DAGCon consensus algorithm [¹⁸] to correct the templates. Originally the DAGCon algorithm achieves higher precision compared to other correction methods and builds a sequence alignment graph (SAG) for the whole template, which consumes a lot of time and memory. MECAT accelerates the consensus process in two steps. First it identifies bases in the template that are covered mainly by the same bases or gaps from supporting reads. In the corrected read, these bases will remain unchanged if they are covered mainly by the same bases, or will be deleted if they are mainly covered by gaps. Second, the DAGCon algorithm is used for correcting the subsequences between two unchanged bases. MECAT originally uses the assembly module from Canu [⁴⁴], which is now replaced by fsa, a string graph [⁴²] based assembler.

ASSEMBLY-AND-CORRECTION METHODS

Flye

Flye [⁴⁷] is developed for accurately resolving genomic repeats in de novo genome assembly. Unlike other mainstream assemblers, Flye starts with creating a disjointigs that represent concatenations of multiple disjoint genomic segments and concatenating all error-prone disjointtigs into a single string. From these concatenate results, Flye creates an accurate assembly graph, untangle it with reads and resolves the bridged repeats. Finally, Flye resolves the unbridged repeats with the repeat graph using small differences between repeat copies and outputs accurate contigs formed by paths in the graph.

Wtdbg2

Wtdbg2 [⁴⁸] proposes a fast pairwise mapping implementation and a layout algorithm based on fuzzy-Bruijn graph (FBG) to speed up the large genome assembly process. To reduce memory usage, the current aligners (such as DALIGN and minimap2) split the input raw reads into small batches and perform pairwise mapping between batches. Wtdbg2 loads all the raw reads into memory and builds a hash table for the \rho-mers in all raw reads. For two reads that share \rho-mers, wtdbg2 performs a Smith-Waterman-like dynamic programming between these two reads in terms of bins. Wtdbg2 encodes each read into a consecutive sequence of bins. Each bin consists of 256 bp. The FBG proposed by wtdbg2 extends the basic ideas from de Bruijn graph. Each base in FBG is a 256 bp bin and a \rho-bin in FBG consists of \rho consecutive bins from reads. After simplifying the FBG, wtdbg2 builds the final consensus with partial order alignment over edge sequences.

SOME BENCHMARKS

To illustrate the efficiency and assembly quality of the assembly workflows described in this review, we present some benchmarks in this section. The benchmark results would also provide clues for use guides to be talked about in the next section. It is noted that there exits an abundance of benchmarks of this sort reported in various references. Therefore we are not going to carry out such tests ourselves but instead quoting two representative reports and listing the results here for easy reference and discussion.

Benchmarks on PacBio data

We list a simplified version of experimental results reported in [⁴⁸] Table 2. The five workflows are used for de novo assembling three moderate genomes (100 M) and one large human genome (3G). More details about the datasets and the experimental setup are found in Table 1 and Supplementary Table 1 of [⁴⁸]. On the three moderate genomes, Canu achieves the highest assembly quality while has the worst efficient performance. FALCON, Flye and MECAT share roughly the same behavior in terms of NG50, NG75 and run time. Wtdgb2 is the fastest assembler and is at least one order of magnitude faster than Canu. NG50 and NG75 of assemblies constructed by Wtdbg2 are generally better than those of FALCON, Flye and MECAT. In assembling the 100X human dataset, FALCON has the best assembly quality. The assembly process of Wtdbg2 is dominated by the polishing step, which consumes about 75% of the computation time. The reference genome cover of the Wtgbg2 assembly is significantly smaller than that of Canu and FALCON.

Benchmarks on Nanopore data

The FALCON and MECAT workflows do not work for Nanopore data. Therefore we assess performance of the other three assemblers here. In Table 3 we quote the results from [⁴⁹] where the three workflows are used for assembling five moderate genomes (100M). Details about the datasets and experiment setup can be found in Table 1 and Supplementary Table 8 of [⁴⁹]. As on PacBio datasets, Canu takes the longest time. Wtdbg2 is still the most efficient workflow and can be two orders of magnitude faster than Canu.

On the other hand, due to the great discrepancy between PacBio and Nanopore sequences, the three assemblers behave quite differently on these two kinds of datasets when considering the assembly contiguity. The Flye workflow has the longest N50 on each dataset. Wtdbg2 obtains better N50 on A. thaliana and D. melanogaster datasets than Canu. While on complex genomes C. reinhardtii, O. sativa and S. pennellii, Canu produces larger N50 than Wtdbg2. This shows the advantage of the Correction-and-Assembly strategies on assembling complex genomes. We also note that Wtgbd2 outputs significantly more contigs than the other two workflows on each dataset.

SOME USE GUIDES

If your task is to carry out a sequence alignment, the preferred alignment tool is generally minimap2 [³⁴], which can meet your most demands of processing. As a general-purpose sequence mapping tool, minimap2 always runs quite satisfactorily in both efficiency and sensitivity. Minimap2 supports both read-to-read pairwise mapping and reference mapping, and a large spectrum of sequences: DNA and mRNA sequences, PacBio and Nanopore noisy long reads, as well as assembly contigs or closely related full chromosomes of hundreds of megabases in length.

If your study (not necessarily be de novo genome assembly) involves raw reads correction as an intermediate step, you don’t bother re-implementing the consensus algorithms mentioned above and then proceeding to carry out the overlap-and-correction pipeline. To correct PacBio noisy long reads, either Canu or MECAT is recommended. However if the dataset is large, MECAT is the preferred option, since MECAT has astounding speed benefits and maintains the same level of accuracy assurance as Canu. To correct Nanopore raw reads, the currently available mainstream tool is Canu.

There exist multiple choices of workflows for de novo genome assembly with both PacBio and Nanopore sequences. For small and moderate datasets, especially for complex genomes, Canu is able to construct contiguous and accurate assemblies solidly. For large PacBio datasets, if the genome is complex, then MECAT is recommended. Otherwise Flye and Wtdbg2 are also appropriate alternatives. For large Nanopore datasets, consider Flye. However if your computing resources is limited you should use Wtdbg2. Note that there are cases in which these workflows can work in a cooperative way. For example, if the quality of contigs assembled by one assembler is stunningly poor, you may wonder if the raw read dataset is corrupt (that is, the data is incomplete or too noisy). To validate your point, you may construct the contigs again with a second fast assembler such as Wtdbg2 or Flye and investigate the quality of the new contigs.

CONCLUSIONS AND DISCUSSIONS

In this survey, we analyze in detail the features of SMS sequences including both the PacBio platform and the Nanopore platform. We review and summarize several important and widely used computation methods that are designed for TGS. We classify these methods into three categories: alignment methods, consensus methods and assembly workflows. The alignment methods and consensus methods are integrated into assembly workflows to carry out intermediate tasks. For each method, we give a general description to its design principle and an in-depth analysis of its core concepts and core ideas. To illustrate the performance of the assembly workflows, we show two benchmark results here and analysis their behavior on each dataset. Based on their design principles and benchmark results, we proceed to provide some general use guides on applying these methods.

Though has made substantial progress, TGS still has defects and the field is still finite that TGS is used in genomic studies. Based on what we have learned, we predict here that in the future TGS will develop towards the following ways. Firstly, as the ultra-long read protocol [⁸] and the circular consensus sequencing continue to evolve, SMS sequences should keep increasing in both length and accuracy. As results, we get more and more complete and contiguous assembled contigs that cover regions remain unresolved today. It will be also feasible to perform polyploid assembly of complex genomes. Besides de novo assembly, it is expected that TGS will promote studies in other fields such as N6-methyladenine DNA modification, transcriptome and structural variation detection. Secondly, it is also expected that the third generation sequencing technology also develop in different areas such as RNA and single-molecule sequencing and to promote the relevant studies with higher resolutions.

We hope that this review can serve as a roadmap for those whose are interested in TGS and serve as a foundation based on which more advanced computation methods are developed.

References

[1] Brown, S. D., Nagaraju, S., Utturkar, S., De Tissera, S., Segovia, S., Mitchell, W., Land, M. L., Dassanayake, A. and Köpke, M. (2014) Comparison of single-molecule sequencing and hybrid approaches for finishing the genome of Clostridium autoethanogenum and analysis of CRISPR systems in industrial relevant Clostridia. Biotechnol. Biofuels, 7, 40

[2] Rhoads, A. and Au, K. F. (2015) PacBio sequencing and its applications. Genom. Proteom. Bioinf., 13, 278–289

[3] Eid, J., Fehr, A., Gray, J., Luong, K., Lyle, J., Otto, G., Peluso, P., Rank, D., Baybayan, P., Bettman, B., (2009) Real-time DNA sequencing from single polymerase molecules. Science, 323, 133–138

[4] Detter, J. C., Johnson, S. L., Bishop-Lilly, K. A., Chain, P. S., Gibbons, H. S., Minogue, T. D., Sozhamannan, S., Van Gieson, E.J. and Resnick, I. G. (2014) Nucleic acid sequencing for characterizing infectious and/or novel agents in complex samples. In: Biological Identification, pp. 3–53. Sawston: Woodhead Publishing

[5] Wenger, A. M., Peluso, P., Rowell, W. J., Chang, P. C., Hall, R. J., Concepcion, G. T., Ebler, J., Fungtammasan, A., Kolesnikov, A., Olson, N. D., (2019) Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome. Nat. Biotechnol., 37, 1155–1162

[6] Chin, C. S., Alexander, D. H., Marks, P., Klammer, A. A., Drake, J., Heiner, C., Clum, A., Copeland, A., Huddleston, J., Eichler, E. E., (2013) Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. Nat. Methods, 10, 563–569

[7] Magi, A., Semeraro, R., Mingrino, A., Giusti, B. and D’Aurizio, R. (2018) Nanopore sequencing data analysis: state of the art, applications and challenges. Brief. Bioinformatics, 19, 1256–1272

[8] Jain, M., Koren, S., Miga, K. H., Quick, J., Rand, A. C., Sasani, T. A., Tyson, J. R., Beggs, A. D., Dilthey, A. T., Fiddes, I. T., (2018) Nanopore sequencing and assembly of a human genome with ultra-long reads. Nat. Biotechnol., 36, 338–345

[9] Magi, A., Giusti, B. and Tattini, L. (2017) Characterization of MinION nanopore data for resequencing analyses. Brief. Bioinformatics, 18, 940–953

[10] Rang, F. J., Kloosterman, W. P. and de Ridder, J. (2018) From squiggle to basepair: computational approaches for improving nanopore sequencing read accuracy. Genome Biol., 19, 90

[11] Loman, N. J., Quick, J. and Simpson, J. T. (2015) A complete bacterial genome assembled de novo using only nanopore sequencing data. Nat. Methods, 12, 733–735

[12] Walker, B. J., Abeel, T., Shea, T., Priest, M., Abouelliel, A., Sakthikumar, S., Cuomo, C. A., Zeng, Q., Wortman, J., Young, S. K., (2014) Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One, 9, e112963

[13] Kingan, S. B., Heaton, H., Cudini, J., Lambert, C. C., Baybayan, P., Galvin, B. D., Durbin, R., Korlach, J. and Lawniczak, M. K. N. (2019) A high-quality de novo genome assembly from a single mosquito using PacBio sequencing. Genes (Basel), 10, 62

[14] Vaser, R., Sović I., Nagarajan, N. and Šikić M. (2017) Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res., 27, 737–746

[15] Bashir, A., Klammer, A., Robins, W. P., Chin, C. S., Webster, D., Paxinos, E., Hsu, D., Ashby, M., Wang, S., Peluso, P., (2012) A hybrid approach for the automated finishing of bacterial genomes. Nat. Biotechnol., 30, 701–707

[16] Koren, S., Schatz, M. C., Walenz, B. P., Martin, J., Howard, J. T., Ganapathy, G., Wang, Z., Rasko, D. A., McCombie, W. R., Jarvis, E. D., (2012) Hybrid error correction and de novo assembly of single-molecule sequencing reads. Nat. Biotechnol., 30, 693–700

[17] Ribeiro, F. J., Przybylski, D., Yin, S., Sharpe, T., Gnerre, S., Abouelleil, A., Berlin, A. M., Montmayeur,A., Shea, T. P., Walker, B. J., (2012) Finished bacterial genomes from shotgun sequence data. Genome Res., 22, 2270–2277

[18] Chin, C. S., Alexander, D. H., Marks, P., Klammer, A. A., Drake, J., Heiner, C., Clum, A., Copeland, A., Huddleston, J., Eichler, E. E., (2013) Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. Nat. Methods, 10, 563–569

[19] Smith, T. F. and Waterman, M. S. (1981) Identification of common molecular subsequences. J. Mol. Biol., 147, 195–197

[20] Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., Miller, W. and Lipman, D. J. (1997) Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res., 25, 3389–3402

[21] Chaisson, M. J. and Tesler, G. (2012) Mapping single molecule sequencing reads using basic local alignment with successive refinement (BLASR): application and theory. BMC Bioinformatics, 13, 238

[22] Ferragina, P. and Manzini, G. (2000) Opportunistic data structures with applications In: Proceedings the 41st Annual Symposium on Foundations of Computer Science, pp. 390–398. IEEE

[23] Li, H. and Durbin, R. (2010) Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics, 26, 589–595

[24] Ben, L. and Salzberg, S. L. (2012) Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359

[25] Chen, Y., Ye, W., Zhang, Y. and Xu, Y. (2015) High speed BLASTN: an accelerated MegaBLAST search tool. Nucleic Acids Res., 43, 7762–7768

[26] Lam, T. W., Sung, W. K., Tam, S. L., Wong, C. K. and Yiu, S. M. (2008) Compressed indexing and local alignment of DNA. Bioinformatics, 24, 791–797

[27] Abouelhoda, M. I. and Ohlebusch, E. A. (2003) Local chaining algorithm and its applications in comparative genomics. In: International Workshop on Algorithms in Bioinformatics, pp. 1–16. Berlin: Springer

[28] Eppstein, D., Galil, Z., Giancarlo, R. and Italiano, G. F. (1992) Sparse dynamic programming I: linear cost functions. J. Assoc. Comput. Mach., 39, 519–545

[29] Myers, G. (2014) Efficient local alignment discovery amongst noisy long reads In: International Workshop on Algorithms in Bioinformatics, pp. 52–67. Berlin: Springer

[30] Cormen, T. H., Leiserson,C. E., Rivest, R. L. and Stein, C. (2009) Introduction to Algorithms (3rd. edition), pp.197–204. Cambridge: MIT Press

[31] Myers, E. W. (1986) AnO (ND) difference algorithm and its variations. Algorithmica, 1, 251–266

[32] Roberts, M., Hayes, W., Hunt, B. R., Mount, S. M. and Yorke, J. A. (2004) Reducing storage requirements for biological sequence comparison. Bioinformatics, 20, 3363–3369

[33] Li, H. (2016) Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences. Bioinformatics, 32, 2103–2110

[34] Li, H. (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics, 34, 3094–3100

[35] Gotoh, O. (1990) Optimal sequence alignment allowing for long gaps. Bull. Math. Biol., 52, 359–373

[36] Suzuki, H. and Kasahara, M. (2018) Introducing difference recurrence relations for faster semi-global alignment of long sequences. BMC Bioinformatics, 19, 45

[37] Šošic, M. and Šikic, M. (2017) Edlib: a C/C ++ library for fast, exact sequence alignment using edit distance. Bioinformatics, 33, 1394–1395

[38] Myers, G. (1999) A fast bit-vector algorithm for approximate string matching based on dynamic programming. J. Assoc. Comput. Mach., 46, 395–415

[39] Hirschberg, D. S. (1975) A linear space algorithm for computing maximal common subsequences. Commun. ACM, 18, 341–343

[40] Lee, C., Grasso, C. and Sharlow, M. F. (2002) Multiple sequence alignment using partial order graphs. Bioinformatics, 18, 452–464

[41] Chin, C. S., Peluso, P., Sedlazeck, F. J., Nattestad, M., Concepcion,G. T., Clum, A., Dunn, C., O’Malley, R., Figueroa-Balderas, R., Morales-Cruz, A., (2016) Phased diploid genome assembly with single-molecule real-time sequencing. Nat. Methods, 13, 1050–1054

[42] Myers, E. W. (2005) The fragment assembly string graph. Bioinformatics, 21(suppl_2), ii79–ii85

[43] Berlin, K., Koren, S., Chin, C. S., Drake, J. P., Landolin, J. M. and Phillippy, A. M. (2015) Assembling large genomes with single-molecule sequencing and locality-sensitive hashing. Nat. Biotechnol., 33, 623–630

[44] Koren, S., Walenz, B. P., Berlin, K., Miller, J. R., Bergman, N. H. and Phillippy, A. M. (2017) Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res., 27, 722–736

[45] Miller, J. R., Delcher, A. L., Koren, S., Venter, E., Walenz, B. P., Brownley, A., Johnson, J., Li, K., Mobarry, C. and Sutton, G. (2008) Aggressive assembly of pyrosequencing reads with mates. Bioinformatics, 24, 2818–2824

[46] Xiao, C. L., Chen, Y., Xie, S. Q., Chen,K.-N., Wang, Y., Han, Y., Luo F., Xie , Z. (2017) MECAT: fast mapping, error correction, and de novo assembly for single-molecule sequencing reads. Nat. Methods, 14, 1072–1074

[47] Kolmogorov, M., Yuan, J., Lin, Y. and Pevzner, P. A. (2019) Assembly of long, error-prone reads using repeat graphs. Nat. Biotechnol., 37, 540–546

[48] Ruan, J. and Li, H. (2020) Fast and accurate long-read assembly with wtdbg2. Nat. Methods, 17, 155–158

[49] Chen, Y., Nie, F., Xie, S.-Q. and Zheng, Y.-F. (2020) Fast and accurate assembly of Nanopore reads via progressive error correction and adaptive read selection. bioRxiv, 930107

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature