Single molecular sequencing (SMS) developed by Pacific Biosciences and Oxford Nanopore Technologies has been increasingly applied in DNA sequencing studies since its emergence (Ono et al., 2012; Rhoads and Au, 2015). Compared to next-generation sequencing (NGS), SMS can produce considerably longer reads with less sequencing bias and lower cost (Yang et al., 2017; Marchet et al., 2019). Recently, the N50 and maximum lengths of the reads generated by the MinION nanopore sequencer can achieve more than 100 kbp and 1 Mbp, respectively (Chen et al., 2021). Such long read lengths offer new solutions to bioinformatic questions that are hard to be resolved by NGS, such as the de novo genome assembly, genome resequencing, structural variant (SV) discovery, and transcriptome analysis (Stöcker et al., 2016; Wei and Zhang, 2018; Liu et al., 2019; Cao et al., 2021).

Read mapping is to find the possible genomic origins of reads by aligning them against a reference genome. It has become one of the most basic and computing-intensive step in downstream pipelines for SMS dataset analysis (Zhang et al., 2018). However, since the SMS reads have a higher error rate (∼15%) than the NGS reads, most of the read mapping methods developed for NGS short read data, such as Bowtie (Langmead et al., 2009), SOAP3 (Liu et al., 2012), CloudBurst (Michael, 2009), and FEM (Zhang et al., 2018), are not suitable to handle with SMS reads. Thus, a growing number of long noisy read mapping approaches specially designed for SMS long reads have been proposed during the past decade.

Broadly speaking, most existing mapping methods or tools for SMS reads adopt the typical seed-chain-align strategy (Liu et al., 2016a). In brief, they first find the matched k-mers (also called seeds) of query reads to the reference genome; then, the candidate regions (also called alignment skeleton) for aligning in the query read and the reference genome are chosen based on the seeds. Finally, the local subsequences surrounding the seeds are base-to-base aligned to compose the final read alignment with the seeds. Based on the index technique of seed searching, SMS mappers can be categorized into three groups: Burrows–Wheeler Transform full-text minute-space (BWT-FM) index-based (Langmead and Salzberg, 2012), hash table index-based, and short read aligner-based methods.

BWT-FM index-based methods find the matched seeds by constructing the BWT-FM index of the reference genome; such methods include BLASR (Chaisson and Tesler, 2012), BWA-MEM (Li, 2013), lordFAST (Haghshenas et al., 2018), and smsMap (Wei et al., 2020). BLASR (Chaisson and Tesler, 2012) initially builds the BWT-FM index of the genome to find short exact matches; then, the candidate aligned region is generated by using sparse dynamic programming (SDP), and the final detailed alignment within the area defined by SDP is performed by dynamic programming. BWA-MEM (Li, 2013) first searches the supermaximal exact matches according to the BWT-FM index of the reference and detects a group of seeds that are colinear and close to each other by a chaining algorithm, which then extends the seed with a banded affine-gap-penalty dynamic programming. lordFAST (Haghshenas et al., 2018) selects the candidate alignment regions in the genome using the longest matches identified by the BWT-FM index; then, the base-to-base alignment between consecutive seeds is obtained by performing dynamic programming. Recently, the smsMap (Wei et al., 2020) mapper proposed by us also utilizes the BWT-FM index to quickly find the matches in the reference genome and then locates the best starting positions in genome and query read by defining a credibility function. Finally, the detailed alignment result is obtained with the column reduction banded alignment method.

