Fast randomized approximate string matching with succinct hash data structures

Definitions

Throughout this paper we will work with the alphabet Σ _{DN A} = {A, C, G, T, N, $} (with N and $ being the undefined base and the contig end-marker, respectively) which, in practice, will be encoded in Σ _{DN A'} = {0, 1, 2, 3} assigning a numerical value in Σ _{DN A'} to N and $ characters. The size of our alphabet is, therefore, σ = |Σ _{DN A'} | = 4, while with n, m, and w we will denote the reference length, the pattern length, and the (fixed) size of a computer memory-word (i.e. the number σ ^w − 1 is assumed to fit in a memory word), respectively. As hash functions we will use functions of the form h : Σ ^m → Σ ^w mapping length-m Σ-strings to length-w Σ-strings. If necessary, we will use the symbol ${}_w^{m}h$ instead of h when we need to be clear on h's domain and codomain sizes. Given a string $P\in \sum ^m$, the value $h\left(P\right)\in \sum ^m$ will be also dubbed the fingerprint of P (in Σ ^w ). With $T\in {\sum }^n$ we will denote the text that we want to index using our data structure. $T_i^j$ will denote $T\left[i,...,i+j-1\right]$, i.e. the length-j prefix of the i-th suffix of T. A hash data structure H for the text T with hash function h, will be a set of ordered pairs (an index) such that $H=\left\{⟨h\left(T_i^m\right),i⟩:0\le i\le n-m\right\}$, that can be used to store and retrieve the positions of length-m substrings of T (m is therefore fixed once the index is built). A lookup operation on the hash H given the fingerprint h(P), will consist in the retrieval of all the positions $0\le i<n$ such that $⟨h\left(P\right),i⟩\in H$ and cases where $⟨h\left(P\right),i⟩\in H$ but $T_i^m\ne P$ will be referred to as false positives.

The symbol ⊕ represents the exclusive OR (XOR) bitwise operator. a ⊕ b where $a,b\in \sum$, will indicate the bitwise XOR among the bits of the binary encoding of a and b. Analogously, x ⊕ y, where $x,y\in {\sum }^m$ will indicate the bitwise XOR operation among the bits of the binary encoding of the two words x and y and V will denote the bitwise OR operator. $d_H\left(x,y\right)$ is the Hamming distance between $x,y\in {\sum }^m$. We point out that the Hamming distance is computed between characters in Σ, and not between the bits of the binary encoding of each of them. Patterns and fingerprints are viewed as bit vectors only when computing bitwise operations such as OR, AND, and XOR.

Succinct representation of hash indexes

We begin by introducing the technique allowing succinct representation of hash indexes. The central property of the class of hash functions we are going to use is given by the following definition:

With the following theorem we introduce the hash function used in the rest of our work and in the implementation of our structure:

is a de Bruijn hash function.

A detailed proof of this theorem is given in Additional file 1. Given a de Bruijn hash function ${}_w^{m}h:{\sum }^m\to {\sum }^w$ we can "extend" it to another de Bruijn hash function ${}_{n-m+w}^{n}h:\sum ^n\to \sum ^{n-m+w}$, operating on input strings of length n greater than or equal to m, as follows:

for every 0 ≤ i ≤ n − m.

It is easy to show that a function enjoying the property in Definition 2 is a homomorphism on de Bruijn graphs (having as sets of nodes Σ ^m and Σ ^w , respectively). Since ${}_w^{m}h$ univocally determines ${}_{n-m+w}^{n}h$ and the two functions coincide on the common part Σ ^m of their domain, in what follows we will simply use the symbol h to indicate both.

From Definitions 1 and 2 we can immediately derive the following important property:

Lemma 1 If h is a de Bruijn hash function, n ≥ m, and P ∈ Σ ^m occurs in T ∈ Σ ⁿ at position i, then h(P) occurs in h(T) at position i. The opposite implication does not (always) hold; we will refer to cases of the latter kind as false positives.

On the ground of Lemma 1 we can propose, differently from standard approaches in the literature, to build an index over the hash value of the text, instead of building it over the text. This can be done while preserving our ability to locate substrings in the text, since we can simply turn our task into that of locating fingerprints in the hash of the text T. We call dB-hash the data structure obtained with this technique. Notice that the underlying hash data structure is simulated by searching the occurrences of h(P) in h(T) during a lookup operation, so the algorithm is transparent to the particular indexing technique used.

Search algorithm

The core of our searching procedure is based on the algorithm rNA (Vezzi et al. [5], Policriti et al. [6]), a hash-based randomized numerical aligner based on the concept of Hamming-aware hash functions (see [5] and [6] for more details). Hamming-aware hash functions are particular hash functions capable to "squeeze" the Hamming ball of radius k around a pattern P to a Hamming ball of the same radius around the hash value of P. This feature allows to search the entire Hamming ball around P much more efficiently. The following theorem holds:

for every $P,P\prime \in {\sum }_{DNA\prime }^m$

See Additional file 1 for a detailed proof of this theorem. Since h_⊕ is a de Bruijn and Hamming aware hash function, we can use it to build our structure and adapt the rNA algorithm to it. We call dB-rNA the new version of the rNA algorithm adapted to the dB hash data structure. More in detail, the Hamming-awareness property of h_⊕ guarantees that, given a pattern P to be searched in the index, the set {h_⊕(P') : d _H (P, P') ≤ k} is small--$\mathcal{ O}$((2σ − 2) ^k w ^k ) = $\mathcal{ O}$(6 ^k w ^k ) elements in our application--and can be computed in time proportional to its size. Notice that, with a generic hash function h, only the trivial upper bound $\mathcal{ O}$((σ − 1) ^k m ^k ) can be given to the size of this set since each different P' such that d _H (P, P') ≤ k could give rise to a distinct fingerprint. Our proposed algorithm is almost the same as the one described in [5, 6], the only difference being that the underlying data structure is a dB-hash instead of a standard hash. Briefly, the search proceeds in 3 steps. For each pattern P:

1 h_⊕(P) is computed;

2 the index is searched for each element in the set $\left\{h\left(P\prime \right):d_H\left(P,P\prime \right)\le k\right\}$,

3 for each occurrence found, the text and the pattern P are compared to determine Hamming distance and discard false positives.

In practice, in our implementation we also split the pattern P in non-overlapping blocks before searching the index. With this strategy we reduce the maximum number of errors to be searched, improving the speed of the tool.

Complexity analysis

here, σ = 4 is the size of the alphabet Σ _{DN A'} . A fully formal proof of (an extended version of) this analysis can be found in [11].

Quality-aware strategy

Our tool implements a quality-aware heuristic that significantly improves search speed, at the price of a small loss in sensitivity. Briefly, we use base qualities to pick up only a small fraction of the elements from the Hamming ball centred on the hash h(B) of the searched block B, following the assumption that a high quality base is unlikely to be a miscall. Since in practice we divide the read in non-overlapping blocks and the heuristic affects only the searched block, with this strategy we lose only a small fraction of single variants like SNPs. More in detail, let $Q\in {\mathbb{N}}^m$ be (e.g.) the Phred quality (see [12]) string associated to the searched block B. We compute a hash value on Q using the following hash function:

q is a quality threshold (in our implementation we use q = 15). The values 0 and 3 have been chosen due to their binary representation (00 and 11, respectively). If $h_{\vee }\left(Q\right)\left[i\right]=3$, then during search we try to insert an error in position i of h_⊕ (B) since at least one of the bases used to compute h_⊕ (B)[i] has a low-quality. The quality-aware strategy is then implemented as follows: let B be the block to be searched and Q its associated Phred quality string. A fingerprint f (representing a block at distance at most k from B) is searched in the structure if and only if $\left(f\oplus h_{\oplus }\left(B\right)\right)\vee h_{\vee }\left(Q\right)=h_{\vee }\left(Q\right)$, i.e. if f differs from h_⊕(B) only in positions corresponding to low quality bases. Since the number of low-quality base pairs in a read is typically low, this strategy allows to drastically reduce search space (which in practice leads approximately to a 10x speedup) if reliable qualities are available. In the results section we show, moreover, that this strategy has only a negligible impact on SNP detection and on the overall precision of our tool.

BW-ERNE: implementation details

We implemented our algorithm and data structure in the short reads aligner BW-ERNE (Burrows-Wheeler Extended Randomized Numerical alignEr), downloadable at http://erne.sourceforge.net. As hash function for our index we use h_⊕. Given a text $T\in {\sum }_{DNA}^n$, we calculate h_⊕ (T) ^{BW T} --the Burrows-Wheeler transform of h_⊕(T)--adding the necessary additional structures needed to perform backward search and to retrieve text positions from h_⊕ (T) ^{BW T} positions. BW-ERNE includes (from its predecessor ERNE) also a simple and fast strategy to allow a single indel in the alignment. This strategy does not affect running times and permits to correctly align a large fraction of short reads that come with indels (see Results section). It is well known that DNA is extremely difficult to compress and for this reason we choose not to introduce compression in our structure. Even if our index is not compressed, experiments show (see Section) that its memory requirements are similar or even smaller than those of other tools based on the FM index such as Bowtie [2], BWA [4] and SOAP2 [3]. Briefly, the structure is composed by three parts: the index, the plain text, and an auxiliary (standard) hash.

Memory footprint of the indexes

BW-ERNE significantly reduces the space requirements of ERNE, requiring approximately 8.7 times less space on the Vitis Vinifera genome (730 MB vs 64 GB) and approximately 4.4 times less space on the Human genome (4.3 GB vs 19 GB). The plot in Figure 1 compares the different indexes sizes on Vitis Vinifera and gives a clear idea of the most important difference between full-text and succinct indexes, the former requiring space proportional to n log n (n pointers to the text) while the latter is using an amount of space close to that necessary for the plain text. More-over, differences among BWT-based tools are minimal (few hundred MB) if compared to the gap between hash-based (ERNE) and BWT-based (Bowtie, BWA, SOAP2, BW-ERNE) aligners. Among the tested tools only Bowtie requires less space than our tool, even if in the dB-hash data structure we do not perform any data compression. This is due to the fact that DNA is a high-entropy string and, therefore, almost incompressible using standard techniques such as the ones implemented in the FM index. Finally, the horizontal green line, marking the size of the reference fasta file, gives an idea of how efficiently can succinct and compressed indexes represent their structures: in particular, Bowtie is able to reach a RAM memory footprint that is smaller than that of the reference file itself, while BWA and BW-ERNE require slightly more space. On the Human genome, BW-ERNE index requires only 4.3 GB of space to be stored and loaded during alignment. This is slightly more than the indexes of Bowtie and BWA (2.7 GB and 3.2 GB of RAM, respectively), but considerably less than ERNE's hash table, which requires 19 GB of space to be stored and loaded during alignment.

Simulated data with reliable base qualities

In this experiment we compared the tools on 5M of 100bp single-end reads with simulated base qualities. This dataset was generated from the Vitis Vinifera genome using the SimSeq simulator, adopting the built-in Illumina error model. The number of correctly mapped reads was estimated comparing the bam files generated by SimSeq and the aligners (reads mapping in multiple locations were evaluated on the unique reported alignment).

The plots in Figure 2 show that BW-ERNE is able to exploit at best quality information without losing accuracy with respect to ERNE while, at the same time, improving significantly its performances. BWERNE was executed with default settings and using 1 thread. The plot on the left hand side of Figure 2 shows that BW-ERNE was 2 times faster than Bowtie, and 4 times faster than BWA. This speed came with no penalties on the number of mapped reads, which was the highest among all tools, with BW-ERNE and ERNE aligning 97% of all reads, BWA 94%, Bowtie 90% and SOAP2 82%. Finally, the plot on the right hand side of Figure 2 shows the accuracy of the tools in terms of correctly mapped reads. The gap between mapped and correctly mapped reads is due to reads mapping in multiple locations, which were judged on the base of the unique reported alignment. The plot shows that ERNE and BW-ERNE were the most accurate tools, correctly aligning the highest number of reads.

Simulated data without reliable base qualitie

Validation of the quality-aware strategy in presence of SNPs

In order to assess the impact of our quality-aware strategy on SNP detection, we simulated (using SimSeq) 10M of 100bp single-end reads, using as reference the Human genome (hg19). After the simulation, each base was randomly substituted with probability 0.005 to simulate SNPs (which, from the aligner's point of view, are simply mismatches with high base-quality). This strategy allowed us to track reads containing SNPs, permitting a separate verification of the alignment's correctness for reads with and without this kind of mutations. BW-ERNE was executed twice on the mutated dataset: with default settings (quality-aware strategy enabled) and with the --sensitive option enabled (quality-aware strategy disabled). An alignment was considered correct if and only if both chromosome and strand coincided with those outputted by SimSeq and if the alignment's position was within 50 bases from the position outputted by Simseq (in order to account for indels and clipped bases). Reads with multiple alignments were judged on the basis of their unique reported alignment. Of the 10M simulated reads, 39.42% contained at least one SNP. BWERNE in sensitive mode (quality-aware strategy disabled) correctly aligned 87.80% of the reads without SNPs and 87.64% of the reads with SNPs, thus showing (as expected) no significant bias towards reads without SNPs (the 0.16% difference can be explained with the fact that reads with SNPs are inherently more difficult to align). BW-ERNE with the quality-aware strategy enabled correctly aligned 86.54% of the reads without SNPs and 85.25% of the reads with SNPs, thus showing only a slight bias towards reads without SNPs.

Real data -- Human genome

To conclude, the experiment on a real high-coverage Human Illumina library, allowed us to validate the results obtained on simulated reads and to assess the performances of BW-ERNE on large datasets. In this experiment we aligned 320M of 100bp pairedend reads, corresponding to a 10x coverage of the Human genome, downloaded from the 1000genomes project's database (top 160M reads in ftp://ftp. http://1000genomes.ebi.ac.uk/vol1/ftp/data/NA12878/sequence_read/SRR622457_1.filt.fastq.gz and top 160M reads in ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data/NA12878/sequence_read/SRR622457_2.filt.fastq.gz). BW-ERNE was executed with default settings (quality-aware strategy enabled) and using 4 threads. Our tool completed the alignment in 3 hours and 15 minutes, for an overall throughput of 98 millions of reads per hour. 15916364 reads (5% of all the reads) were automatically discarded by the builtin trimmer due to low base quality. Of the remaining 304083636 reads, 300717664 (98.9%) were successfully aligned to the reference and 3365972 (1.1%) were not found. Among the aligned reads, 284650317 (94.6%) were aligned in only one position and 16067347 (5.4%) in multiple positions.