S-conLSH: alignment-free gapped mapping of noisy long reads

Abstract

Background

The advancement of SMRT technology has unfolded new opportunities of genome analysis with its longer read length and low GC bias. Alignment of the reads to their appropriate positions in the respective reference genome is the first but costliest step of any analysis pipeline based on SMRT sequencing. However, the state-of-the-art aligners often fail to identify distant homologies due to lack of conserved regions, caused by frequent genetic duplication and recombination. Therefore, we developed a novel alignment-free method of sequence mapping that is fast and accurate.

Results

We present a new mapper called S-conLSH that uses Spaced context based Locality Sensitive Hashing. With multiple spaced patterns, S-conLSH facilitates a gapped mapping of noisy long reads to the corresponding target locations of a reference genome. We have examined the performance of the proposed method on 5 different real and simulated datasets. S-conLSH is at least 2 times faster than the recently developed method lordFAST. It achieves a sensitivity of 99%, without using any traditional base-to-base alignment, on human simulated sequence data. By default, S-conLSH provides an alignment-free mapping in PAF format. However, it has an option of generating aligned output as SAM-file, if it is required for any downstream processing.

Conclusions

S-conLSH is one of the first alignment-free reference genome mapping tools achieving a high level of sensitivity. The spaced-context is especially suitable for extracting distant similarities. The variable-length spaced-seeds or patterns add flexibility to the proposed algorithm by introducing gapped mapping of the noisy long reads. Therefore, S-conLSH may be considered as a prominent direction towards alignment-free sequence analysis.

Background

Single molecule real time (SMRT) sequencing developed by Pacific Biosciences [1] and Oxford nanopore technologies [2] have started to replace previous short length next generation sequencing (NGS) technologies. These new technologies have enabled us to address many unsolved problems regarding genetic variations. With the increase in read length to around 20 KB [3], SMRT reads can be used to resolve ambiguities in read mapping caused by repetitive regions. Low GC bias and the ability to detect DNA methylation [1] from native DNA made SMRT data appealing for many real life applications. However, the high sequencing error rate of 13-15% per base [3] poses a real challenge in sequence analysis. Specialized methods like BWA-MEM [4], BLASR [5], rHAT [6], Minimap2 [7], lordFAST [8], etc., have been designed to align noisy long reads back to the respective reference genomes. BLASR [5] clusters the matched words from the reads and genome after indexing using suffix arrays or BWT-FM [9]. It uses a probability-based error optimization technique to find the alignment. BWA-MEM [4], originally designed for short read mapping, has been extended for PacBio and Oxford nanopore reads (with option -x pacbio and -x ont2d respectively) by efficient seeding and chaining of short exact matches. However, both methods are too slow to achieve a desired level of sensitivity [6]. This issue was addressed by rHAT [6] using a regional hash table where windows from the reference genome with the highest k-mer matches are chosen as candidate sites for further extension using a direct acyclic graph. Unfortunately, this method has a large memory footprint if used with the default word length of $k=13$, and it fails to accommodate longer k-mers to resolve repeats. Minimap2 [7], a recently developed method, uses concave gap cost, efficient chaining and fast implementation using SSE or NEON instructions to align reads with high sensitivity and speed. Another new method lordFAST [8] has been introduced to align PacBio’s continuous long reads with improved accuracy. MUMmer4 [10], a versatile genome alignment system, also has an option for PacBio read alignment (-l 15 -c 31), although it is less sensitive and accurate than the specialized aligners.

However, all the above mentioned methods come with large computational costs. Here, time and memory consumption are dominated by the alignment overhead. On top of that, alignment algorithms are often unable to correctly align distant homologs in the “twilight zone” with 20–35% sequence identity, as such weak similarities are difficult to distinguish from random similarities. For these reasons, alignment-free methods have become popular in recent years. See [11,12,13] for recent review papers and [14] for a systematic evaluation of these approaches. An alignment-free method Minimap [15] has been developed in 2016 for mapping of reads to the appropriate positions in the reference genome. Minimap groups approximate colinear hits using single linkage clustering to map the reads. However, Minimap suffers from low specificity. In this article, a new alignment-free method called S-conLSH has been proposed to overcome the above mentioned problems. Being suitable for low conserved areas and less computationally expensive, S-conLSH is sensitive as well as very fast at the same time.