Hash-based mapping methods search the matched seeds through a hash table index of the genome. YAHA (Faust and Hall, 2012) utilizes a hash table index to store the k-mer locations in the reference and then extends seeds to generate fragments that contain contiguous matching bases between the query sequence and the reference genome. Next, it combines the fragments to select potential regions for alignment; lastly, it completes the full alignment by applying a modified version of Smith–Waterman algorithm to the unmatched regions. rHAT (Liu et al., 2015) starts to build the regional hash table (RHT) of the genome; then, the matches of the k-mers within the query read are retrieved through RHT and a direct acyclic graph is built to compose the skeleton of the alignment. At the end, unaligned pairs of segments in the skeleton are aligned with the banded Smith–Waterman algorithm. GraphMap (Ivan et al., 2016) implements the q-gram seeding strategy that allows for fast lookup of inexact matches by constructing a hash index of the reference sequence, which then generates alignment anchors through a fast graph-based ordering of seeds and extends anchors to achieve final alignments. minimap2 (Li, 2018) indexes minimizers in a hash table that allows for fast lookup of exact matches and then identifies colinear anchor sets; final base-level alignments are obtained by applying dynamic programming to regions between adjacent anchors. conLSH (Chakraborty and Bandyopadhyay, 2020) computes context-based locality-sensitive hashing values of genomes to facilitate seed search and then generates a series of sites for candidate alignment after extension; finally, it produces the best possible alignment results by applying the sparse dynamic programming-based approach. Later, S-conLSH (Chakraborty et al., 2021), an improved version of conLSH, was developed by introducing the spaced context-based locality-sensitive hashing for mapping long noisy SMS reads.

Short read aligner-based methods retrieve matches by using extant mappers of short reads. For example, LAMSA (Liu et al., 2016b) first looks for long approximate matches through short read aligner GEM (Marco-Sola et al., 2012), finds a set of possible alignment skeletons, and fills the gaps within the skeletons to generate valid alignments for the whole read. NGMLR (Sedlazeck et al., 2018) identifies similar segments in read and genome by short aligner of NextGenMap (Sedlazeck et al., 2013); then, it extracts the sub-sequences in read and reference to compute a pairwise sequence alignment using a convex gap cost model. Lastly, it reports the set of linear alignments with the highest joint score as the mapping results.

With the rapid development of SMS sequencing technologies, read length is continuously increasing; these longer reads could more frequently span the breakpoints of SVs than those of shorter reads (Kolmogorov et al., 2019; Ren et al., 2021). This may greatly influence read alignment since most state-of-the-art mappers do not consider the SV events, or few methods were designed for handling relatively small variants. Meanwhile, when the matched seeds are diversely distributed in some read parts caused by high sequencing errors, most extant aligners are incapable of obtaining the alignment of the read, leading to many unaligned reads or reads that are partially aligned. As a result, these methods usually focus on producing local mapping results for the query read rather than obtaining the whole end-to-end alignment, leading to a low mapping sensitivity.

To address the aforementioned issues, in this study, we propose a k-mer neighborhood graph mapper (named kngMap), a novel long-read mapping algorithm which is specifically designed to improve mapping sensitivity and deal with SV events. Overall, kngMap works in four main stages. It initially constructs a searching index for the reference genome to quickly find matched k-mers for query reads. Such matches are then used to construct a k-mer d-neighborhood graph where matched k-mers are viewed as vertices and each pair of matched k-mers is connected by a direct edge. Third, a high quality of alignment skeleton is identified by designing a chaining approach. At the end, each unaligned gap in the alignment skeleton is classified into several categories of SV events and filled with a specific alignment method according to its category. The whole read alignment is accomplished by integrating the skeletons and the alignments of the gaps. We benchmarked kngMap on simulated dataset reads with different types of SVs and real-life datasets generated from PacBio SMRT and Oxford Nanopore platforms. The experimental results demonstrated that kngMap has superior ability in terms of base-level sensitivity and end-to-end alignment, which can produce consecutive alignments for the whole read; meanwhile, it can deal with different categories of SVs in the reads.

The main motivation of kngMap is to efficiently improve mapping sensitivity and simultaneously have a superior ability of dealing with SVs for long reads generated by SMS sequencing technologies. The underlying design principle is to effectively find one high-quality alignment skeleton in the reference genome for each read, even though the read contains SVs and more sequencing errors, before the costly procedure of base-to-base alignment to the reference genome. An overview of the kngMap algorithm is depicted in Figure 1. kngMap mainly includes four stages: 1) building a searching index of the reference genome in advance, which is used to quickly find matched k-mers for a query read Figure 1A; 2) constructing a k-mer d-neighborhood graph where matched k-mers are viewed as vertices and each pair of matched k-mers is connected by a direct edge based on the positions in genome and the query read Figure 1B; 3) building and refining the high quality of the alignment skeleton by designing a chaining approach Figure 1C; and 4) classifying the unaligned gaps between pairs of consecutive matched k-mers in the alignment skeleton into several categories of SV events and handling each of them with a specific alignment method according to its category Figure 1D. A more detailed illustration of each step is provided later.

