Document text
ii
“main” — 2017/4/19 — 0:38 — page 1 — #1ii
i
ii
i-
xxxxxx/xxxxxx/xxxxxx
Advance Access Publication Date: Day Month Y ear
Manuscript Category
Subject Section
FMtree: A fast locating algorithm of FM-indexes for
genomic data
Haoyu Cheng1,2,3, Ming Wu1,2,3and Yun Xu1,2,3,
1School of Computer Science, University of Science and T echnology of China, Heifei, Anhui, 230027, China
2Key Laboratory on High Performance Computing, Anhui Province
2Collaborative Innovation Center of High Performance Computing, National University of Defense T echnology, Changsha, 410073,
China
T o whom correspondence should be addressed.
Associate Editor: XXXXXXX
Received on XXXXX; revised on XXXXX; accepted on XXXXX
Abstract
Motivation: As a fundamental task in bioinformatics, searching for massive short patterns over a long
text has been accelerated by various compressed full-text indexes. These indexes are able to provide
similar searching functionalities to classical indexes, e.g., suffix trees and suffix arrays, while requiring less
space. For genomic data, a well-known family of compressed full-text index, called FM-indexes, presents
unmatched performance in practice. One major drawback of FM-indexes is that their locating operations,
which report all occurrence positions of patterns in a given text, are particularly slow, especially for the
patterns with many occurrences.
Results: In this paper, we introduce a novel locating algorithm, FMtree, to fast retrieve all occurrence
positions of any pattern via FM-indexes. When searching for a pattern over a given text, FMtree organizes
the search space of the locating operation into a conceptual quadtree. As a result, multiple occurrence
positions of this pattern can be retrieved simultaneously by traversing the quadtree. Compared with the
existing locating algorithms, our tree-based algorithm reduces large numbers of redundant operations and
presents better data locality. Experimental results show that FMtree is one order of magnitude faster than
the state-of-the-art algorithms, and still memory-efficient.
Availability: FMtree is freely available at https://github.com/chhylp123/FMtree .
Contact: [email protected]
Supplementary information: Supplementary data are available online.
1 Introduction
The string matching problem is to identify the occurrence positions of
a short string P(called pattern ) in a given long string T(called text).
For genomic data, both the pattern Pand the textTare the sequences
over a small alphabet =fa;c;g;tg. In order to speed up the string
matching process, a well-known approach is indexed matching . It first
builds an index data structure for the text in advance, and then searches
for the pattern via the index. If the text is static, this approach is usually
much faster than online matching approach, which directly searches for
the pattern over the text. Many bioinformatics applications have adopted
the indexed string matching approach, such as read mapping (Langmeadand Salzberg, 2012; Li, 2013), genome assembly (Li, 2012; Simpson and
Durbin, 2012) and read error correction (Schulz et al., 2014).
A number of full-text indexes have been proposed for several decades.
Classical indexes like suffix arrays (Manber and Myers, 1993), efficiently
support two basic functions: count andlocate . Given a pattern Pand a
textT, the count function is to report the number of occurrences of P
inT, while the locate function is to retrieve all occurrence positions of
PinT. A serious problem of these classical indexes is that their space
usage is relatively large, especially for a very long text. To address this
problem, various compressed full-text indexes have been developed in
recent years (Ferragina et al., 2009). Generally, most of them are able to be
classified into three families: FM-indexes (Ferragina and Manzini, 2000;
Grabowski et al ., 2004), compressed suffix arrays (CSAs) (Grossi and
Vitter, 2005; Sadakane, 2003) and Lempel-Ziv compression based indexes
- 1arXiv:1704.04615v2 [cs.DS] 18 Apr 2017
ii
“main” — 2017/4/19 — 0:38 — page 2 — #2ii
i
ii
i2 Cheng et al.
(LZ-indexes) (Arroyuelo, 2006; Arroyuelo et al., 2012). These indexes are
designed to provide similar count and locate functionalities to classical
indexes, while requiring less space. For genomic data, the most efficient
family of compressed full-text indexes is FM-indexes (Gog et al., 2017;
Vyverman et al., 2012), especially when searching for short patterns (Hon
et al., 2004). Given a human genome with about 3.15 billion characters,
FM-indexes usually require less than 3GB RAM, while classical suffix
arrays require about 12GB RAM. In addition, thanks to the small alphabet
size of genomic data, the counting time of FM-indexes is comparable
to that of classical indexes (Deorowicz and Grabowski, 2013; Gog and
Petri, 2014). Thus, FM-indexes have become the essential data structure
in many bioinformatics algorithms (Langmead and Salzberg, 2012; Li,
2013; Marco-Sola et al., 2012).
The major bottleneck of FM-indexes is that their locating operation is
several orders of magnitude slower than that of classical indexes (Ferragina
et al., 2009). To reduce space usage, FM-indexes only save a small fraction
of text positions, called sampled positions , rather than all of them. When
locating a pattern via FM-indexes, the sampled positions of the pattern can
be directly retrieved, while the non-sampled positions have to be calculated
one-by-one exploiting the expensive LF-mapping operations (see details
in Section 2.2). For short patterns with many occurrence positions, FM-
indexes need to perform large numbers of LF-mapping operations until
all non-sampled positions have been obtained. Unfortunately, searching
for short patterns in a long text is an important task in bioinformatics
(Langmead and Salzberg, 2012; Ahmadi et al., 2012; Cheng et al., 2015;
Hach et al ., 2010; Liu et al ., 2015), and these short patterns are very
frequent in practice (Xin et al., 2016). In this case, the cost of locating
operations dominates the overall string matching time.
Unfortunately, although there are various studies about compressed
full-text indexes, only a small fraction of them focus on accelerating the
locating operations. González et al . propose locally compressed suffix
array (LCSA) (González and Navarro, 2007; González et al ., 2015)
to improve the data locality of the locating operations. For popular
compressed full-text indexes such as FM-indexes and CSAs, a serious
problem is that locating patterns via them results in many random memory
accesses. To solve this problem, LCSA directly compresses suffix array
exploiting the repetitions of suffix array. Since all occurrence positions
of a pattern are saved consecutively in suffix array, LCSA can obtain
these positions by decompressing consecutive elements of suffix array. In
this case, its memory accesses are highly local. However, compared with
other compressed full-text indexes, LCSA requires much more space when
indexing genomic data.
Besides, Ferragina et al . develop a distribution-aware algorithm
(Ferragina et al., 2013), which adjusts FM-indexes or CSAs according
to the distribution of the occurrence positions of query patterns. This
algorithm assumes that the distribution of the patterns’ occurrence
positions has been known in advance, so that it inclines to sample the
text positions which have high probability to be located. However, in most
cases, it is impossible to know this distribution during the index building
phase. Besides, the distribution-aware algorithm cannot achieve good
performance unless the distribution of the patterns’ occurrence positions
is very skewed. This requirement also limits the usage of this algorithm.
In theory, LZ-indexes are more efficient than FM-indexes and CSAs
when performing locating operations (Ferragina et al ., 2009). But for
genomic data, it is difficult to develop a highly optimized LZ-index
like the existing sophisticated implementations of FM-indexes in many
bioinformatics algorithms (Langmead and Salzberg, 2012; Li, 2013;
Marco-Sola et al ., 2012). Moreover, there does not exist a practical
implementation of LZ-indexes which is able to process large texts with
billions of characters.
Here we introduce a novel locating algorithm, FMtree, to significantly
accelerate the locating operations of FM-indexes for genomic data. Whenlocating a pattern via FM-indexes, the search space of locating operation
is organized into a quadtree. By utilizing this quadtree, FMtree is able
to calculate the non-sampled positions block-by-block, while current
algorithms have to calculate these positions one-by-one. Thus, our tree-
based locating algorithm is cache-friendly and avoids many unnecessary
operations. Another advantage of FMtree is that it can be applied
to any implementation of FM-indexes without modification. Overall
experimental results show that FMtree significantly outperforms previous
algorithms for genomic data.
2 Background
2.1 Definitions and notation
A stringSis a sequence of characters over the alphabet , and the size of
is ||. We letjSjdenote the length of S,S[i]denote thei-th character
ofS(0ijSj 1), andS[i;j]denote the substring that starts at
S[i]and ends atS[j]. Besides, consider a character sand a stringS,sS
denotes the concatenation of sandS. We also let sndenote a string of
lengthnsuch thats0is an empty string, s1=sandsn=ssn 1.
To solve the string matching problem, existing compressed full-text
indexes need to support the following two basic operations:
count( P,T): Return the number of occurrences of pattern Pin textT.
locate( P,T): Return all occurrence positions of pattern Pin textT.
For the convenience of further discussion, we assume that a special
character $ is at the end of text T, where $ is lexicographically smaller
than other characters in . Since we focus on genomic data, all characters
in textTand patternPbelong to =fa;c;g;tgexceptT[jTj 1] = $ .
2.2 Overview of FM-indexes
As a family of compressed full-text indexes, FM-indexes are first proposed
to emulate classical suffix arrays (Ferragina and Manzini, 2000). Given a
textT, its suffix array SA(Manber and Myers, 1993) saves the positions
of all suffixes of Tin lexicographic order. For a pattern P, it is obvious
that the positions of all suffixes prefixed by Pare saved consecutively in
an interval of SA, calledSA[sp;ep ]. In fact,SA[sp;ep ]consists of all
occurrence positions of PinT.
Compared with suffix arrays, FM-indexes provide similar searching
functionalities, while requiring less space. The critical data structure of a
FM-index is Burrows Wheeler Transform (Burrows and Wheeler, 1994)
ofT, calledBWT (T), which permutes the characters of Treversibly.
Conceptually, BWT (T)can be constructed in the following two steps:
Building a conceptual matrix M(T)including all cyclic rotations of
Tin lexicographic order. Each row in M(T)is a cyclic rotation of T.
LetBWT (T)be the last column of M(T).
BWT(T) $ acgtaacc a
a $acgtaac c
a acca$acg t
a cca$acgt a
a cgtaacca $
c a$acgtaa c
c ca$acgta a
c gtaacca$ a
g taacca$a c
t aacca$ac g Matrix M(T)
9
8
4
5
0
7
6
1
2
3SA
F
Fig. 1. An example of Burrows Wheeler Transform BWT (T)and suffix array SA for
stringT=“acgtaacca$".
ii
“main” — 2017/4/19 — 0:38 — page 3 — #3ii
i
ii
iFMtree 3
An example is presented in Fig. 1. For FM-indexes, there are two
important properties of M(T)andBWT (T):
LetFbe the first column of M(T),BWT (T)[i]precedesF[i]
in textT. It is obvious that the i-th row ofM(T)corresponds to
SA[i]. Specifically, BWT (T)[i] =T[SA[i] 1]ifSA[i]6= 0, and
BWT (T)[i] = $ ifSA[i] = 0 . In addition, F[i] =T[SA[i]].
For any character s2, thei-thsinFand thei-thsinBWT (T)
correspond to the same sinT. For example in Fig. 1, F[5]and
BWT (T)[1]are the firstcinFandBWT (T), respectively. In fact,
both of them correspond to T[7].
Based on the above properties, a core operation LF (Last-to-First
mapping or LF-mapping) of FM-indexes is defined as:
LF(l) =C[BWT (T)[l]] +rank BWT (T)[l](BWT (T);l) (1)
whereC[BWT (T)[l]]denotes the number of characters in Twhich are
smaller thanBWT (T)[l], andrank BWT (T)[l](BWT (T);l)is a rank
operation that reports the number of BWT (T)[l]inBWT (T)[0;l
1]. Essentially, LF operation scans text Tbackward, namely, F[l] =
T[SA[l]],F[LF[l]] =BWT (T)[l] =T[SA[l] 1]andSA[LF[l]] =
SA[l] 1. For instance, in Fig. 1, SA[0] = 9 ,LF[0] = 1 ,SA[LF[0]] =
SA[1] =SA[0] 1 = 8 .
Counting via FM-indexes. When counting a pattern Pin textT, it is
obvious that all M(T)’s rows prefixed by Pare saved consecutively in
an intervalM(T)[sp;ep ]. Thus, to answer count (P;T), FM-indexes first
perform the backward search algorithm to determine the range [sp;ep ].
After that, count (P;T) =ep sp+1. Algorithm 1 presents the procedure
of the backward search algorithm. More precisely, this algorithm searches
P[0;jPj 1]backward injPjsteps. Ini-th step (i=jPj 1;jPj
2;:::; 0), this algorithm updates [sp;ep ]such thatM(T)[sp;ep ]includes
allM(T)’s rows prefixed by P[i;jPj 1]. Note that in first step ( i=
jPj 1),[sp;ep ]is the range of a single character P[jPj 1], so that it can
be directly obtained via array C(line 2 in Algorithm 1). For any character
s2,C[s]saves the number of characters which are lexicographically
smaller thansinT. Finally, when i= 0, the range [sp;ep ]ofP[0;jPj 1]
is obtained and count (P;T) =ep sp+ 1. SinceM(T)[i]corresponds
toSA[i],SA[sp;ep ]actually saves all occurrence positions of PinT.
Locating via FM-indexes. As mentioned above, the backward search
algorithm of FM-indexes determines the range [sp;ep ]such that
SA[sp;ep ]consists of all occurrence positions of PinT. If the whole
Algorithm 1: backward_search( P)
Input: the patternP[0;jPj 1]
Output: theSArange [sp;ep ]ofP
1:i=jPj 1;s=P[i];i=i 1
2:sp=C[s];ep=C[s+ 1] 1
3:whilei>0 &&sp6epdo
4:s=P[i];i=i 1
5:sp=C[s] +rank s(BWT (T);sp)
6:ep=C[s] +rank s(BWT (T);ep+ 1) 1
7:end while
8:ifsp>ep then
9:return "not found"
10:else
11: return [sp;ep ]
12:end ifAlgorithm 2: locate(sp;ep )
Input: theSArange [sp;ep ]
Output: the position set Rconsists of all positions in SA[sp;ep ]
1:fori=sptoepdo
2:j=i;m= 0
3:whileSA[j]is not sampled do
4:j=LF(j),m=m+ 1
5:end while
6: AddSA[j] +mtoR
7:end for
8:returnR
SA is saved, locate (P;T)can be answered by retrieving SA[sp;ep ]
directly. However, for a long text, its SA is very space-consuming.
Thus, FM-indexes only save a fraction of positions in SA, called the
sampled positions. This strategy reduces the space usage of FM-indexes,
but of course comes at the expense of additional computational overhead.
Algorithm 2 presents the most practical locating algorithm of FM-indexes.
To obtain position SA[i](spiep), this algorithm first scans text T
backward by performing LF operation mtimes until a sampled position
SA[j]is reached. After that, SA[i] =SA[j] +m.
Implementations of FM-indexes. As we can see, both the counting
operation and the locating operation of FM-indexes can be reduced to
rank operations. To support rank operations, practical implementations of
FM-indexes break BWT (T)into small blocks. For the beginning line of
each block and each character s2, the rank value is precomputed and
saved. When calculating rank s(BWT (T);l), these implementations
first retrieve the precomputed rank value of the block which includes
BWT (T)[l], and then add the number of rest sin this block before l-
th line. Besides, the locating algorithm of FM-indexes requires a sampled
suffix arraySSA , which saves all sampled positions in suffix array order.
Several implementations of FM-indexes also need a bitmap B, where
B[i] = 1 denotes that SA[i]is saved inSSA (0ijTj 1). If
B[i] = 1 ,SA[i]can be found in SSA[rank 1(B;i)]. We refer to the
surveys (Ferragina et al., 2009; Navarro and Mäkinen, 2007) for more
information about the implementations of FM-indexes.
3 Methods
3.1 Analysis of Existing Locating Algorithm
The first performance bottleneck is that existing locating algorithm needs
to perform a large number of LF operations. As mentioned in Section
2.2, FM-indexes sample the positions in SAto reduce the space usage. A
popular sampling strategy is to sample every SA[i]ifSA[i]modD= 0,
whereDis the regular sampling distance. Here we refer to this strategy
asvalue sampling strategy. This strategy guarantees that any occurrence
position of a pattern can be obtained in at most D 1steps of LF operation.
Thus, to locate a pattern with occoccurrence positions, the number of LF
operations is (D 1)occin worst case. Unfortunately, short and frequent
patterns with large value of occare widely used in practice, which results
in massive LF operations. Apart from value sampling strategy, many FM-
index-based bioinformatics algorithms adopt another sampling strategy,
called subscript sampling strategy. This strategy samples every SA[i]ifi
modD= 0, whereDis the regular sampling distance. Compared with
value sampling strategy, subscript sampling strategy leads to even worse
locating performance. The reason is that for any position, value sampling
strategy can guarantee to obtain it in at most D 1steps of LF operation,
while subscript sampling strategy cannot.
ii
“main” — 2017/4/19 — 0:38 — page 4 — #4ii
i
ii
i4 Cheng et al.
Besides, the poor data locality is another bottleneck of existing locating
algorithm. To locate a position SA[i0], existing locating algorithm needs
to performmsteps of LF operation until a sampled position SA[im]is
reached. Essentially, this procedure scans text Tbackward from SA[i0]to
SA[im], so thatSA[im] =SA[i0] m. Int-th step (t= 1;2;:::;m ),
the aim of existing locating algorithm is to calculate itsuch thatSA[it] =
SA[i0] t=SA[it 1] 1. As shown in Algorithm 2, itis obtained
by calculating LF(it 1), which needs to access BWT (T)[it 1]. In
addition, if a FM-index is sampled by value sampling strategy, in t-th step,
we need to access B[it 1]to check ifSA[it 1]is a sampled position.
Thus, the memory access addresses to BWT (T)andBint-th step are
determined by it 1. Similarly, in (t+ 1) -th step, the memory access
addresses are determined by it. If the suffix T[SA[it 1];jTj 1]is not
lexicographically similar to the suffix T[SA[it];jTj 1],it 1andit
would be very different. In this case, the memory accesses of these two
steps are non-contiguous. For example, consider a non-sampled position
SA[1003] of patternP=“acgagt” and BWT (T)[1003] =c. To obtain
SA[1003] , existing locating algorithm checks B[1003] and calculates
LF(1003) in first step. Since BWT (T)[1003] =T[SA[1003] 1] =c
andSA[LF(1003)] = SA[1003] 1,SA[LF(1003)] is actually
an occurrence position of cP=“cacgagt”. Obviously, “c acgagt” is
significantly lexicographically larger than “acgagt”, so that LF(1003) is
very different to 1003. In second step, existing locating algorithm of FM-
indexes needs to check B[LF(1003)] and calculate LF(LF(1003)) .
Therefore, the memory access addresses in first step and second step to
BWT (T)andBare non-contiguous.
3.2 Our Proposed Algorithm: FMtree
The key idea of FMtree is to organize the search space of the locating
operation into a conceptual quadtree, so that multiple locations can be
located simultaneously by traversing the quadtree. This idea is based
on the observation that, different occurrence positions of a pattern P
may be obtained by performing similar LF operations. For example,
consider two non-sampled positions SA[2136] andSA[2137] of pattern
P=“acgagt”, and BWT (T)[2136] =BWT (T)[2137] =t. In first
step, existing locating algorithm calculates LF(2136) forSA[2136] ,
and calculates LF(2137) forSA[2137] . SinceBWT (T)[2136] =
BWT (T)[2137] =t, bothSA[LF(2136)] andSA[LF(2137)] are
the occurrence positions of tP=“tacgagt”. As a result, LF(2136)
andLF(2137) are very similar. In second step, for SA[2136] and
SA[2137] , existing locating algorithm calculates LF(LF(2136)) and
LF(LF(2137)) , respectively. Thus, if we locate SA[2136] and
SA[2137] jointly, namely, we calculate LF(2136) andLF(2137)
together in first step, and calculate LF(LF(2136)) andLF(LF(2137))
together in second step, the memory accesses in each step to BWT (T)
andBwould be contiguous.
More precisely, given two positions SA[i]andSA[j]with
BWT (T)[i] =BWT (T)[j] =s, if there does not exist sin
BWT (T)[i+1;j 1], it is obvious that LF(i) =LF(j) 1. Thus, for
the positions in SA[sp;ep ]with same character sinBWT (T), the results
of one step of LF operation belong to the range [LF(sps);LF(eps)],
whereSA[sps]andSA[eps]are the first position and last position in
SA[sp;ep ]with character sinBWT (T). IfSA[sp;ep ]consists of
all occurrence positions of pattern P,SA[LF(sps);LF(eps)]actually
includes all occurrence positions of sP. Therefore, [LF(sps);LF(eps)]
= backward_search (sP) = [C[s] +rank s(BWT (T);sp);C[s] +
rank s(BWT (T);ep+1) 1], as shown in Algorithm 1. In other words,
to obtain all positions in SA[sp;ep ], large numbers of LF operations for
all positions can be reduced to a few rank operations only for spandep.
Formally, this is based on the following theorem:Theorem 1. Given a text Tover alphabet and its FM-index which
is sampled by value sampling strategy with sampling distance D. Let
FM(P;T)be the sampled position set including all sampled occurrence
positions ofPin FM-index, and L(P;T)be the position set including
all occurrence positions of PinT. ThenL(P;T)can be calculated as
follows:
L(P;T) =D 1[
i=0fxjx=y+i;y2FM(iP;T)g (2)
where * is a wildcard of .
Proof. Generally, all positions in Tare classified into Dsets:
g0;g1;:::;g (D 1). Each setgiconsists of every position SA[j]ofT,
whereSA[j]modD=i. LetLgi(P;T)be the position set including
all occurrence positions of Pingi. Note thatFM(P;T) =Lg0(P;T).
Obviously,Lgi(P;T)can be obtained as:
Lgi(P;T) =fxjx=y+i;y2FM(iP;T)g (3)
SinceL(P;T)consists of all occurrence positions of PinDsets(g0;
g1;:::;g (D 1)),L(P;T)is:
L(P;T) =D 1[
i=0Lgi(P;T) (4)
Basic algorithm of FMtree. According to Theorem 1, we propose the core
algorithm of FMtree. To utilize FMtree, FM-indexes must be sampled by
value sampling strategy with regular sampling distance D. When locating
a patternPvia FM-indexes, there are total Dsteps in FMtree. Specifically,
ini-th step (i= 0;1;:::;D 1), this algorithm consists of the following
three stages:
FMtree first searches iPvia FM-indexes to obtain their
corresponding SA ranges. In total, there are | |iSA ranges ini-
th step, sinceiPrepresents | |idifferent strings. For each string
S[0;jSj 1]ini-th step (i= 1;:::;D 1), itsSArange [sp;ep ]
is updated from the SArange [sp1;ep1]of stringS[1;jSj 1]in
(i 1)-th step by two rank operations to BWT (T). More precisely,
sp=C[S[0]] +rank S[0](BWT (T);sp1), andep=C[S[0]] +
rank S[0](BWT (T);ep1+ 1) 1, as shown in Algorithm 1. Note
that theSArange ofPin 0-th step has been calculated in advance
exploiting the backward search algorithm.
Then FMtree needs to retrieve the sampled positions in these SA
ranges. As shown in Section 2.2, FM-index saves all sampled positions
ofSAinSSA in suffix array order, and utilizes bitmap Bsuch that
B[i] = 1 denotes that SA[i]is saved inSSA . Therefore, given
aSArange [sp;ep ], all sampled positions in SA[sp;ep ]are saved
consecutively in SSA[ssp;sep ], wheressp=rank 1(B;sp )and
sep=rank 1(B;ep ).
Finally, once all sampled positions of iPhave been obtained, FMtree
addsito these positions to obtain the occurrence positions of P.
With the above three stages, in i-th step, FMtree is able to obtain all
positions inLgi(P;T). Thus,L(P;T)can be obtained in total Dsteps
of FMtree. For genomic data with alphabet size jj= 4, the search space
of FMtree is actually a quadtree of height D, as shown in Fig. 2. Indeed,
thei-th step of FMtree corresponds to the i-th layer of this quadtree. We
observe that FMtree locates all occurrence positions of Pblock-by-block,
while existing locating algorithm has to locate these positions one-by-one.
ii
“main” — 2017/4/19 — 0:38 — page 5 — #5ii
i
ii
iFMtree 5
aP P
cP gP tP
aaP caP gaP taP a g t
a c g tc
acP ccP gcP tcP a c g t
agP cgP ggPtgP a c g t
atP ctP gtP ttP a c g t
Fig. 2. An example illustrates the search space of the basic FMtree when the sampling
distance D= 3 and the pattern is P.
Further optimizations of FMtree. The key problem of FMtree is that with
the increasing value of sampling distance D, the number of rank operations
increases exponentially. In i-th step, the number of rank operations to
BWT (T)is2jji= 24i. In addition, an equal number of rank
operations to Bare also required to determine the SSA ranges. As such
there are many rank operations when Dis large.
To solve this problem, we first propose an optimization to FMtree,
called early leaf nodes calculation . For the basic FMtree, we observe that
the cost of calculating its leaf nodes (i.e., the (D 1)-th step of FMtree)
dominates its overall locating time. The reason is that in (D 1)-th step, the
number of rank operations to both BWT (T)andBis 24D 1, which
is larger than the total number of rank operations in the rest D 1steps
of FMtree. Early leaf nodes calculation is proposed to avoid the expensive
(D 1)-th step of FMtree. It is based on the following theorem:
Theorem 2. Given a textT[0;jTj 1]and a pattern P[0;jPj 1], the
position setL(P;T), which includes all occurrence positions of PinT,
can be calculated as follows:
L(P;T) =D[
i=1fxjx=y i;y2FM(P[i;jPj 1];T);yi;
P[0;i 1] =T[y i;y 1]g(5)
Proof. Since FM(P[i;jPj 1];T) =Lg0(P[i;jPj 1];T), for each
positionSA[j]inFM(P[i;jPj 1];T),(SA[j] i)modD=D i.
ThenLg(D i)(P[0;jPj 1];T)can be obtained as:
Lg(D i)(P;T) =fxjx=y i;y2FM(P[i;jPj 1];T);
yi;P[0;i 1] =T[y i;y 1]g(6)
Therefore,L(P;T)is:
L(P;T) =D[
i=1Lg(D i)(P;T) (7)
Actually, the aim of the (D 1)-th step in basic FMtree is to
obtainLg(D 1)(P;T). According to Theorem 2, this position set can
be obtained using early leaf nodes calculation in following three stages:
Early leaf nodes calculation first searches P[1;jPj 1]via FM-
indexes to obtain its corresponding SArange [sp1;ep1].
Then for every SA[j](sp1jep1) ifB[j] = 1 andT[SA[j]
1] =P[0], early leaf nodes calculation adds it to a position set R. Note
that sinceT[SA[j] 1] =BWT (T)[j], early leaf nodes calculation
actually checks all elements in B[sp1;ep1]andBWT (T)[sp1;ep1].
Thus, the memory accesses in this stage are highly local. In contrast,
the(D 1)-th step of basic FMtree results in many random memory
accesses. A more practical and efficient implementation about this
stage can be found in Supplementary Section S1.Algorithm 3: FMtree(P;sp;ep;sp 1;ep1;D)
Input: the patternP[0;jPj 1]; theSArange [sp;ep ]ofP[0;jPj 1];
theSArange [sp1;ep1]ofP[1;jPj 1]; the sampling distance D
Output: the position set Rconsists of all positions in SA[sp;ep ]
1:total _num =ep sp+ 1,num = 0
2: perform early leaf nodes calculation exploiting [sp1;ep1]andP[0];
add the obtained mpositions toR;num =num +m
3:tree _height =D 1
4:node:sp =sp;node:ep =ep;node:layer = 0
5:Queue:EnQueue (node)
6:whileQueue is not empty &&num<total _num do
7:Queue:DeQueue (node)
8:sp=node:sp;ep =node:ep;layer =node:layer
9:ifep sp+ 1<threshold then
10: calculate positions in SA[sp;ep ]one-by-one in at most
tree _height layer 1steps of LF operation; add the
obtainedmpositions toR;num =num +m.
11: else
12:ssp=rank 1(B;sp );sep =rank 1(B;ep )
13:num =num +sep ssp+ 1
14: fork=ssptosepdo
15: Add SSA[k] +layer toR
16: end for
17: iflayer + 1<tree _height then
18:spchild [0] =C[a] +rank a(BWT (T);sp)
19:spchild [1] =C[c] +rank c(BWT (T);sp)
20:spchild [2] =C[g] +rank g(BWT (T);sp)
21:spchild [3] =C[t] +rank t(BWT (T);sp)
22:epchild [0] =C[a] +rank a(BWT (T);ep+ 1) 1
23:epchild [1] =C[c] +rank c(BWT (T);ep+ 1) 1
24:epchild [2] =C[g] +rank g(BWT (T);ep+ 1) 1
25:epchild [3] =C[t] +rank t(BWT (T);ep+ 1) 1
26: fort= 0 to3do
27:node:sp =spchild [t];node:ep =epchild [t]
28:node:layer =layer + 1
29:Queue:EnQueue (node)
30: end for
31: end if
32: end if
33:end while
34:returnR
Finally, early leaf nodes calculation subtracts 1 from the positions in
Rwhich has been obtained in second stage.
Apart from early leaf nodes calculation, two simple branch-cut
strategies are proposed to further improve the performance of FMtree (from
0-th step to (D 2)-th step). First, for an interval SA[spi;epi]ini-th step
(0iD 2), ifepi spi+ 1is smaller than a predefined threshold,
FMtree calculates all positions in SA[spi;epi]one-by-one in at most
D i 2steps of LF operation. Thus, the number of LF operations for
SA[spi;epi]in followingD i 2steps is (epi spi+1)(D i 2)in
worst case. Note that the rank operation to BWT (T)is the dominant cost
of LF operation. When Dis large, if FMtree does not adopt this branch-cut
strategy, the number of rank operations to BWT (T)in followingD i 2
steps is 2(41+42+:::+4D i 2) = 8(4D i 2 1)=3, which is
much larger than (epi spi+1)(D i 2). Similarly, the number of
rank operations to Bcan also be reduced. Second, when locating a pattern
withoccoccurrence positions, FMtree terminates once occoccurrence
positions have been obtained.
ii
“main” — 2017/4/19 — 0:38 — page 6 — #6ii
i
ii
i6 Cheng et al.
Full algorithm of FMtree. Algorithm 3 presents the full algorithm of
FMtree. By utilizing a queue data structure Queue , the conceptual quadtree
of FMtree is traversed in breadth-first order. In fact, the height of this
quadtree isD 1instead ofD(line 3 in Algorithm 3). This is because
early leaf nodes calculation (line 2 in Algorithm 3) is used to avoid the
(D 1)-th step in basic FMtree. For any node in quadtree, the SAranges
of its four children are calculated jointly to improve the data locality (line
18-25 in Algorithm 3). A detailed analysis is presented in Supplementary
Section S2.
4 Results
In our experiments, we used the following three datasets:
Dna.200MB consists of 209.72 million characters from Pizza&Chili
corpus ( http://pizzachili.dcc.uchile.cl/ ), which is
the standard benchmark in compressed full-text indexes (Ferragina
et al., 2009).
For practical bioinformatics algorithms, their indexes must be able to
process large texts with billions of characters. Thus, the human genome
including 3.16 billion characters was used in our experiments.
Another large text is the mouse genome. It consists of 2.73 billion
characters.
These datasets were used as texts in our experiments. Like popular
bioinformatics algorithms (Li, 2013; Li et al., 2009), character nin these
three datasets was converted to one of a,c,gandtrandomly.
We first compare FMtree with two state-of-the-art methods, including
locally compressed suffix array (LCSA) (González et al., 2015; González
and Navarro, 2007) and LZ-index (Arroyuelo, 2006). LCSA is designed
specifically to accelerate the locating operation of compressed full-text
indexes, and LZ-index has been proven that it is very competitive in
locating speed (Ferragina et al ., 2009). Besides, we implemented two
FM-index-based locating algorithms: Original_v and Original_s. The only
difference between Original_v and Original_s is their sampling strategies.
Original_v locates patterns via the FM-index sampled by value sampling
strategy, while Original_s locates patterns via the FM-index sampled by
subscript sampling strategy (see Section 3.1). For FMtree, Original v and
Original s, we implemented a highly optimized FM-index for genomic
data. Its main data structures ( CandBWT (T)) are similar to those in (Li
et al., 2009), which is a well-established FM-index-based bioinformatics
algorithm (see Supplementary Section S3). Note that FMtree, Original_v
and Original_s are independent on any particular implementation of FM-
indexes, so that we did not test these locating algorithms with other
implementations of FM-indexes. Another family of compressed full-text
indexes CSAs was not tested in our experiments. The reason is that for
genomic data, previous studies (Ferragina et al ., 2009; Gog and Petri,
2014) have shown that CSAs cannot outperform FM-indexes, LCSA
and LZ-indexes for locating operation. For detailed description about the
experimental setting, please see Supplementary Section S3.
4.1 Comparison on Small Text
In the first experiment, dna.200MB was used as text to evaluate the
performance of different methods. As previous studies about the locating
operations (Ferragina et al., 2009; Gog and Petri, 2014), patterns were
generated by randomly selecting 10 short substrings of length 5 from the
text. By utilizing these short patterns with many occurrence positions,
we could focus on the performance of locating operations and ignore the
influence of counting operations. The reason is that compared with the
locating time of different methods, their counting time was negligible in
this experiment. For LZ-index, a parameter trades the locating time for
space usage. With the increasing value of , the space usage of LZ-index
decreases, but its locating time increases. We set =f1;2;3;4gin this
/s49/s48/s48 /s49/s53/s48 /s50/s48/s48 /s50/s53/s48 /s51/s48/s48 /s51/s53/s48 /s52/s48/s48 /s52/s53/s48 /s53/s48/s48 /s53/s53/s48 /s54/s48/s48 /s54/s53/s48 /s55/s48/s48 /s55/s53/s48 /s56/s48/s48 /s56/s53/s48/s49/s48/s45/s49/s49/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s32/s83/s101/s97/s114/s99/s104/s105/s110/s103/s32/s116/s105/s109/s101/s32/s40/s109/s115/s32/s112/s101/s114/s32/s112/s97/s116/s116/s101/s114/s110/s41
/s83/s112/s97/s99/s101/s32/s117/s115/s97/s103/s101/s32/s40/s77/s66/s41/s32/s76/s90/s45/s105/s110/s100/s101/s120
/s32/s76/s67/s83/s65
/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77/s116/s114/s101/s101Fig. 3. Time/space tradeoffs of different methods for dna.200MB.
experiment. For LCSA, its default parameters were directly used. For FM-
indexes used in FMtree, Original_v and Original_s, we set the sampling
distanceD=f2;3;4;5;6;7;8g, to make the space usage of FM-indexes
similar to that of LZ-index and LCSA.
As shown in Fig. 3, with similar space usage, FMtree is one or two
order of magnitude faster than other methods. Apart from FMtree, LCSA
outperforms other methods in locating speed. However, it requires much
more space than the others. We also observe that the space usage of
Original_s is slightly less than that of FMtree and Original_v. Unlike
Original_s which uses subscript sampling strategy, FMtree and Original_v
adopt value sampling strategy. Thus, FMtree and Original_v need an extra
bitmapBto mark all sampled positions in SA. In exchange, Original_s is
about 1.5 times slower than Original_v, and at least 40 times slower than
FMtree.
4.2 Comparison on Large and Practical Texts
In the second experiment, we studied the performance of different methods
on two large texts: human genome and mouse genome. To generate
patterns, we randomly extracted short substrings from both human genome
and mouse genome. For each text, we generated four datasets including
100k patterns of length 12, 16, 20 and 25, respectively. In fact, short
patterns of length 10 to 25 are widely used in existing bioinformatics
/s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48 /s52/s50/s48/s48 /s52/s52/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48 /s52/s50/s48/s48 /s52/s52/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48 /s52/s50/s48/s48 /s52/s52/s48/s48/s49/s48/s49/s49/s48/s50
/s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48 /s52/s50/s48/s48 /s52/s52/s48/s48/s49/s48/s49/s49/s48/s50/s40/s97/s41
/s32/s32/s83/s101/s97/s114/s99/s104/s105/s110/s103/s32/s116/s105/s109/s101/s32/s40 /s115/s32/s112/s101/s114/s32/s112/s97/s116/s116/s101/s114/s110/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s40/s98/s41
/s32/s32
/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s40/s99/s41
/s32/s32/s83/s101/s97/s114/s99/s104/s105/s110/s103/s32/s116/s105/s109/s101/s32/s40 /s115/s32/s112/s101/s114/s32/s112/s97/s116/s116/s101/s114/s110/s41
/s83/s112/s97/s99/s101/s32/s117/s115/s117/s97/s103/s101/s40/s77/s66/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s83/s112/s97/s99/s101/s32/s117/s115/s117/s97/s103/s101/s40/s77/s66/s41
/s32/s32
/s40/s100/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
Fig. 4. Time/space tradeoffs of different methods for human genome. For subgraph (a),
(b), (c) and (d), the length of pattern is 12, 16, 20 and 25, respectively.
ii
“main” — 2017/4/19 — 0:38 — page 7 — #7ii
i
ii
iFMtree 7
/s49/s56/s48/s48 /s50/s48/s48/s48 /s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s49/s56/s48/s48 /s50/s48/s48/s48 /s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s49/s56/s48/s48 /s50/s48/s48/s48 /s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51
/s49/s56/s48/s48 /s50/s48/s48/s48 /s50/s50/s48/s48 /s50/s52/s48/s48 /s50/s54/s48/s48 /s50/s56/s48/s48 /s51/s48/s48/s48 /s51/s50/s48/s48 /s51/s52/s48/s48 /s51/s54/s48/s48 /s51/s56/s48/s48 /s52/s48/s48/s48/s49/s48/s49/s49/s48/s50/s49/s48/s51/s40/s97/s41
/s32/s32/s83/s101/s97/s114/s99/s104/s105/s110/s103/s32/s116/s105/s109/s101/s32/s40 /s115/s32/s112/s101/s114/s32/s112/s97/s116/s116/s101/s114/s110/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s40/s98/s41
/s32/s32
/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s40/s99/s41
/s32/s32/s83/s101/s97/s114/s99/s104/s105/s110/s103/s32/s116/s105/s109/s101/s32/s40 /s115/s32/s112/s101/s114/s32/s112/s97/s116/s116/s101/s114/s110/s41
/s83/s112/s97/s99/s101/s32/s117/s115/s117/s97/s103/s101/s40/s77/s66/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
/s83/s112/s97/s99/s101/s32/s117/s115/s117/s97/s103/s101/s40/s77/s66/s41
/s32/s32
/s40/s100/s41/s32/s79/s114/s105/s103/s105/s110/s95/s118
/s32/s79/s114/s105/s103/s105/s110/s95/s115
/s32/s70/s77 /s116/s114/s101/s101
Fig. 5. Time/space tradeoffs of different methods for mouse genome. For subgraph (a), (b),
(c) and (d), the length of pattern is 12, 16, 20 and 25, respectively.
algorithms (Hach et al., 2010; Langmead and Salzberg, 2012). Note that in
this experiment, we did not tested LZ-index and LCSA due to two reasons.
First, their implementations cannot process large texts like human genome
and mouse genome. Second, the results in first experiment have shown that
they are significantly slower than FMtree.
Fig. 4 and Fig. 5 present the results of FMtree, Original_v and
Original_s with sampling distance D=f4;5;6;7;8g. With the
increasing value of D, all methods require less space. Note that in this
experiment, the space usage of FM-index is close to that of the input human
genome and mouse genome, which require 3GB RAM and 2.6GB RAM,
respectively. For human genome (see Fig. 4), FMtree is up to 62 times
faster than Original_v and Original_s. And for mouse genome (see Fig. 5),
FMtree is at most 86 times faster than other methods. We also present the
locating time of different methods in Table S1 and Table S2. With respect
to Original_v and Original_s, FMtree achieves highest gain when we focus
on the locating time instead of the overall searching time. This is because
in addition to the locating operation, the searching algorithm of FM-index
also consists of the counting operation.
5 Conclusion and future work
In this paper we propose a novel algorithm, FMtree, to accelerate the
locating operations of FM-indexes for genomic data. When searching
for a pattern via FM-indexes, FMtree builds a conceptual quadtree,
so that multiple occurrence positions of the pattern can be obtained
simultaneously by traversing this quadtree. In contrast, existing locating
algorithm has to calculate all occurrence positions one-by-one. Therefore,
FMtree reduces massive unnecessary operations and presents better data
locality. We also introduce several strategies to further speed up FMtree.
For genomic data with small alphabet size, FMtree is significantly
faster than state-of-the-art methods in our experiments. However, the
performance of FMtree decreases rapidly with the increasing size of
alphabet. In the future, it would be interesting to optimize FMtree for
the applications with large alphabet size.
Acknowledgements
This work was partially supported by the National Nature Science
Foundation of China under the grant No. 61672480 and the Program forExcellent Graduate Students in Collaborative Innovation Center of High
Performance Computing.
References
Ahmadi, A. et al . (2012) Hobbes: optimized gram-based methods for
efficient read alignment. Nucleic acids research ,40(6) , e41-e41.
Arroyuelo, D. et al. (2006) Reducing the space requirement of LZ-index.
Annual Symposium on Combinatorial Pattern Matching , 318-329.
Arroyuelo, D. et al. (2012) Stronger Lempel-Ziv based compressed text
indexing. Algorithmica ,62, 54-101.
Burrows, M. and Wheeler, D. J. (1994) A block-sorting lossless data
compression algorithm. Technical Report 124 , Digital Equipment
Corporation, California.
Cheng, H. et al. (2015) BitMapper: an efficient all-mapper based on bit-
vector computing. BMC bioinformatics ,16(1) , 192.
Deorowicz, S. and Grabowski, S. (2013) Data compression for sequencing
data. Algorithms for Molecular Biology ,8(1), 25.
Ferragina, P. and Manzini, G. (2000) Opportunistic data structures with
applications. In Foundations of Computer Science, 2000. Proceedings.
41st Annual Symposium , 390-398.
Ferragina, P. et al . (2009) Compressed text indexes: From theory to
practice. Journal of Experimental Algorithmics (JEA) ,13, 12.
Ferragina, P. et al. (2013) Distribution-aware compressed full-text indexes.
Algorithmica ,67(4) , 529-546.
Gog, S. and Petri, M. (2014) Optimized succinct data structures for massive
data. Software: Practice and Experience ,44(11) , 1287-1314.
Gog, S. et al. (2015) Improved and extended locating functionality on
compressed suffix arrays. Journal of Discrete Algorithms ,32, 53-63.
Gog, S. et al. (2017) CSA++: Fast Pattern Search for Large Alphabets.
2017 Proceedings of the Ninteenth Workshop on Algorithm Engineering
and Experiments (ALENEX) , 73-82.
González, R. and Navarro, G. (2007) Compressed text indexes with fast
locate. Annual Symposium on Combinatorial Pattern Matching , 216-
227.
González, R. et al. (2015) Locally compressed suffix arrays. Journal of
Experimental Algorithmics (JEA) ,19, 1-1.
Grabowski, S. et al . (2004) First Huffman, then Burrows-Wheeler: A
simple alphabet-independent FM-index. International Symposium on
String Processing and Information Retrieval , 210-211.
Grossi, R. and Vitter, J. S. (2005) Compressed suffix arrays and suffix trees
with applications to text indexing and string matching. SIAM Journal on
Computing ,35, 378-407.
Hach, F. et al. (2010) mrsFAST: a cache-oblivious algorithm for short-read
mapping. Nature methods ,7(8), 576-577.
Hon, W.K. et al. (2004) Practical aspects of Compressed Suffix Arrays
and FM-Index in Searching DNA Sequences. ALENEX/ANALC , 31-38.
Langmead, B. and Salzberg, S. L. (2012) Fast gapped-read alignment with
Bowtie 2. Nature methods ,9(4), 357-359.
Li, H. (2012) Exploring single-sample SNP and INDEL calling with
whole-genome de novo assembly. Bioinformatics ,28(14) , 1838-1844.
Li, H. (2013) Aligning sequence reads, clone sequences and assembly
contigs with BWA-MEM. arXiv preprint ,arXiv:1303.3997 .
Li, R. et al . (2009) SOAP2: an improved ultrafast tool for short read
alignment. Bioinformatics ,25(15) , 1966-1967.
Liu, B. et al. (2015) rHAT: fast alignment of noisy long reads with regional
hashing. Bioinformatics ,32(11) , 1625-1631.
Manber, U. and Myers, G (1993) Suffix arrays: a new method for on-line
string searches. SIAM Journal on Computing ,22(5) , 935-948.
Marco-Sola, S. et al. (2012) The GEM mapper: fast, accurate and versatile
alignment by filtration. Nature methods ,9, 1185-1188.
ii
“main” — 2017/4/19 — 0:38 — page 8 — #8ii
i
ii
i8 Cheng et al.
Navarro, G and Mäkinen, V . (2007) Compressed full-text indexes. ACM
Computing Surveys (CSUR) ,39, 2.
Sadakane, K. et al . (2003) New text indexing functionalities of the
compressed suffix arrays. Journal of Algorithms ,48, 294-313.
Schulz, M. H. et al. (2014) Fiona: a parallel and automatic strategy for
read error correction. Bioinformatics ,30(17) , i356-i363.
Simpson, J. T. and Durbin, R. (2012) Efficient de novo assembly of large
genomes using compressed data structures. Genome research ,22(3) ,549-556.
Vyverman, M. et al. (2012) Prospects and limitations of full-text index
structures in genome analysis. Nucleic acids research ,40, 6993-7015.
Xin, H. et al. (2016) Optimal seed solver: optimizing seed selection in
read mapping. Bioinformatics ,32(11) , 1632-1642.