A large proportion of the sequencing errors in SMRT data are indels rather than mismatches [3]. This makes it even more complicated to differentiate genomic variations from sequence errors. To resolve this issue, a concept of ‘context-based’ Locality Sensitive Hashing (conLSH) has been introduced by Chakraborty and Bandyopadhyay [16]. Locality Sensitive Hashing (LSH) [17, 18] has been successfully applied in many real life-science applications, ranging from genome comparison [19, 20] to large scale genome assembly [21]. In LSH, points close to each other in the feature space are hashed into localized slots of the hash table. However, in practice, the neighborhood or context of an object plays a key role in measuring its similarity to another object. Chakraborty and Bandyopadhyay [16] have shown that contexts of symbols (a base in reference to DNA) are important to decide the closeness of strings. They proposed conLSH to group sequences in localized slots of the hash table if they share a common context. However, a match for the entire context is a stringent criterion, considering the error profile of SMRT data. Even a mismatch or indel of length one, caused by a sequencing error, may mislead the aligner.

Therefore, to address this problem, an idea of spaced-context is introduced in this article. Unlike conLSH [16] which produces base-level alignments of the sequences using Sparse Dynamic Programming (SDP) based algorithm, the proposed method (S-conLSH) is an alignment-free tool. It employs multiple spaced-seeds or patterns to find gapped mappings of noisy SMRT reads to reference genomes. The spaced-seeds are strings of 0’s and 1’s where ‘1’ represents the match position and ‘0’ denotes don’t care position where matching in the symbols is not mandatory. The substring formed by extracting the symbols corresponding to the ‘1’ positions in the pattern is defined as the spaced-context of a sequence. Therefore, a spaced-context can minimize the effect of erroneous bases and, thereby, enhances the quality of mapping because it does not check all the bases for a match. This differentiates the proposed method from conLSH which looks into the entire context to compute the hash values.

A pattern-based approach was originally proposed by [22] when they developed PatternHunter, a fast and sensitive homology search tool. Later, multiple patterns or “spaced seeds” were proposed by the same authors [23]. Efficient algorithms to find optimal sets of patterns have been introduced by [24] and [25]. A fast alignment-free sequence comparison method using multiple spaced seeds has been described in [26], see also [27] and [28].

The algorithm, S-conLSH, described in this article is an alignment-free tool designed for mapping of noisy and long reads to the reference genome. The concept of Spaced-context has been elaborated in the Methods section along with a description of the proposed algorithm.

Results

Six different real and simulated datasets of E.coli, A.thaliana, O.sativa, S.cerevisiae and H.sapiens have been used to benchmark the performance of S-conLSH in comparison to other state-of-the art aligners, viz., Minimap2 [7], lordFAST [8], Minimap [15], conLSH [16] and MUMmer4 [10]. All these methods are executed in a setting designed for PacBio read alignment (see Table 1). The default parameter settings used for S-conLSH are $K=2$, context size ($2\lambda +1$) $=7$, $L=2$ and $z=5$ (refer to the Methods section of this article for details of the S-conLSH parameters). The two patterns used in our experiment in the default set up of $L=2$ are 11111110011111110000 and 111111100000111111100. The datasets used in the experiment have been summarized in Table 2. For the sake of simplicity, the results have been demonstrated by executing different aligners in a single thread. Please refer to the Tables [S-1] to [S-5] of Additional file 1: Note 1 for detailed review of their performance in multi-threaded mode.

Table 1 Parameter settings and commands used by different methods for mapping of PacBio SMRT reads

The aligner, rHAT [6] has been excluded from the study, as it has been reported to malfunction in certain scenarios [7]. The PacBio read alignment module of BWA-MEM [4] has been replaced by Minimap2, as it retains all the main features of BWA-MEM, while being  50$\times$faster and more accurate. Therefore, the results of BWA-MEM are not shown separately in the tables. Moreover, BLASR [5] has also not been used in the comparative study, as Minimap2 and lordFAST have been found to outperform it in all respects.

Table 2 The summary of real and simulated datasets used in the experiment along with the corresponding reference genome links