The performance of kngMap was evaluated against 10 state-of-the-art long read aligners: BLASR (v) (Chaisson and Tesler, 2012), minimap2 (v2.12-r829) (Li, 2018), BWA-MEM (v0.7.17-r1194) (Li, 2013), GraphMap (v0.5.2) (Ivan et al., 2016), rHAT (v0.1.2) (Liu et al., 2015), NGMLR (v0.2.7) (Sedlazeck et al., 2018), lordFAST (v0.0.9) (Haghshenas et al., 2018), smsMap (Wei et al., 2020), conLSH (v0.0.1) (Chakraborty and Bandyopadhyay, 2020), and S-conLSH (v0.0.1) (Chakraborty et al., 2021). All methods were benchmarked on simulated and real-life SMS sequencing datasets. All the experiments were performed on a Linux machine server running Ubuntu 14.04 system equipped with two twelve-core (two threads per core) Intel Xeon Gold 5118 CPU @ 2.30GHz and 256 gigabytes (GB) RAM. The run command lines and parameters of each mapping tools are available in Supplementary Table S1.

It is to be noted that with the rapid development of SMS sequencing technology, more than 99% of the raw sequencing data produced by the SMS sequencing platform has a read length of 1,000 bp or longer (Haghshenas et al., 2018). Herein, reads with lengths shorter than 1,000 bp were filtered so that the remaining reads that have 1,000 bp or longer are aligned in the following experiments.

For SMS read mapping, the seed-chain-align procedure (Lin and Hsu, 2020), as is used by most mapping methods, is a classical and effective strategy to obtain the superior alignment results since it can help to seek the alignment skeleton consisting of matched segment pairs between the long SMS read and the reference genome. Generally, the seed-chain-align procedure can be summarized as three phases: first, matched k-mers are found as seeds in one read and in the reference sequence. Then, a group of seeds that are colinear or close to each other is chained as the candidate alignment skeleton. Finally, the non-seed fragments within the alignment skeleton are being extended to generate the base-level alignment. In the stage of generating the alignment skeleton, these methods are usually difficult to capture the matches with long distance, especially for the read parts with SVs and high sequencing errors since that there are few or no matched k-mers in the parts. Thus, most methods failed to generate the alignment skeleton of the query read or just build an alignment skeleton that covers a small part of the read. As a result, these methods cannot obtain the alignment result or just report a local mapping result for the query read, other than obtaining the whole end-to-end alignment, leading to low mapping sensitivity and aligned coverage.

To address the aforementioned challenge, here, we developed kngMap to increase the mapping sensitivity and with 100% aligned coverage for long noisy alignment. kngMap is also a seed-chain-align method using the hashing index technique, compared with other methods. There are three main key features of kngMap: 1) kngMap proposes a scoring strategy in the chaining procedure to choose a group of anchors to form the initial alignment skeleton for each read, even in the situation that the read contains SVs or some regions that matched seeds are dispersedly distributed caused by high sequencing errors; 2) defining an increased credibility function to refine the alignment skeleton, and then a high quality of alignment skeleton is subsequently obtained; and 3) for each of the gaps within the skeleton, kngMap classifies it into one of three categories and implements a specific alignment strategy to fill the corresponding gaps. The scoring strategy and increased credibility function ensures that kngMap can effectively find the alignment skeleton for every query read, even in the situation that the matched seeds are dispersedly distributed in the reference genome. Thus, kngMap can obtain higher mapping sensitivity, that is, align more reads and bases. The last feature can guarantee that the whole end-to-end alignment is obtained, not local alignment achieved by other methods. Therefore, the aligned coverage of kngMap is higher than that of other methods. In addition, in order to further evaluate the performance of the kngMap for a higher error rate up to ∼15%, we applied the PBSIM2 simulator to generate another two simulated datasets with 20% (parameter: --accuracy-mean 20) and 25% (parameter: --accuracy-mean 25) error rates; the other parameter settings are similar to those given in Supplementary Table S2. Supplementary Table S14 shows the evaluation result of kngMap, rHAT, smsMap, lordFAST, BLASR, BWA-MEM, GraphMap, minimap2, and NGMLR on the simulated datasets with 20% error rate. As can be seen from Supplementary Table S14, kngMap correctly mapped the maximum number of reads and bases to the reference genome. Importantly, kngMap achieved 100% aligned coverage, demonstrating that all mapped reads by kngMap are completely aligned. For the base sensitivity and precision in Supplementary Table S14, we can see that kngMap still performed the best, indicating that most bases are correctly mapped in the alignments generated by kngMap. Similar mapping results can be found in Supplementary Table S15 with simulated dataset with 25% error rate. Therefore, these results in Supplementary Tables S14 and S15 demonstrate that our kngMap is more robust to sequencing errors, and it can obtain better mapping results for a higher error rate up to ∼15%.

