i

i
“main” — 2017/4/19 — 0:38 — page 1 — #1

i

i

xxxxxx/xxxxxx/xxxxxx
Advance Access Publication Date: Day Month Year
Manuscript Category

Subject Section

FMtree: A fast locating algorithm of FM-indexes for
genomic data
Haoyu Cheng 1,2,3 , Ming Wu 1,2,3 and Yun Xu 1,2,3,∗
1

School of Computer Science, University of Science and Technology of China, Heifei, Anhui, 230027, China
Key Laboratory on High Performance Computing, Anhui Province
2
Collaborative Innovation Center of High Performance Computing, National University of Defense Technology, Changsha, 410073,
China
2

∗ To

whom correspondence should be addressed.

Associate Editor: XXXXXXX

arXiv:1704.04615v2 [cs.DS] 18 Apr 2017

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: xuyun@ustc.edu.cn
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 P and the text T are the sequences
over a small alphabet Σ = {a, c, g, t}. 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 (Langmead

-

and 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 and locate. Given a pattern P and a
text T , the count function is to report the number of occurrences of P
in T , while the locate function is to retrieve all occurrence positions of
P in T . 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
1

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 2 — #2

i

i

2
(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, FMindexes 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. When

Cheng et al.

locating 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 treebased 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 string S is a sequence of characters over the alphabet Σ, and the size of
Σ is |Σ|. We let |S| denote the length of S, S[i] denote the i-th character
of S (0 ≤ i ≤ |S| − 1), and S[i, j] denote the substring that starts at
S[i] and ends at S[j]. Besides, consider a character s and a string S, sS
denotes the concatenation of s and S. We also let sn denote a string of
length n such that s0 is an empty string, s1 = s and sn = 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 P in text T .
• locate(P,T): Return all occurrence positions of pattern P in text T .

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 text T and pattern P belong to Σ = {a, c, g, t} except T [|T | − 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
text T , its suffix array SA (Manber and Myers, 1993) saves the positions
of all suffixes of T in lexicographic order. For a pattern P , it is obvious
that the positions of all suffixes prefixed by P are saved consecutively in
an interval of SA, called SA[sp, ep]. In fact, SA[sp, ep] consists of all
occurrence positions of P in T .
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)
of T , called BW T (T ), which permutes the characters of T reversibly.
Conceptually, BW T (T ) can be constructed in the following two steps:
• Building a conceptual matrix M (T ) including all cyclic rotations of

T in lexicographic order. Each row in M (T ) is a cyclic rotation of T .
• Let BW T (T ) be the last column of M (T ).

Matrix M(T)

F

$
a
a
a
a
c
c
c
g
t

acgtaacc
$acgtaac
acca$acg
cca$acgt
cgtaacca
a$acgtaa
ca$acgta
gtaacca$
taacca$a
aacca$ac

a
c
t
a
$
c
a
a
c
g

SA

BWT(T)

9
8
4
5
0
7
6
1
2
3

Fig. 1. An example of Burrows Wheeler Transform BW T (T ) and suffix array SA for
string T =“acgtaacca$".

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 3 — #3

i

FMtree

i

3

An example is presented in Fig. 1. For FM-indexes, there are two
important properties of M (T ) and BW T (T ):
• Let F be the first column of M (T ), BW T (T )[i] precedes F [i]

in text T . It is obvious that the i-th row of M (T ) corresponds to
SA[i]. Specifically, BW T (T )[i] = T [SA[i] − 1] if SA[i] 6= 0, and
BW T (T )[i] = $ if SA[i] = 0. In addition, F [i] = T [SA[i]].
• For any character s ∈ Σ, the i-th s in F and the i-th s in BW T (T )
correspond to the same s in T . For example in Fig. 1, F [5] and
BW T (T )[1] are the first c in F and BW T (T ), respectively. In fact,
both of them correspond to T [7].

Algorithm 2: locate(sp, ep)
Input: the SA range [sp, ep]
Output: the position set R consists of all positions in SA[sp, ep]
1: for i = sp to ep do
2: j = i, m = 0
3: while SA[j] is not sampled do
4:
j = LF (j), m = m + 1
5: end while
6: Add SA[j] + m to R
7: end for
8: return R

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[BW T (T )[l]] + rankBW T (T )[l] (BW T (T ), l)

(1)

where C[BW T (T )[l]] denotes the number of characters in T which are
smaller than BW T (T )[l], and rankBW T (T )[l] (BW T (T ), l) is a rank
operation that reports the number of BW T (T )[l] in BW T (T )[0, l −
1]. Essentially, LF operation scans text T backward, namely, F [l] =
T [SA[l]], F [LF [l]] = BW T (T )[l] = T [SA[l] − 1] and SA[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 P in text T , it is
obvious that all M (T )’s rows prefixed by P are saved consecutively in
an interval M (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, |P | − 1] backward in |P | steps. In i-th step (i = |P | − 1, |P | −
2, . . . , 0), this algorithm updates [sp, ep] such that M (T )[sp, ep] includes
all M (T )’s rows prefixed by P [i, |P | − 1]. Note that in first step (i =
|P |−1), [sp, ep] is the range of a single character P [|P |−1], so that it can
be directly obtained via array C (line 2 in Algorithm 1). For any character
s ∈ Σ, C[s] saves the number of characters which are lexicographically
smaller than s in T . Finally, when i = 0, the range [sp, ep] of P [0, |P |−1]
is obtained and count(P, T ) = ep − sp + 1. Since M (T )[i] corresponds
to SA[i], SA[sp, ep] actually saves all occurrence positions of P in T .
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 P in T . If the whole

Algorithm 1: backward_search(P )
Input: the pattern P [0, |P | − 1]
Output: the SA range [sp, ep] of P
1: i = |P | − 1, s = P [i], i = i − 1
2: sp = C[s], ep = C[s + 1] − 1
3: while i > 0 && sp 6 ep do
4: s = P [i], i = i − 1
5: sp = C[s] + ranks (BW T (T ), sp)
6: ep = C[s] + ranks (BW T (T ), ep + 1) − 1
7: end while
8: if sp > ep then
9: return "not found"
10: else
11: return [sp, ep]
12: end if

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] (sp ≤ i ≤ ep), this algorithm first scans text T
backward by performing LF operation m times 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 BW T (T ) into small blocks. For the beginning line of
each block and each character s ∈ Σ, the rank value is precomputed and
saved. When calculating ranks (BW T (T ), l), these implementations
first retrieve the precomputed rank value of the block which includes
BW T (T )[l], and then add the number of rest s in this block before lth line. Besides, the locating algorithm of FM-indexes requires a sampled
suffix array SSA, 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 in SSA (0 ≤ i ≤ |T | − 1). If
B[i] = 1, SA[i] can be found in SSA[rank1 (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 SA to reduce the space usage. A
popular sampling strategy is to sample every SA[i] if SA[i] mod D = 0,
where D is the regular sampling distance. Here we refer to this strategy
as value sampling strategy. This strategy guarantees that any occurrence
position of a pattern can be obtained in at most D−1 steps of LF operation.
Thus, to locate a pattern with occ occurrence positions, the number of LF
operations is (D−1)×occ in worst case. Unfortunately, short and frequent
patterns with large value of occ are widely used in practice, which results
in massive LF operations. Apart from value sampling strategy, many FMindex-based bioinformatics algorithms adopt another sampling strategy,
called subscript sampling strategy. This strategy samples every SA[i] if i
mod D = 0, where D is 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 − 1 steps of LF operation,
while subscript sampling strategy cannot.

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 4 — #4

i

i

4

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 perform m steps of LF operation until a sampled position SA[im ] is
reached. Essentially, this procedure scans text T backward from SA[i0 ] to
SA[im ], so that SA[im ] = SA[i0 ] − m. In t-th step (t = 1, 2, . . . , m),
the aim of existing locating algorithm is to calculate it such that SA[it ] =
SA[i0 ] − t = SA[it−1 ] − 1. As shown in Algorithm 2, it is obtained
by calculating LF (it−1 ), which needs to access BW T (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 if SA[it−1 ] is a sampled position.
Thus, the memory access addresses to BW T (T ) and B in t-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 ], |T | − 1] is not
lexicographically similar to the suffix T [SA[it ], |T | − 1], it−1 and it
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 pattern P = “acgagt” and BW T (T )[1003] = c. To obtain
SA[1003], existing locating algorithm checks B[1003] and calculates
LF (1003) in first step. Since BW T (T )[1003] = T [SA[1003]−1] = c
and SA[LF (1003)] = SA[1003] − 1, SA[LF (1003)] is actually
an occurrence position of cP = “cacgagt”. Obviously, “cacgagt” is
significantly lexicographically larger than “acgagt”, so that LF (1003) is
very different to 1003. In second step, existing locating algorithm of FMindexes needs to check B[LF (1003)] and calculate LF (LF (1003)).
Therefore, the memory access addresses in first step and second step to
BW T (T ) and B are non-contiguous.

Theorem 1. Given a text T over alphabet Σ and its FM-index which
is sampled by value sampling strategy with sampling distance D. Let
F M (P, T ) be the sampled position set including all sampled occurrence
positions of P in FM-index, and L(P, T ) be the position set including
all occurrence positions of P in T . Then L(P, T ) can be calculated as
follows:

3.2 Our Proposed Algorithm: FMtree

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 pattern P via FM-indexes, there are total D steps in FMtree. Specifically,
in i-th step (i = 0, 1, . . . , D − 1), this algorithm consists of the following
three stages:

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] and SA[2137] of pattern
P = “acgagt”, and BW T (T )[2136] = BW T (T )[2137] = t. In first
step, existing locating algorithm calculates LF (2136) for SA[2136],
and calculates LF (2137) for SA[2137]. Since BW T (T )[2136] =
BW T (T )[2137] = t, both SA[LF (2136)] and SA[LF (2137)] are
the occurrence positions of tP = “tacgagt”. As a result, LF (2136)
and LF (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) and LF (2137)
together in first step, and calculate LF (LF (2136)) and LF (LF (2137))
together in second step, the memory accesses in each step to BW T (T )
and B would be contiguous.
More precisely, given two positions SA[i] and SA[j] with
BW T (T )[i] = BW T (T )[j] = s, if there does not exist s in
BW T (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 s in BW T (T ), the results
of one step of LF operation belong to the range [LF (sps ), LF (eps )],
where SA[sps ] and SA[eps ] are the first position and last position in
SA[sp, ep] with character s in BW T (T ). If SA[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] + ranks (BW T (T ), sp), C[s] +
ranks (BW T (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 sp and ep.
Formally, this is based on the following theorem:

L(P, T ) =

D−1
[

{x|x = y + i, y ∈ F M (∗i P, T )}

(2)

i=0

where * is a wildcard of Σ.
Proof. Generally, all positions in T are classified into D sets:
g0 , g1 , . . . , g(D−1) . Each set gi consists of every position SA[j] of T ,
where SA[j] mod D = i. Let Lgi (P, T ) be the position set including
all occurrence positions of P in gi . Note that F M (P, T ) = Lg0 (P, T ).
Obviously, Lgi (P, T ) can be obtained as:
Lgi (P, T ) = {x|x = y + i, y ∈ F M (∗i P, T )}

(3)

Since L(P, T ) consists of all occurrence positions of P in D sets (g0 ,
g1 , . . . , g(D−1) ), L(P, T ) is:

L(P, T ) =

D−1
[

Lgi (P, T )

(4)

i=0

• FMtree first searches ∗i P

via FM-indexes to obtain their
corresponding SA ranges. In total, there are |Σ|i SA ranges in ith step, since ∗i P represents |Σ|i different strings. For each string
S[0, |S| − 1] in i-th step (i = 1, ..., D − 1), its SA range [sp, ep]
is updated from the SA range [sp1 , ep1 ] of string S[1, |S| − 1] in
(i − 1)-th step by two rank operations to BW T (T ). More precisely,
sp = C[S[0]] + rankS[0] (BW T (T ), sp1 ), and ep = C[S[0]] +
rankS[0] (BW T (T ), ep1 + 1) − 1, as shown in Algorithm 1. Note
that the SA range of P in 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
of SA in SSA in suffix array order, and utilizes bitmap B such that
B[i] = 1 denotes that SA[i] is saved in SSA. Therefore, given
a SA range [sp, ep], all sampled positions in SA[sp, ep] are saved
consecutively in SSA[ssp, sep], where ssp = rank1 (B, sp) and
sep = rank1 (B, ep).
• Finally, once all sampled positions of ∗i P have been obtained, FMtree
adds i to 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 in Lgi (P, T ). Thus, L(P, T ) can be obtained in total D steps
of FMtree. For genomic data with alphabet size |Σ| = 4, the search space
of FMtree is actually a quadtree of height D, as shown in Fig. 2. Indeed,
the i-th step of FMtree corresponds to the i-th layer of this quadtree. We
observe that FMtree locates all occurrence positions of P block-by-block,
while existing locating algorithm has to locate these positions one-by-one.

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 5 — #5

i

FMtree

5

Algorithm 3: FMtree(P, sp, ep, sp1 , ep1 , D)

P
c

a
aP
a
aaP

c
caP

i

g

t

cP
g
gaP

t

a
taP

acP

c
ccP

gP
g
gcP

t

a
tcP

agP

c
cgP

tP

g

t

ggP

a
tgP

atP

c
ctP

g
gtP

t
ttP

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
BW T (T ) is 2 × |Σ|i = 2 × 4i . In addition, an equal number of rank
operations to B are also required to determine the SSA ranges. As such
there are many rank operations when D is 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 BW T (T ) and B is 2 ×4D−1 , which
is larger than the total number of rank operations in the rest D − 1 steps
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 text T [0, |T | − 1] and a pattern P [0, |P | − 1], the
position set L(P, T ), which includes all occurrence positions of P in T ,
can be calculated as follows:

L(P, T ) =

D
[

{x|x = y − i, y ∈ F M (P [i, |P | − 1], T ), y ≥ i,
(5)

i=1

P [0, i − 1] = T [y − i, y − 1]}
Proof. Since F M (P [i, |P |−1], T ) = Lg0 (P [i, |P |−1], T ), for each
position SA[j] in F M (P [i, |P | − 1], T ), (SA[j] − i) mod D = D − i.
Then Lg(D−i) (P [0, |P | − 1], T ) can be obtained as:
Lg(D−i) (P, T ) = {x|x = y − i, y ∈ F M (P [i, |P | − 1], T ),

(6)

y ≥ i, P [0, i − 1] = T [y − i, y − 1]}

Input: the pattern P [0, |P | − 1]; the SA range [sp, ep] of P [0, |P | − 1];
the SA range [sp1 , ep1 ] of P [1, |P | − 1]; the sampling distance D
Output: the position set R consists of all positions in SA[sp, ep]
1: total_num = ep − sp + 1, num = 0
2: perform early leaf nodes calculation exploiting [sp1 , ep1 ] and P [0];
add the obtained m positions to R; num = num + m
3: tree_height = D − 1
4: node.sp = sp, node.ep = ep, node.layer = 0
5: Queue.EnQueue(node)
6: while Queue is not empty && num < total_num do
7: Queue.DeQueue(node)
8: sp = node.sp, ep = node.ep, layer = node.layer
9: if ep − sp + 1 < threshold then
10:
calculate positions in SA[sp, ep] one-by-one in at most
tree_height − layer − 1 steps of LF operation; add the
obtained m positions to R; num = num + m.
11: else
12:
ssp = rank1 (B, sp), sep = rank1 (B, ep)
13:
num = num + sep − ssp + 1
14:
for k = ssp to sep do
15:
Add SSA[k] + layer to R
16:
end for
17:
if layer + 1 < tree_height then
18:
spchild[0] = C[a] + ranka (BW T (T ), sp)
19:
spchild[1] = C[c] + rankc (BW T (T ), sp)
20:
spchild[2] = C[g] + rankg (BW T (T ), sp)
21:
spchild[3] = C[t] + rankt (BW T (T ), sp)
22:
epchild[0] = C[a] + ranka (BW T (T ), ep + 1) − 1
23:
epchild[1] = C[c] + rankc (BW T (T ), ep + 1) − 1
24:
epchild[2] = C[g] + rankg (BW T (T ), ep + 1) − 1
25:
epchild[3] = C[t] + rankt (BW T (T ), ep + 1) − 1
26:
for t = 0 to 3 do
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: return R

Therefore, L(P, T ) is:

L(P, T ) =

D
[

Lg(D−i) (P, T )

(7)

i=1

Actually, the aim of the (D − 1)-th step in basic FMtree is to
obtain Lg(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, |P | − 1] via FM-

indexes to obtain its corresponding SA range [sp1 , ep1 ].
• Then for every SA[j] (sp1 ≤ j ≤ ep1 ) if B[j] = 1 and T [SA[j] −

1] = P [0], early leaf nodes calculation adds it to a position set R. Note
that since T [SA[j] − 1] = BW T (T )[j], early leaf nodes calculation
actually checks all elements in B[sp1 , ep1 ] and BW T (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.

• Finally, early leaf nodes calculation subtracts 1 from the positions in

R which 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 ] in i-th step
(0 ≤ i ≤ D − 2), if epi − spi + 1 is smaller than a predefined threshold,
FMtree calculates all positions in SA[spi , epi ] one-by-one in at most
D − i − 2 steps of LF operation. Thus, the number of LF operations for
SA[spi , epi ] in following D−i−2 steps is (epi −spi +1)×(D−i−2) in
worst case. Note that the rank operation to BW T (T ) is the dominant cost
of LF operation. When D is large, if FMtree does not adopt this branch-cut
strategy, the number of rank operations to BW T (T ) in following D−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 B can also be reduced. Second, when locating a pattern
with occ occurrence positions, FMtree terminates once occ occurrence
positions have been obtained.

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 6 — #6

i

i

6

Cheng et al.
3

10

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 is D − 1 instead of D (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 SA ranges
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.

LZ-index

Searching time (ms per pattern)

LCSA

4 Results

Origin_v
2

Origin_s

10

FMtree

1

10

0

10

-1

In our experiments, we used the following three datasets:

10

100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850

Space usage (MB)

• Dna.200MB consists of 209.72 million characters from Pizza&Chili

These datasets were used as texts in our experiments. Like popular
bioinformatics algorithms (Li, 2013; Li et al., 2009), character n in these
three datasets was converted to one of a, c, g and t randomly.
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 (C and BW T (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 FMindexes, 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.

Fig. 3. Time/space tradeoffs of different methods for dna.200MB.

experiment. For LCSA, its default parameters were directly used. For FMindexes used in FMtree, Original_v and Original_s, we set the sampling
distance D = {2, 3, 4, 5, 6, 7, 8}, 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
bitmap B to 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

10

Searching time ( s per pattern)

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.

10

3

10

2

10

1

Origin_v

Origin_s

Origin_s

FMtree

FMtree

10

2

10

1

2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400

4.1 Comparison on Small Text

3

Origin_v

2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400

(a)

(b)

Origin_v

FMtree

Searching time ( s per pattern)

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  = {1, 2, 3, 4} in this

Origin_v

Origin_s

10

2

10

1

2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400

10

2

10

1

Origin_s
FMtree

2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400

(c)

(d)

Space usuage(MB)

Space usuage(MB)

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.

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 7 — #7

i

FMtree

7
Origin_v

Origin_v

Origin_s

Searching time ( s per pattern)

FMtree

10

3

10

2

10

1

10

3

10

2

10

1

10

3

1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000

Origin_s
FMtree

Searching time ( s per pattern)

3

10

2

1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000

10

1

(b)

Origin_v

Origin_v

Origin_s

Origin_s

FMtree

FMtree

1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000

10

2

10

1

Excellent Graduate Students in Collaborative Innovation Center of High
Performance Computing.

References

(a)

10

i

1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000

(c)

(d)

Space usuage(MB)

Space usuage(MB)

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 = {4, 5, 6, 7, 8}. 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 for

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 bitvector 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, 216227.
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.

i

i
i

i

i

i
“main” — 2017/4/19 — 0:38 — page 8 — #8

i

8
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),

i

Cheng et al.

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.

i

i
i

i