By default, S-conLSH produces output in pairwise read mapping format (PAF) ( [15]). There are scripts available to convert the popular SAM [29] alignment formats to PAF ( [15]). If a base-to-base alignment is requested, S-conLSH provides an option (--align 1), where the target locations are aligned using ksw alignment library (https://github.com/attractivechaos/klib) to produce the SAM file. The entire experiment has been conducted on an Intel Core i7-6200U CPU @ 2.30 GHz $\times$ 16(cores), 64-bit machine with 32GB RAM.

The results demonstrated in this article are organized into three categories: (1) performance on simulated datasets, (2) study on real PacBio reads, and (3) Robustness of S-conLSH for different parameter settings.

Experiment on simulated dataset

To study the accuracy of SMRT read mapping, a total of 146,932 noisy long reads have been simulated from hg38 human genome using PBSIM [30] command “pbsim --data-type CLR --depth 1 --length-min 1 --length-max 200000 --seed 0 --sample-fastq real.fastq hg38.fa. The error profile has been sampled from three real human PacBio RS II P5/C3 reads listed below, concatenated as real.fastq.

• m130929_024849_42213_c100518541*_s1_p0.1.subreads.fastq

• m130929_024849_42213_c100518541*_s1_p0.2.subreads.fastq

• m130929_024849_42213_c100518541*_s1_p0.3.subreads.fastq

The simulated reads from 5 different Human chromosomes are used to test the performance of S-conLSH in comparison to the other standard aligners. The sensitivity and precision have been computed based on the ground truth as obtained from the .maf files of PBSIM. A read is considered to be mapped correctly (as defined by [8]) if (1) it gets mapped to the correct chromosome and strand; and (2) the target subsequence of reference genome where the read maps to, must overlap with the true mapping by at least 1bp. The sensitivity is measured as a fraction of correctly mapped reads out of the total number of reads. Precision is defined, in the same way, as the fraction of correctly mapped reads out of the total number of mapped reads.

Table 3 Comparative study of the number of correct mappings, sensitivity, precision, and running time by different methods, Minimap, Minimap2, lordFAST, MUMmer4, conLSH, and S-conLSH, for a total of 146,932 reads simulated from five different human chromosomes

Table 3 summarizes the number of correct mappings, sensitivity, precision, and running time by different methods, Minimap, Minimap2, lordFAST, MUMmer4, conLSH, and S-conLSH for a total of 146,932 reads simulated from five different human chromosomes. The number of reads extracted from each chromosome is listed in Table 3. The result shows that S-conLSH produces the highest number of correct mappings among all five aligners for different chromosomes of Human-sim dataset. S-conLSH maps 32,111 reads out of total 32,290 reads of Chr#1, among which 31,964 mappings are found to be correct when compared with the ground truth. Minimap2 is the second highest in producing the correct mappings in this case. A similar scenario has been generally observed for the four other chromosomes as well. It is clear that Minimap2 always aligns all the reads to some location in the reference genome, but produces more incorrect mappings when compared to S-conLSH. Evidently, S-conLSH provides the highest sensitivity for all the chromosomes considered. Minimap, on the other hand, exhibits higher precision but lower sensitivity as it leaves a large number of reads unaligned. The number of unaligned reads by Minimap increases for large and complicated real datasets.

As can be seen, S-conLSH takes 38 CPU seconds to map the reads of chromosome 1, which is slightly slower than Minimap. The speed of Minimap is achieved as it maps a smaller number of reads compared to other aligners. Interestingly, S-conLSH has been found to have smaller mapping time than all the remaining algorithms, while having the maximum number of correctly mapped reads. As there was no separate indexing and aligning time available for MUMmer4, the total time is mentioned as “Mapping time”. MUMmer4 has been found to consume a large amount of time to achieve a desired level of sensitivity. It is evident from Table 3 that the indexing time is quite low for both Minimap and Minimap2. Indexing time for S-conLSH is relatively higher, though it is much smaller as compared to lordFAST. Here, it may be noted that indexing is performed only once for a given reference genome, while the read mapping will need to be performed every time a different individual is sequenced. The compressed and memory-efficient B-tree indexing of conLSH makes it the fastest in processing of reference genomes. However, the mapping time of conLSH is large as it performs base-to-base alignments using Sparse Dynamic Programming. The stringent ungapped matching requirement of the aligner over the entire context of the sequences results in lower sensitivity, after lordFAST. The proposed alignment-free tool, S-conLSH, has been found to be useful in such cases as it obtains the gapped mapping of the noisy reads using spaced-contexts.

While this section reports results of single-threaded execution, the Tables [S-1]–[S-5] of the Additional file 1: Note 1 exhibit the performance boost-up of S-conLSH in multi-threaded systems. S-conLSH achieves more than $50\%$ reduction in mapping time when run in 4 concurrent threads over the single-threaded version of itself. The indexing time of S-conLSH also improves with higher degree of parallelism and becomes comparable with that of Minimap2 when the number of threads is equal to 8. Moreover, this performance achievement comes with almost no additional burden of memory requirement. Please refer to Additional file 1: Note 1 for a detailed report.

Experiment on real PacBio datasets

This section demonstrates the performance of S-conLSH in comparison to other state-of-the-art aligners on five different SMRT datasets of E.coli-real, A.thaliana-real, O.sativa-real, S.cerevisiae-real and H.sapiens-real (refer Table 2 for details). A comparative study of running time, percentage of reads aligned and coverage by different aligners has been detailed in Table 4 for real human SMRT subread named m130929_024849_42213_c100518541* _s1_p0.1.subreads.fastq consisting of 23,235 reads . Results on MUMmer4 are excluded since it takes inordinately long.

Table 4 Comparative study of running time, percentage of reads aligned and coverage by different aligners for H.sapiens-real SMRT dataset of 23,235 reads

Table 5 Comparative study of running time, percentage of reads aligned by different aligners for four datasets of E. coli-real, A. thaliana-real, O. sativa-real and S. cerevisiae-real

It can be seen that S-conLSH provides the highest coverage value among the five standard methods used in the experiment. Minimap does not have any coverage statistics, as it is unable to produce alignment as SAM file. The performance in terms of indexing and mapping time, as shown in Table 4, is similar to that has already been observed for simulated datasets. The percentage of read alignment is the highest by Minimap2 and lordFAST. This is similar to the scenario obtained on simulated datasets where Minimap2 and lordFAST align all the reads against the reference genome, even though it may contain some incorrect mappings. S-conLSH, on the other hand, has a mapping ratio of $99.9\%$, which is lower than Minimap2 and lordFAST. This is due to the fact that S-conLSH gives higher priority to the mapping accuracy and it leaves a few reads unaligned if potential target locations are not found. S-conLSH has a higher memory footprint of about 13GB for indexing the entire human genome.

Similar results are observed for E.coli-real, A.thaliana-real, O.sativa-real, and S.cerevisiae-real real PacBio datasets as can be seen in Table 5. It is clear that S-conLSH is among the fastest in terms of mapping time, after Minimap. However, Minimap fails to align a good portion of the reads for large datasets like A.thaliana and S.cerevisiae. The aligner, conLSH, on the other hand, requires lower indexing time but higher mapping time to align a reasonable amount of reads to the reference genome. However, the alignment quality of conLSH is often compromised as studied in previous subsections. The percentage of reads mapped by S-conLSH is generally lower than Minimap2 and lordFAST, as it tries to ensure the best of the mapping quality. It is, however, difficult to measure the quality of read mappings for real datasets, as there is no ground truth available for such cases.

Robustness of S-conLSH for different parameter settings

An exhaustive experiment with different values of K, $\lambda$, L, and z has been carried out to study the robustness of the proposed method S-conLSH.

Table 6 Performance of conLSH with change of K, L, z, and $\lambda$ for real human SMRT dataset

Table 7 Performance of conLSH with change of z for chromosome 1 of H. sapiens-sim dataset consisting of 32,290 reads

Table 6 summarizes the study of indexing and mapping time along with the percentage of reads aligned for different values of S-conLSH parameters on real human SMRT dataset m130929_024849_42213_c100518541*_s1_p0.1.subreads.fastq. As can be observed the best performance (highest percentage of read mapping in minimum time) is achieved with the settings $K=2$, $(2\lambda +1)=7$, $L=2$ and $z=5$. The mapping time increases with L as it directly corresponds to the number of hash tables (one for each of the L different spaced-seeds) used to retrieve the target locations. Therefore the search becomes more rigorous as it considers all spaced-contexts obtained from L different patterns. This is, however, useful for highly sensitive applications at the expense of a few more seconds of mapping time. An efficient solution could be the use of S-conLSH with higher values of L, distributed over multiple concurrent threads. Please refer to Additional file 1: Note 2 regarding the performance of S-conLSH for different values of L in a multi-threaded system.

Indexing time, on the other hand, is proportional to the product $(2\lambda +1)K$. It seems that z has little effect on the performance of S-conLSH as the running time and the percentage of reads aligned mostly stay invariant with z. However, the parameter z is important to enhance the sensitivity of the method. This is reflected in Table 7 when studied on simulated reads. The highest number of correct mappings is obtained with the default settings (shown in italic) when $z=5$. The zeros in the spaced-seed help to find the distant similarities as it encompasses a larger portion of the sequence while the weight ($(2\lambda +1)K$) of the pattern remains the same. However, a very large value of z may degrade the accuracy as it joins unrelated contexts together. It is evident from Tables 6 and 7 that the performance of S-conLSH remains reasonably good irrespective of the variation of the parameter values. Therefore, it can be concluded that the algorithm S-conLSH is quite robust even though it requires some tuning of different parameters for the best performance.

Discussion and conclusions

S-conLSH is one of the first alignment-free reference genome mapping tools achieving a high level of sensitivity. Earlier, Minimap was designed to map reads against the reference genome without performing an actual base-to-base alignment. However, the low sensitivity of Minimap precluded its applications in real-life domains. Minimap2 is one of the best performing state-of-the-art alignment-based methods which provides an excellent balance of running time and sensitivity. The method described in this article, S-conLSH, has been observed to outperform Minimap2 in respect of sensitivity, precision, and mapping time. However, it has a longer indexing time and a higher memory footprint. Nevertheless, sequence indexing is a one-time affair, and memory is inexpensive nowadays.

The spaced-context in S-conLSH is especially suitable for extracting distant similarities. The variable-length spaced-seeds or patterns add flexibility to the proposed algorithm. Multiple patterns (with higher values of L) increase the sensitivity but at the cost of more time. Moreover, with the introduction of don’t care positions, the patterns become longer, thus providing better performance in resolving conflicts that occur due to the repetitive regions. The provision of rehashing for chimeric read alignment and reverse strand mapping make S-conLSH ideal for applications in the real-life sequence analysis pipeline.

A memory-efficient version of the S-conLSH can be developed in the future. The algorithm, at its current stage, can not conclude on the optimal selection of the patterns. A study on finding the optimal set of spaced-seeds can be carried out in future to improve the performance of the algorithm. Though the experiment demonstrated in this article is confined to the noisy long reads of PacBio datasets, it can be further extended on ONT reads as well. Finally, we would like to conclude with a strong expectation that the proposed method S-conLSH will draw the attention of the peers as one of the best performing reference mapping tools designed so far.

Methods

The algorithm S-conLSH for mapping noisy long reads to the reference genome essentially consists of two steps, reference genome indexing and read mapping. The complete workflow of S-conLSH is provided in Fig. 1 and the entire procedure is detailed below.

• 1 Reference Genome Indexing

  The reference genome is sliced into overlapping windows, and these windows are hashed into hash tables using suitably designed S-conLSH functions (see Definition 5) as shown in Fig. 1. S-conLSH uses two hash tables ‘$h\_index$’ and ‘Hashtab’. An entry in $h\_index$ has two fields (f, n): f stores an offset to the table Hashtab, where sequences are clustered according to their hash values, and n is the total number of sequences hashed at a particular value. Therefore, $Hashtab[h\_index[$x].f] to $Hashtab[h\_index[$x$].f+h\_index[$x].n] are the sequences hashed at value x.

• 2 Read mapping

  For each noisy long read, S-conLSH utilizes the same hash function for computing the hash values and retrieves sequences of the reference genome that are hashed in the same position as the read. Finally, the locations of the sequences with the highest hits are chained and reported as an alignment-free mapping of the query read (see Fig. 1).

By default, S-conLSH provides alignment-free mappings of the SMRT reads to the reference genome. If a base level alignment is required, S-conLSH provides an option (--align 1) to generate alignment in SAM format using ksw library (https://github.com/attractivechaos/klib). Some key aspects of S-conLSH are detailed in the following subsections.

Fig. 1

A schematic workflow of indexing and mapping using S-conLSH

Context based locality sensitive hashing

Locality Sensitive Hashing [17, 18] is an approximate near-neighbor search algorithm, where the points having a smaller distance in the feature space, will have a higher probability of making a collision. Under this assumption, a query is compared only to the objects having the same hash value, rather than to all the items in the database. This makes the algorithm work in sublinear time. In the definitions below, we use the following notations:

For a string x of length d over some set $\Sigma$ of symbols and $1\le i\le j \le d$, x[i] denotes the ith symbol of x, and x[i..j] denotes the (contiguous) substring of x from position i to position j. If ${\mathcal {H}}$ is a finite set of functions defined on some set X, for any $h\in {\mathcal {H}}$, randomly drawn with uniform probability, and $x,y\in X$, $Pr_{{\mathcal {H}}}[h(x)=h(y)]$ denotes the probability that $h(x) = h(y)$.

The definition of Locality Sensitive Hashing as introduced in [17, 18] is given below:

Definition 1

(Locality Sensitive Hashing) [17, 18] Let (X, D) be a metric space, let ${\mathcal {H}}$ be a family of hash functions mapping X to some set U, and let $R, c, P_1, P_2$ be real numbers with $c>1$ and $0\le P_2 < P_1 \le 1$. ${\mathcal {H}}$ is said to be $(R,cR,P_1,P_2)$-sensitive if for any $x,y\in X$ and $h\in {\mathcal {H}}$

• $\begin{aligned}Pr_{{\mathcal {H}}}[h(x)=h(y)]\ge P_1\end{aligned}$ whenever $D(x,y)\le R$, and

• $\begin{aligned} Pr_{{\mathcal {H}}}[h(x)=h(y)]\le P_2 \end{aligned}$ whenever $D(x,y)\ge cR$.

To illustrate the concept of locality sensitive hashing for DNA sequences, let us consider a finite set $\Sigma =\{A,T,C,G\}$ called the alphabet, together with an integer $d>0$. Let X be the set of all length-d words over $\Sigma$, endowed with the Hamming distance, and let U be the alphabet $\Sigma$. For $1\le i \le d$, let the function $h_i:X\rightarrow U$ be defined by $h_i(x) = x[i]$, $\forall x\in X$. Next, let R and cR be real numbers with $c>1$ and $0\le R < cR \le d$, and define $P_1 = \frac{d-R}{d}$ and $P_2 = \frac{d - cR}{d}$. Then the set ${\mathcal {H}}=\{h_i: 1\le i\le d$} is $(R,cR,P_1,P_2)$-sensitive. To see this, observe that for any two words $p,q \in X$, the probability $Pr_{{\mathcal {H}}}[h(p)=h(q)]$ is same as the fraction of positions i with $p[i] = q[i]$. Therefore,

$$\begin{aligned} Pr_{{\mathcal {H}}}[h(p)=h(q)] = \frac{d - D(p,q)}{d} \ge \frac{d-R}{d} = P_1 \end{aligned}$$

if $D(p,q) \le R$, and

$$\begin{aligned} Pr_{{\mathcal {H}}}[h(p)=h(q)] = \frac{d - D(p,q)}{d} \le \frac{d-cR}{d} = P_2 \end{aligned}$$

if $cR \le D(p,q)$.

Therefore, $P_1>P_2$ as $cR>R$. This proves that the family of hash functions ${\mathcal {H}}=\{h_i: 1\le i\le d$} is locality sensitive.

In biological applications, it is often useful to consider the local context of sequence positions and to consider matching subwords, as shown in conLSH [16]. It groups similar sequences in the localized slots of the hash tables considering the neighborhoods or contexts of the data points. A context in connection to sequence analysis can be formally defined as:

Definition 2

(Context) Let $x:(x_1x_2\dots x_d)$ be a sequence of length d. A context at the i-th position of x, for $i \in \{\lambda +1, \dots ,d-\lambda \}$, is a subsequence $x[i-\lambda \dots i\dots i+\lambda ]$ of length $2\lambda +1$, formed by taking $\lambda$ characters from each of the right and left sides of x[i]. Here, $\lambda$ is a positive constant termed the context factor.

To define context based locality sensitive hashing, the above example is generalized such that, for a given subword length $(2\lambda +1)< d$, each hash function in ${\mathcal {H}}$ will map words containing the same length-$(2\lambda +1)$ subwords at some position to the same bucket in U. The subword length $(2\lambda +1)$ is called the context size, where $\lambda$ is the context factor.

Definition 3

(Context based Locality Sensitive Hashing (conLSH)) Let $\Sigma$ be a set called the alphabet. Let $\lambda$ and d be integers with $(2\lambda +1) < d$. Let X be the set of all length-d words over $\Sigma$ and U be the set of all length-$(2\lambda +1)$ words over $\Sigma$. For $R,cR, P_1$, and $P_2$ as above, a $(R,cR,P_1,P_2)$-sensitive family ${\mathcal {H}}$ of functions mapping X to U is called $(R,cR,\lambda , P_1,P_2)$-sensitive, if for each $h\in {\mathcal {H}}$, there are positions $i_h$ and $j_h$ with $\lambda +1 \le i_h, j_h \le d - \lambda$ such that for all $p,q\in X$ one has $h(p) = h(q)$ whenever

$$\begin{aligned} p[i_h-\lambda \dots i_h\dots i_h+\lambda ] = q[j_h-\lambda \dots j_h\dots j_h+\lambda ] \end{aligned}$$

holds.

Gapped read mapping using spaced-context based locality sensitive hashing

The proposed method S-conLSH, uses spaced-seeds or patterns of 0’s and 1’s in connection with S-conLSH function. For a pattern ${\mathcal {P}}$, the spaced-context of a DNA sequence can be defined as:

Definition 4

(Spaced-context) Let ${\mathcal {P}}$ be a binary string or pattern of length $\ell$, where ‘1’ represents match position and ‘0’ represents don’t-care position. Let $\ell _w$ denote the weight of ${\mathcal {P}}$ which is equal to the number of ‘1’s in the pattern. Evidently, $\ell _w \le \ell$. Let x be a sequence of length d over alphabet $\{A,T,G,C\}$ such that $\ell \le d$. Then, a string sw over $\{A,T,G,C\}$ of length $\ell _w$ is called a spaced-context of x with respect to ${\mathcal {P}}$, if $sw[i] = x[j]$ holds if and only if ${\mathcal {P}}[j]=1$, where $i\le j$, $1\le i \le \ell _w$ and $1\le j\le \ell$.

Sequences sharing a similar spaced-context with respect to a pre-defined pattern ${\mathcal {P}}$, are hashed together in S-conLSH.

The concept of gap-amplification is used in locality sensitive hashing to ensure that the dissimilar items are well separated from the similar ones. To do this, gap between the probability values $P_1$ and $P_2$ needs to be increased. This is achieved by choosing L different hash functions, ${g_1,g_2,\dots , g_L}$, such that $g_j$ is the concatenation of K randomly chosen hash functions from ${\mathcal {H}}$, i.e., $g_j=(h_{1,j},h_{2,j}, \dots , h_{K,j})$, for $1 \le j \le L$. This procedure is known as “gap amplification” and K is called the “concatenation factor” [18]. For every hash function $g_j$, $1 \le j \le L$, there is a pattern ${\mathcal {P}}_j$ associated with it. The spaced-context based Locality Sensitive Hashing is now defined as follows:

Definition 5

(Spaced-context based Locality Sensitive Hashing (S-conLSH)) Let $sw_j(x)$ be the spaced-context of sequence x with respect to the binary pattern ${\mathcal {P}}_j$ of length $\ell$, $1 \le j \le L$. Let ${\mathcal {P}}_j$ be defined by the regular expression $\left( 0^*(1)^{(2\lambda +1)}\right) ^K0^*$. Therefore, the weight of ${\mathcal {P}}_j$, i.e., $\ell _w=(2\lambda +1)K$. The maximum value of $\ell$ would be $(2\lambda +1)K+z(K+1)$ assuming that at most z zeros are present between two successive contexts of 1’s in ${\mathcal {P}}_j$, where $z\ge 0$ is an integer parameter. Let d be an integer with $\ell \le d$, X be the set of all length-d words over $\Sigma$ and U be the set of all length-$\ell _w$ words over $\Sigma$. For $R,cR, P_1$, $P_2$, and $\lambda$ as introduced in Definition 3, a $(R,cR,\lambda , P_1,P_2)$-sensitive hash function $g_j=(h_{1,j},h_{2,j}, \dots , h_{K,j})$, where $h_{i,j}\in {\mathcal {H}}, 1\le i\le K$, mapping X to U is called $(R,cR,\lambda ,z, P_1,P_2)$-sensitive, if for any $p,q\in X$ one has $g_j(p) = g_j(q)$ whenever $sw_j(p)=sw_j(q)$ holds with respect to the pattern ${\mathcal {P}}_j$.

Fig. 2

A schematic illustration of gapped-mapping using S-conLSH. a Multiple patterns having context size $=3$ and $K=2$. b, c Hashing of the strings “ATTCGGTAA” and “TTCTAAGTA” respectively using different patterns. d Final hash table and gapped-mapping of the two strings due to the collision at “TTCTAA”

Therefore, instead of restricting similarity over the $(2\lambda +1)K$ consecutive bases as was done for conLSH [16], S-conLSH incorporates greater flexibility by checking only the positions which correspond to a 1 in the pattern. For example, the binary string “011100111” is a pattern for a system having $K=2$, $z=2$ and context size $(2\lambda +1)=3$. The hash value or the spaced-context of the string “ATTCGGTAA” for the above pattern will be “TTCTAA” (see Fig. 2(b)). In S-conLSH, noisy long reads are hashed using L functions corresponding to L different patterns generated using Algorithm 1. Multiple pattern based functions enable gapped-mapping of the reads as illustrated in Fig. 2. Consider a scenario of two patterns ${\mathcal {P}}_1=$“011100111” and ${\mathcal {P}}_2=$“111111” having context size $=3$, $L=2$ and $K=2$. The string $p=$“ATTCGGTAA” generates two hash values $sw_1(p)=$“TTCTAA” and $sw_2(p)=$“ATTCGG” for the patterns ${\mathcal {P}}_1$ and ${\mathcal {P}}_2$ respectively (see Fig. 2b). Similarly, $sw_1(q)=$“TCTGTA” and $sw_2(q)=$“TTCTAA” are the hash values for string $q=$“TTCTAAGTA” (Fig. 2c). As shown in the hash table of Fig. 2d, the two strings collide to the same bucket of the hash table due to the common hash value “TTCTAA”. This results in mapping with three gaps or indels, corresponding to the three 0’s of “011100111”, in the second string. This gapped-mapping is a powerful feature of S-conLSH which is quite uncommon in standard spaced-seed based methods (refer Additional file 1: Note 3 for details).

To obtain an integer hash value from the Spaced-context, an encoding function $f: S\mapsto \{0,1,\dots ,(4^{K(2\lambda +1)}-1)\}$, $f(sw)=\sum _{i=1}^{(2\lambda +1)K} f(sw[i])\times 4^{(2\lambda +1)K-i}$, $\forall sw\in S$, has been defined assuming $f(A)=0$, $f(C)=1$, $f(G)=2$ and $f(T)=3$, where S is the set of all spaced-contexts of length $(2\lambda +1)K$ defined over the alphabet $\{A,T,C,G\}$. A pattern produces hash values of length equal to its weight. Keeping the weight same, the pattern length is increased in S-conLSH by introducing don’t care positions (or, zeros). This allows S-conLSH to look at a larger portion of the sequences without increasing the computational overhead. Consequently, S-conLSH is able to find distant homologs that might otherwise be overlooked. Not only that, it provides better sensitivity in resolving repeats because of the consideration of the neighborhood (or, contexts) when measuring the similarity between the sequences. S-conLSH has a provision of split mapping for chimeric reads as follows. If a read fails to get associated with end-to-end mapping, it is split into a series of non-overlapping segments and re-hashed to find target location(s) for each segment.