The massive amount of long noisy read data produced by SMS technologies brings some challenges to existing mapping approaches. In addition to accuracy, computational complexity is another important issue that needs to be considered. The time complexity of kngMap has four main components: 1) In the phase of constructing the reference genome index, it needs to extract all k-mers of the reference genome sequences to build the hash table index; thus, the maximum complexity is in the order of O(G), where G is the length of the genome sequences. 2) In the phase of generating k-mer d-neighborhood graph, it needs to retrieve all the matches of the k-mers for each query read, so the maximum complexity is O(N*L), where N is the number of reads, and L is the average read length. 3) In the phase of generating the skeleton of alignment, it needs to recurrently calculate the score of each node, so the maximum complexity is O(N*M), where M is the average number of matched k-mers for all reads. 4) In the phase of filling the gaps between anchors, it needs to obtain the detailed base-to-base alignment, so the maximum complexity is O(K*Q ² ), where K is the average number of gaps, Q is the average length of gaps, and Q<<L. In summary, the total complexity for kngMap is O(G+N*L+N*M+K*Q ² ). Since Q<<L, kngMap has a time complexity of the order of G and L. In order to graphically evaluate the computational efficiency of our kngMap, we compared kngMap with other mapping tools on the real H. sapiens datasets. Figure 5 shows the relative running time (wall-time) by using the nine tools. In terms of computational efficiency, we can see that the speed of kngMap is almost as fast as minimap2 and is several folds faster than other mapping methods. For example, kngMap is overall about 2-folds faster than rHAT and 6- to 9-folds faster than smsMap, lordFAST, BLASR, BWA-MEM, and GraphMap mappers. The speed comparison result described in Figure 5 indicates that kngMap is efficient to align SMS reads.

Since the infancy of SMS technologies (e.g., PacBio and Oxford Nanopore MinION) that produce longer but higher sequencing error reads, mapping these reads to a reference genome is often the most basic and computing-expensive step for downstream genome sequence analysis. Developing novel long read mapping tools is on demand for improving the mapping sensitivity and effectiveness of SMS read alignment.

In this study, we present kngMap, a new, fast, and highly sensitive mapping algorithm for long noisy reads. Mainly, kngMap contains three key characteristics: 1) kngMap proposes a scoring strategy in the chaining procedure to choose a group of anchors to form the initial alignment skeleton for each read. 2) kngMap designs an increased credibility function to refine the alignment skeleton. The scoring strategy and increased credibility function progressively ensure that kngMap can locate the aligned region for every query read, even in the situation that the matched seeds are dispersedly distributed in the reference genome. 3) For each of the gaps within the skeleton, kngMap classifies it into one of three categories and implements a specific alignment strategy to fill the corresponding gaps, which can help to robustly handle potentially different types of SVs in the reads. kngMap was benchmarked on simulated and real datasets across various genomes with other state-of-the-art mappers. The experimental results demonstrated that kngMap has higher accuracy and sensitivity that can correctly map more sequences and bases to the reference genome and achieves 100% aligned read coverage ratio; meanwhile, it also has good ability to span different types of SVs within the reads.