Indexing Variation Graphs∗
Jouni Sirén†

arXiv:1604.06605v4 [cs.DS] 12 Jan 2017

Abstract
Variation graphs, which represent genetic variation within
a population, are replacing sequences as reference genomes.
Path indexes are one of the most important tools for working
with variation graphs. They generalize text indexes to graphs,
allowing one to find the paths matching the query string. We
propose using de Bruijn graphs as path indexes, compressing
them by merging redundant subgraphs, and encoding them
with the Burrows-Wheeler transform. The resulting fast,
space-efficient, and versatile index is used in the variation
graph toolkit vg.

1

Introduction

Sequence analysis pipelines typically start with mapping
the reads from the sequenced genome to a reference
genome of the same species. As reference genomes are
usually assembled from the genomes of a small number
of individuals, they are biased towards those individuals.
This reference bias may affect the results of subsequent
analysis, especially when the sequenced individuals are
from different populations than the ones behind the
reference genome.
Variation graphs (also graph genomes, genome
graphs, graph references, or reference graphs), which
encode the genetic variation within a population as a
graph, have been proposed as a solution to the reference
bias [37, 28, 6, 8, 29]. These graphs are expected to
replace sequences as reference genomes. The shift to
graphs will likely take years, as new methods and tools
are needed to replace those based on linear references.
The variation graph toolkit vg [14] is a community
effort to develop such tools. This paper describes
GCSA2, the path index used in vg. A path index is a
generalization of text indexes for labeled graphs. Given
a query string, the index finds the paths with a label
matching the query. Indexing graphs is inherently hard,
as the number of paths increases exponentially with
path length. The design of a path index is hence a tradeoff between maximum query length, index size, query
performance, and pruning complex regions of the graph.
Mapping reads to a graph was first investigated
by Schneeberger et al. [42]. Sirén et al. developed
∗ Supported
† Wellcome

by the Wellcome Trust grant 098051.
Trust Sanger Institute, UK. jouni.siren@iki.fi.

GCSA [43], which generalized the FM-index [10] (a text
index based on the Burrows-Wheeler transform [4]) to
directed acyclic graphs. GCSA depends on pruning the
complex regions, as there are no limits on query length.
Kim et al. [19] combined GCSA with the HISAT read
aligner [20]. The result was the first practical graphbased read aligner, though it uses the graph for more
accurate mapping to a linear reference.
De Bruijn graphs can be used as k-mer indexes of
other graphs. There are already read aligners based on
them [23, 25]. The succinct de Bruijn graph of Bowe
et. al. [3] encodes the graph with a generalization of
the FM-index. Rødland [39] proposed another similar
generalization. Succinct de Bruijn graphs can simulate
order-k de Bruijn graphs for multiple values of k [2],
but they still need to store the graph explicitly for
the largest value of k. Other representations include
compacted [5] and compressed [28] de Bruijn graphs,
which represent unary paths in the graph as single
nodes or edges. Probabilistic de Bruijn graphs [38] use
Bloom filters to support faster queries, at the expense
of producing false positives. All these space-efficient
representations require several bits per k-mer for the
graph, and more for mapping back to the indexed graph.
Some path indexes store the graph as a collection of
sequences. BWBBLE [17] uses the powerset alphabet for
encoding substitutions and creates new sequences with
a sufficient amount of context for other variants. vBWT
[26] encodes variant sites explicitly in the sequence as
X(A|B|C)Y , using distinct separator symbols for each
site. Queries in both BWBBLE and vBWT are slower
than in ordinary FM-indexes, as variant sites force the
search to branch. The hypertext index [44] works with
graphs that have string labels on the nodes. The strings
are indexed using an FM-index. Partial matches in the
strings are combined into full matches with range queries
in the edge matrix. While matches crossing one edge
are easy to find, the approach becomes impractical with
matches crossing multiple edges.
There are also structures (e.g. [18, 46, 7, 32, 33])
using graphs as a space-efficient way of indexing similar
sequences. While the problem is different from indexing
the paths in a graph, the techniques used are similar.
The above methods can be classified in three
categories based on the data models they use. Many deal

with graphs arising from aligned sequences, assuming
a shared global sequence with local variation. GCSAbased methods can index directed acyclic graphs, as well
as cyclic graphs that are sufficiently similar to de Bruijn
graphs. GCSA2, the hypertext index, and indexes based
on de Bruijn graphs work with arbitrary graphs.
GCSA2 combines ideas from the original GCSA
and from succinct de Bruijn graphs. Conceptually it
uses a de Bruijn graph as a k-mer index of a variation
graph. The de Bruijn graph is pruned (compressed
structurally) by using strings shorter than k characters
as nodes, if the shorter strings identify the start nodes of
the corresponding paths uniquely. The pruned graph is
encoded with a generalization of the FM-index. GCSA2
often uses less space (e.g. less than 1 bit per k-mer)
than other de Bruijn graph-based indexes, which have
to store some information for each k-mer explicitly. The
index also includes extensions based on suffix trees.
The extensions are used for e.g. finding maximal exact
matches in the vg read aligner.
The main differences to the original GCSA are:
1. The graph encoding in GCSA2 has been optimized
for small alphabets, improving query performance
by up to an order of magnitude.
2. The construction algorithm stores the graphs on
disk, reducing the memory requirements of building
a whole-genome human index from hundreds of
gigabytes to tens of gigabytes.

and refer to substrings of length k as k-mers. Substring
S[i, j] is a proper substring of string S, if S 6= S[i, j].
We say that string S 0 is a substring of string collection
S, if it is a substring of a string S ∈ S.
Sometimes we consider infinite character sequences
S = (si )i∈Z , where set Z is a contiguous infinite subset
of Z. The notion of substring generalizes to infinite
sequences in a natural way. A substring of an infinite
sequence S is left-infinite if it extends infinitely to the left,
and right-infinite if it extends infinitely to the right. A
substring of a finite or infinite sequence S is left-maximal
if it is left-infinite or a prefix, and right-maximal if it is
right-infinite or a suffix.
We are primarily interested in sequences over the
DNA alphabet {$, A, C, G, T, N}. Characters A, C, G, and T
are called bases, while character N represents an arbitrary
or unknown base. The alphabet may contain other
technical characters in addition to the endmarker $.
Each character c of the DNA alphabet has a complement
c defined as A = T, C = G, G = C, T = A, and
c = c for other characters c. Given a DNA sequence
←
−
S, its reverse complement is the sequence S obtained
by reversing the non-technical parts of the sequence
and replacing each character with its complement. For
←−−−−−−
example, GATTACA$ = TGTAATC$.
Given a string S[0, n − 1], we define S.rank(i, c) as
the number of occurrences of character c in the prefix
S[0, i − 1]. We also define

S.select(i, c) = max{j ≤ n | S.rank(j, c) < i}
3. GCSA2 can index denser graphs, including cyclic
as the position of the occurrence of character c with
graphs, by limiting maximum query length.
rank i > 0.1 A bitvector is a binary sequence supporting
4. GCSA2 extends the functionality of the FM-index efficient rank/select queries. Wavelet trees [16] are spacewith suffix tree operations.
efficient data structures that use bitvectors to support
rank/select queries on arbitrary strings.
The first two improvements can also be used with the
Let S be a string and S 0 be a string or an infinite
original GCSA. The third point represents a different character sequence over alphabet Σ. We say that
approach to indexing graphs. It was not possible the sequences S and S 0 prefix-match, if S is a prefix of S 0 or
with the original GCSA, as the construction algorithm S 0 is a prefix of S. Set S of strings is prefix-free, if no
required a prefix-range-sorted graph. The extended two strings S, S 0 ∈ S (with S 6= S 0 ) prefix-match.
functionality depends on the new construction algorithm
and on limiting maximum query length.
2.2 Text Indexes The suffix tree [47] is the most
This paper describes the GCSA2 data structure. Its fundamental full-text index supporting substring queries.
uses in the vg toolkit will be discussed in the vg paper. It is formed by taking the suffixes of the text, storing
them in a trie, and compacting unary paths in the trie
into single edges. Although fast and versatile, suffix
2.1 Strings A string S[0, n−1] = s0 · · · sn−1 of length trees are impractical with large texts, as they require
|S| = n is a sequence of characters over an alphabet much more space than the text itself.
Σ = {0, . . . , σ−1}. Text strings T [0, n−1] are terminated
Suffix arrays (SA) [27] were introduced as a spaceby an endmarker T [n − 1] = $ = 0 not found anywhere efficient alternative to the suffix tree. The suffix array
else in the text. A substring of string S is a sequence of
1 These definitions are used in the SDSL library [15]. We assume
the form S[i, j] = si · · · sj . We call substrings of the type
S[0, j] and S[i, n − 1] prefixes and suffixes, respectively, for convenience that S.select(0, c) = −1.
2

Background

Suffixes
$
A$
ATA$
ATCATA$
CATA$
CATCATA$
GCATCATA$
TA$
TCATA$

LF BWT Suffixes
1
7
4
5
8
6
0
2
3

A
T
C
C
T
G
$
A
A

$
A$
ATA$
ATCATA$
CATA$
CATCATA$
GCATCATA$
TA$
TCATA$

SA LCP
8
7
5
2
4
1
0
6
3

0
0
1
2
0
3
0
0
1

pointer. If we find a sample at SA[LFk (i)], we know that
SA[i] = (SA[LFk (i)] + k) mod n.

If we have sampled one out of d suffix array pointers at
regular intervals, finding each occurrence takes O(d) rank
queries. If we also sample one out of d0 inverse suffix
array pointers2 , we can extract an arbitrary substring
X of the text using O(|X| + d0 ) rank queries.
The longest-common-prefix array (LCP array) [27]
is an integer array LCP[0, n − 1], where each value LCP[i]
tells the length of the longest common prefix of suffixes
Figure 1: LF-mapping, BWT, suffix array, and LCP T [SA[i−1], n−1] and T [SA[i], n−1] (with LCP[0] = 0). If
array for text GCATCATA$.
we have the FM-index, the LCP array, and the topology
of the suffix tree, we get the compressed suffix tree, which
supports the full functionality of the suffix tree in a spaceof text T [0, n − 1] is an array of pointers SA[0, n − 1]
efficient manner [40].
to the suffixes of the text in lexicographic order. We
find the occurrences of pattern string X in the text in
2.3 Graphs A graph G = (V, E) consists of a set of
O(|X| log n) time by using binary search in the suffix
nodes V = {0, . . . , |V |−1} and a set of edges E ⊆ V ×V .
array. The suffix array requires n log n bits of space
We say that (u, v) ∈ E is an edge from node u to node
in addition to the text, while its functionality is more
v, and assume that the edges are directed : (u, v) 6= (v, u)
limited than that of the suffix tree. See Figure 1 for an
for u 6= v. The indegree G.in(v) of node v is the number
example of the suffix array and related structures.
of incoming edges to v, while the outdegree G.out(v) is
The Burrows-Wheeler transform (BWT) [4] of text
the number of outgoing edges from v.
T [0, n − 1] is a permutation BWT[0, n − 1] such that
The graphs we use are labeled with alphabet Σ:
BWT[i] = T [(SA[i] − 1) mod n]. Given the lexicographic
each node v ∈ V has a label G.label(v) ∈ Σ. A path in
rank i of suffix T [SA[i], n−1], we can use LF-mapping on
a graph is a sequence of nodes P = v0 · · · v|P |−1 such
the BWT to find the lexicographic rank of the previous
that (vi , vi+1 ) ∈ E for 0 ≤ i < |P | − 1. We say that v0
suffix T [(SA[i] − 1) mod n, n − 1]. Let
is the start node and v|P |−1 is the end node of the path.
The label of path P is the concatenation of node labels
LF(i) = C[BWT[i]] + BWT.rank(i, BWT[i]),
G.label(P ) = G.label(v0 ) · · · G.label(v|P |−1 ).
If the graph has nodes with indegree or outdegree
where C[c] is the number of occurrences of characters
c0 < c in the BWT. Then SA[LF(i)] = (SA[i] − 1) mod n. 0, we add a source node s and a sink node t to it.
To distinguish these technical nodes from the actual
We generalize the definition to any character c ∈ Σ:
nodes, we label them with characters G.label(s) = #
LF(i, c) = C[c] + BWT.rank(i, c).
and G.label(t) = $, which are not used anywhere else in
the graph. We add an edge (s, v) to all nodes v ∈ V \ {s}
0
Let X be a string. If there are i suffixes S of text T with no incoming edges, and an edge (v, t) from all nodes
0
such that S < X in lexicographic order, then there are v ∈ V \ {t} with no outgoing edges. We also add edge
LF(i, c) suffixes S 0 such that S 0 < cX.
(t, s) to guarantee that G.in(v) ≥ 1 and G.out(v) ≥ 1 for
We can use the BWT as a space-efficient text index. all nodes v ∈ V . However, this edge is not considered a
The FM-index [10] combines a representation of the real edge, and no path can cross it.
BWT supporting rank/select queries, the C array, and
We also consider infinite paths P = (vi )i∈Z , genera set of sampled pointers from the suffix array. It uses alizing the definition in a similar way as we did with
backward searching to find the lexicographic range of infinite character sequences in Section 2.1. We say that
suffixes matching pattern X (having X as a prefix). If path P is left-maximal if it starts at the source node or
the lexicographic range matching suffix X[i + 1, |X| − 1] extends infinitely to the left; right-maximal if it ends
of the pattern is SA[sp, ep], the range matching suffix at the sink node or extends infinitely to the right; and
X[i, |X| − 1] is SA[LF(sp, X[i]), LF(ep + 1, X[i]) − 1]. maximal if it is both left-maximal and right-maximal.
Matching the entire pattern takes O(|X|) rank queries.
We use the sampled suffix array pointers for finding
2 The suffix array is a permutation of {0, . . . , n − 1}, and the
the text positions containing the occurrences. If SA[i]
is not sampled, we iterate LF(i) until we find a sampled inverse suffix array is the inverse permutation.

Given a graph G = (V, E), we may want to reason
about the predecessors of a node with the given label.
Let v ∈ V be a node and c ∈ Σ be a character. We write
G.pred(v, c) to denote the set of nodes u ∈ V such that
G.label(u) = c and there is an edge (u, v) ∈ E.
We work with de Bruijn graphs and their generalizations. For that purpose, we define collections of (finite
or infinite) sequences suitable for constructing order-k
de Bruijn graphs.

that includes directed acyclic graphs and de Bruijn
graphs. Before indexing, we transform the graph into
an equivalent graph, where the nodes can be sorted
unambiguously by the labels of the right-maximal paths
starting from them. The transformation increases the
size of the graph exponentially in the worst case. In
addition to sequences BWT and IN, GCSA also uses an
outdegree bitvector OUT, which is encoded in the same
way as IN. LF-mapping uses select queries on bitvector
IN to map nodes to BWT ranges, ordinary LF-mapping
Definition 2.1. (k-collection) Let S be a collection with BWT to map incoming edges to the corresponding
of character sequences over alphabet Σ, and let k > 0 be outgoing edges, and rank queries on bitvector OUT to
a parameter value. We say that S is a k-collection, if map the outgoing edges to the predecessor nodes.
each sequence S ∈ S (a) is left-infinite or begins with #k ;
(b) is right-infinite or ends with $k ; and (c) contains no 3 Path Indexes
other occurrences of characters # and $.
A path index is a generalization of text indexes for labeled
Definition 2.2. (de Bruijn graph) Let k > 0, and graphs. Given a path index for input graph G = (V, E),
let S be a k-collection. The order-k de Bruijn graph of we use the index to find the start nodes v0 ∈ V of the
paths P = v0 · · · v|X|−1 matching pattern X (paths P
S is a graph G = (V, E) such that
with G.label(P ) = X).
• each node vX ∈ V represents a distinct k-mer X
The proofs of the lemmas can be found in Apoccurring in S, with G.label(vX ) = X[0];
pendix A.
• each node vX ∈ V has a key G.key(vX ) = X; and

3.1 De Bruijn Graphs as Path Indexes The
• each edge (vX , vY ) ∈ E represents a (k + 1)-mer k-mer index is the simplest path index. It consists
X[0]Y = XY [k − 1] occurring in S.
of a set of key-value pairs (X, VX ), where X ∈ Σk is a
k-mer and VX ⊆ V is the set of the start nodes of the
We use #k as the source node s and $k as the sink node
paths matching the k-mer. If we store the pairs in a hash
t, adding the edge (t, s) in the usual way.
table, we can quickly search for patterns of length k. If
De Bruijn graphs have several properties that make we use binary search in a sorted list of pairs, queries
them useful for indexing purposes. Node keys prefix- become slower, but we gain the ability to search for
match the labels of all paths starting from the node. patterns shorter than k characters. The main drawback
This makes it possible to sort the nodes unambiguously of these basic k-mer indexes is their size, as they store
by path labels. Every substring of the k-collection is the key-value pairs explicitly.
We can represent k-mer indexes as de Bruijn graphs.
the label of a path in the de Bruijn graph, and every
path label of length at most k + 1 is a substring of the For that purpose, we define the de Bruijn graph of graph
collection. In Section 3, we develop an index structure G = (V, E) by using the collection S of the labels of the
maximal paths in the graph. If sequence S ∈ S is the
based on a generalization of de Bruijn graphs.
label of path P = (vi )i∈Z , we set S.node(S, i) = vi for
2.4 FM-Index for Graphs Suffix trees, suffix ar- all positions i ∈ Z. We transform S into a k-collection
rays, and the FM-index can be generalized to index by inserting characters # to the beginning of each nonmultiple texts. There are also generalizations to other left-infinite sequence when necessary, and characters $
combinatorial structures. The XBW transform [9] is an to the end of each non-right-infinite sequence. If S[i] is a
FM-index for labeled trees. The nodes of the tree are # we inserted, we set S.node(S, i) = s : j, where s ∈ V
sorted by the path labels from the node to the root. BWT is the source node and j is the distance to the nearest
stores the labels of the children of each node, while leaf non-inserted # in S.
nodes are marked with special characters. If a node has
k children, we encode that as a binary sequence 0k−1 1.
We concatenate these sequences to form the indegree
bitvector IN. The labels of the children of the ith node
are found in BWT[IN.select(i, 1) + 1, IN.select(i + 1, 1)].
The generalized compressed suffix array (GCSA)
[43] extends the XBW transform to a class of graphs

Definition 3.1. (de Bruijn graph of a graph)
Let G be a labeled graph, and let S be the k-collection of
maximal path labels in G. The order-k de Bruijn graph
of S is the order-k de Bruijn graph of graph G.
Let G0 = (V 0 , E 0 ) be a de Bruijn graph of graph
G = (V, E), and let S be the k-collection used to define

it. We attach a set of nodes of graph G to each node
v 0 ∈ V 0 as a value G0 .value(v):
{S.node(S, i) | S ∈ S, S[i, i + k − 1] = G0 .key(v)}.
Apart from some technicalities near the source/sink
nodes, G0 .value(v) is the set of the start nodes of the
paths matching pattern G0 .key(v) in graph G. The index
produces no false negatives (path labels that exist in
the input graph but not in the index). There may be
false positives (path labels that exist in the index but
not in the input graph) with patterns longer than k, but
we can avoid them by verifying the results in the input
graph. See Figure 2 for an example.
We can encode the de Bruijn graph as a GCSA
using (|V 0 | + o(|V 0 |))(log σ + 2) bits. By using a similar
sampling scheme for the values as in the FM-index, we
get a k-mer index that uses a couple of bytes per k-mer
for typical variation graphs (see Section 5). While this
is much less than with the basic k-mer indexes, it is still
too much for large variation graphs.
3.2 Path Graphs When a path of length k < k has
a unique label, its extensions become redundant nodes
in the order-k de Bruijn graph, if the path branches
after the first k 0 characters. By using shorter keys when
possible, we can build a smaller graph that is equivalent
to the de Bruijn graph as a path index.
0

Definition 3.2. (Path graph) Let k > 0, let S be the
k-collection of the labels of the maximal paths in graph
G = (V, E), and let K be a prefix-free set of substrings of
length k or less from S. Assume that each right-maximal
substring S ∈ S prefix-matches a string K ∈ K and that
$k ∈ S. The order-k path graph of graph G with key
set K is a graph G0 = (V 0 , E 0 ), where

• Pattern X ∈ Σ∗ matches node v ∈ V 0 , if there
is a path P 0 in G0 with G0 .label(P 0 ) = X. We
use G0 .find(X) to denote the set of nodes VX0 ⊆ V 0
matching the pattern.
• If VX0 ⊆ V 0 is the set of nodes matching pattern X, the set of Soccurrences for the pattern
is G0 .locate(VX0 ) = v0 ∈V 0 G0 .value(v 0 ). We use
X
G0 .locate(X) to denote G0 .locate(G0 .find(X)).
Lemma 3.1. (No false negatives) Let G0 = (V 0 , E 0 )
be a path graph of G = (V, E), and let X ∈ (Σ \ {#, $})∗
be a pattern. Set G0 .locate(X) contains the start nodes
of all paths P in graph G with G.label(P ) = X.
Lemma 3.2. (Context length) Let G0 = (V 0 , E 0 ) be
a path graph, and let X ∈ (Σ \ {#, $})∗ be a pattern.
Set G0 .f ind(X) consists of all nodes v 0 ∈ V 0 such that
X[0, m − 1] is a prefix of G0 .key(v 0 ), for a context length
m, which depends on the graph and the pattern.
Lemma 3.3. (Short keys) Let G0 = (V 0 , E 0 ) be a path
graph with |G0 .key(u0 )| ≤ |G0 .key(v 0 )| + 1 for all edges
(u0 , v 0 ) ∈ E 0 . Then (a) |G0 .pred(v 0 , c)| ≤ 1 for all nodes
v 0 ∈ V 0 and characters c ∈ Σ; and (b) key G0 .key(v 0 )
prefix-matches pattern X ∈ (Σ \ {#, $})∗ for all nodes
v 0 ∈ G0 .find(X).
A path graph may produce false positives with
patterns longer than k 0 characters, where k 0 is the length
of the shortest key. In the next section, we define a class
of path graphs that can be proven to be equivalent to
de Bruijn graphs.

3.3 Pruned de Bruijn Graphs We can compress
de Bruijn graphs structurally by merging keys sharing a
common prefix, if the corresponding values are identical.
• each node vK ∈ V 0 represents a distinct key K ∈ K, These pruned de Bruijn graphs, which arise naturally
with G0 .label(vK ) = K[0];
from GCSA construction, are similar to manifold de
Bruijn graphs [24]. As path indexes, they are equivalent
0
• each node has a key G .key(vK ) = K and a value
to de Bruijn graphs with patterns of length up to k
0
G .value(vK ) = VK , where VK is the set of nodes
characters.
S.node(S, i) ∈ V for S ∈ S and positions i such
that S[i, i + |K| − 1] = K; and
Definition 3.4. (Equivalent path graphs) Let G0

00
let k > 0 be a parameter
• each edge (vK , vK 0 ) ∈ E 0 represents the occurrence and G be two path graphs, and
0
value.
We
say
that
graphs
G
and
G00 are k-equivalent,
0
of substring K[0]K in S such that strings K and
0
00
if we have G .locate(X) = G .locate(X) for all patterns
K[0]K 0 prefix-match.
X ∈ (Σ \ {#, $})∗ with 1 ≤ |X| ≤ k.
We use #k as the source node s and $k as the sink node
Definition 3.5. (Pruned de Bruijn graph) Let G
t, adding the edge (t, s) in the usual way.
be a graph, and let G0 be an order-k path graph of G.
Path graph G0 is an order-k pruned de Bruijn graph, if
Definition 3.3. (Path graph as an index)
0
0
0
Let G = (V, E) be a graph, and let G = (V , E ) be a it is k-equivalent to the order-k de Bruijn graph of G.
path graph of G.

ATC
3
GCA
1
3:A
0:#

1:G

2:C

6:C
5:T

4:T

CAT
2, 6

TCA
5

ATA
7

ATG
3

7:A
9:T

10:A

11:$

###
0:2

##G
0:1

#GC
0

TA$
9

A$$
10

$$$
11

TTC
4

8:G
GCT
1

CTT
2

TGT
5

GTA
8

TTG
4

Figure 2: Left: Input graph G = (V, E), with each node v ∈ V labeled with v : G.label(v). Right: The order-3
de Bruijn graph G0 = (V 0 , E 0 ) of graph G, with each node v 0 ∈ V 0 labeled with G0 .key(v 0 ) and G0 .value(v 0 ). Both:
Edges (t, s) are not shown. The highlighted path in the de Bruijn graph is a false positive, as it consists of two
disjoint paths in the input graph.
Lemma 3.4. (No short false positives)
Let G = (V, E) be a graph, let G0 = (V 0 , E 0 ) be an
order-k pruned de Bruijn graph of graph G, and let
X ∈ (Σ \ {#, $})∗ be a pattern with 1 ≤ |X| ≤ k. Then
G0 .locate(X) is a set of start nodes v ∈ V of paths
matching the pattern in graph G.

and outdegrees in bitvectors IN and OUT, and store
the predecessor labels in BWT. See Figure 3 for an
example. If lexicographic range [spi+1 , epi+1 ] matches
suffix X[i + 1, |X| − 1] of pattern X, we can find the
range [spi , epi ] matching suffix X[i, |X| − 1] as

Lemma 3.5. (Pruning) Let G = (V, E) be a graph, let
G0 = (V 0 , E 0 ) be the order-k pruned de Bruijn graph
of G with key set K, let K ∈ Σ∗ be a string of length
|K| > 0, and let VK0 be the set of nodes v 0 ∈ V 0 having
string K as a proper prefix of G0 .key(v 0 ). If |VK0 | > 0 and
G0 .value(u0 ) = G0 .value(v 0 ) for all u0 , v 0 ∈ VK0 , the path
graph with key set (K \ {G0 .key(v 0 ) | v 0 ∈ VK0 }) ∪ {K} is
an order-k pruned de Bruijn graph of G.

[spin , epin ] = [IN.select(spi+1 , 1) + 1,

We can compress a de Bruijn graph structurally
by merging sets of nodes sharing a common prefix of
their keys, as long as the conditions of Lemma 3.5
hold. Let G0 = (V 0 , E 0 ) be an order-k pruned de Bruijn
graph, and let G00 = (V 00 , E 00 ) be the same graph after
further pruning. Each node v 00 ∈ V 00 is an equivalence
class of nodes V 0 (v 00 ) ⊆ V 0 corresponding to a shared
prefix G00 .key(v 00 ) of keys. For all v 0 ∈ V 0 (v 00 ), we have
G0 .value(v 0 ) = G00 .value(v 00 ). See Figure 3 for an example
of a pruned de Bruijn graph.
Definition 3.6. (Maximally pruned graph)
Let G0 be a a pruned de Bruijn graph of graph G. We
say that G0 is maximally pruned, if we cannot prune it
any further using Lemma 3.5.
Lemma 3.6. (Maximal pruning) Let G0 = (V 0 , E 0 )
be a maximally pruned de Bruijn graph of G = (V, E).
Then |G0 .key(u0 )| ≤ |G0 .key(v 0 )| + 1 for all (u0 , v 0 ) ∈ E 0 .
4

GCSA2

As in the original GCSA, we sort the nodes of the
path graph in lexicographic order, encode the indegrees

IN.select(epi+1 + 1, 1)];
[spout , epout ] = [LF(spin , X[i]),
LF(epin + 1, X[i]) − 1];
[spi , epi ] = [OUT.rank(spout , 1),
OUT.rank(epout , 1)].
In order to support locate queries, we sample the
values of a node, if the node has multiple incoming edges
or the values cannot be derived from the predecessor.
We may also sample other values to improve query
performance. The sampled nodes are marked in bitvector
BS , the number of values in each sample is encoded in
unary in bitvector BV , and the sampled values are stored
in integer array VS . A detailed description of the data
structure can be found in Appendix B.
GCSA2 construction starts from paths of length k in
the input graph. We build a maximally pruned order-2k,
order-4k, or order-8k de Bruijn graph using a similar
prefix-doubling algorithm as in the original GCSA [43],
and encode the result as a GCSA. To avoid excessive
memory usage, we keep the paths and the graphs on
disk, and read only a single chromosome at a time into
memory. The details of the construction algorithm can
be found in Appendix C.
We can improve the query performance with (maximally pruned) de Bruijn graphs by using a simplified
encoding (Appendix D). We replace BWT and IN with
bitvectors Bc for all c ∈ Σ, where Bc [j] = 1 if and only
if node j in lexicographic order has a predecessor with

ATC
3
CA
2, 6

TC
5

ATA
7

ATG
3
###
0:2

##G
0:1

#G
0

GC
1

TA
9
CT
2

TT
4

TG
5

A$
10

$$$
11

GT
8

key

OUT

BWT

IN

key

BS

$$$
A$
ATA
ATC
ATG
CA
CT
GC
GT
TA
TC
TG
TT
#G
##G
###

1
1
1
1
1
0
0
1
1
0
1
1
1
1
1
0
1
1
1
1

A
T
C
C
C
G
T
G
#
T
A
G
A
T
A
T
C
#
#
$

1
1
1
1
1
0
1
1
1
1
0
1
0
1
0
1
1
1
1
1

$$$
A$
ATA
ATC
ATG
CA
CT
GC
GT
TA
TC
TG
TT
#G
##G
###

0
0
1
1
1
1
0
0
1
1
1
1
1
0
0
1

BV

VS

1
1
1
0
1
1
1
1
1
1
1

7
3
3
2
6
8
9
5
5
4
0:2

Figure 3: Left: An order-3 pruned de Bruijn graph G00 3-equivalent to the de Bruijn graph in Figure 2. Right:
GCSA for graph G00 . Leftward arrows illustrate backward searching, with the red arrows showing it from T to AT.
Rightward arrows mark the samples belonging to each node, with the blue ones showing them for node CAT.
label c. This simplifies backward searching to
spout = C[X[i]] + BX[i] .rank(spi+1 , 1);
epout = C[X[i]] + BX[i] .rank(epi+1 + 1, 1) − 1;
[spi , epi ] = [OUT.rank(spout , 1), OUT.rank(epout , 1)].
The compacted trie of keys resembles a suffix tree.
We can simulate it space-efficiently by using the LCP
array [1, 11], and thus extend GCSA2 to support many
suffix tree operations. For example, we can search for
maximal exact matches by using LF-mapping and parent
queries [35], and use that as a basis for a read aligner
similar to BWA-MEM [21]. We can also use document
counting techniques [41] to quickly count the number
of distinct matches in a lexicographic range. Further
details of these extensions can be found in Appendix E.
5

Implementation and Experiments

GCSA2 is the path indexing library of vg [14]. The
implementation is written in C++, and the source code
is available on GitHub.3 It depends on SDSL [15] and
libstdc++ parallel mode. We use the simplified encoding
(Appendix D) with fast non-compressed bitvectors
in most index components. Bitvectors Bc for rare
characters (N, #, and $) are compressed as sparse
bitvectors [36].
We used a system with two 16-core AMD Opteron
6378 processors and 256 gigabytes of memory for the
experiments, and stored all files on a distributed Lustre
file system. The system was running Ubuntu 12.04
on Linux kernel 3.2.0. We used vg version 1.3.0 for
processing the graphs and GCSA2 version 0.8 using
SDSL version 2.1.1 for the benchmarks. All code was
compiled with gcc/g++ version 4.9.2.
3 https://github.com/jltsiren/gcsa2

5.1 Construction Variation graphs, as defined in vg,
use strings as node labels. A node can be traversed
in both forward and reverse complement orientations,
and edges may cross between the orientations. For
indexing, the graph is implicitly converted into an input
graph with single-character labels. We always sample
the input graph nodes corresponding to the initial offsets
of variation graph node labels.
We built vg graphs from the human reference genome
(GRCh37) and 1000 Genomes Project variation [45]. To
avoid excessive growth, we removed paths where 16-mers
crossed more than 4 nontrivial edges with vg mod -p
-l -e 4, and subgraphs shorter than 100 bases with vg
mod -S -l 100. We extracted all paths of length 16
from the forward strand of the graph. There were a total
of 4.80 billion paths with 1.53 billion distinct labels. We
then built GCSA with 1–3 doubling steps, producing
order-32, order-64, and order-128 indexes.
Tables 1 and 2 show construction requirements
and index sizes, respectively. We can build a wholegenome index overnight using less than 96 gigabytes
of memory, including disk cache. The index contains
1.031k ·2.348 billion k-mers, but the path graph only uses
4.4–5.7 billion nodes to represent them. For k = 128,
GCSA2 requires 0.63 bits per k-mer, out of which
0.28 bits is used for the path graph. Extensions based
on suffix trees increase the size to 1.08 bits per k-mer.
Index construction uses more memory with k = 32
than with larger values of k. The order-32 path graph
has more nodes, where we cannot derive the values from
the predecessor node. As we sample more values, we
need more memory in the final phase of construction.
With larger values of k, the path graph resembles the
input graph better, and we sample less values. For the
same reason, index size decreases with larger values of
k, even though the graph requires more space.

k k-mers Nodes
Graph
Index
With extensions
32 6.20G 4.37G 2.89 GB / 4.00 bits 9.50 GB / 13.2 bits 13.2 GB / 18.2 bits
64 16.7G 5.24G 3.46 GB / 1.78 bits 8.64 GB / 4.46 bits 13.6 GB / 6.99 bits
5.73G 3.78 GB / 0.28 bits 8.58 GB / 0.63 bits 14.6 GB / 1.08 bits
128 116G
Table 2: GCSA2 index sizes. Order of the path graph; number of k-mers and nodes in the path graph in billions;
index size in gigabytes and in bits per k-mer for the graph (Bc and OUT), the index (the graph, BS , BV , and VS ),
and the index with the extensions from Appendix E.

k
32
64
128

Time Memory Disk
Read
Write
7.44 h 59.8 GB 387 GB 1.37 TB 0.88 TB
10.4 h 51.9 GB 415 GB 2.03 TB 1.51 TB
14.1 h 52.3 GB 478 GB 2.78 TB 2.25 TB

Table 1: GCSA2 construction. Order of the path graph;
construction time in hours; peak memory and disk usage
in gigabytes; and disk I/O volume for reading and writing
in terabytes.

5.2 Queries We compared the query performance of
the order-128 GCSA2 to several FM-indexes for the
reference sequence. SSA is the SDSL implementation
(csa_wt<>) of the succinct suffix array [30], using a
Huffman-shaped wavelet tree on top of the BWT. As
the default FM-index in SDSL, it prioritizes query
performance over compression. We used SSA with SA
sample period 17 for good locate performance. BWA
is the FM-index in the Burrows-Wheeler Aligner [22]
(version 0.7.15 with the default SA sample period 32).
Optimized for DNA sequences, BWA indexes both the
reference and its reverse complement.
As building the original GCSA requires around 65n
bytes of memory for a path graph with n nodes, we could
not compare GCSA and GCSA2 directly on a system
with 256 gigabytes of memory. Instead, we used RLCSA
[31] (May 2016 version) as a proxy. The RLCSA is an
FM-index for repetitive sequence collections using the
same basic components as the original GCSA. Under a
mixed query load, RLCSA with SA sample period 32 is
1.5x to 3x faster than GCSA, depending on algorithmic
overhead and the mix of find and locate queries [43].
We extracted k-mers for k ∈ {16, 32, 64, 128} from
the (non-pruned) vg graphs by using vg sim, filtered
out k-mers consisting entirely of Ns, and queried for the
remaining k-mers using a single thread. The results can
be seen in Table 3.
Backward searching in an FM-index stops early if
there are no matches. In order to compare the find
performance of the indexes reliably, we must hence
concentrate on the 16-mers, where the fraction of

matching patterns is similar for all indexes. GCSA2 and
the fast FM-indexes (SSA and BWA) all have similar
performance, while RLCSA is several times slower. As a
result, we can estimate that find queries in GCSA2 are
an order of magnitude faster than in GCSA.
When comparing the locate performance of different
FM-indexes, the distribution of the query positions
should be close to uniform. Otherwise the biases from
e.g. different suffix array sampling strategies or the
variation in the number of distinct occurrences per node
in GCSA2 can make the results unreliable. As the k-mers
have been sampled uniformly from the variation graph,
we get the best results with the 16-mers, where all indexes
can match most of the patterns.
GCSA2 uses denser SA sampling than the other
indexes, with effective sample period 10.6. On the
average, GCSA2 calls locate for 2.84 nodes per distinct
value, making the amount of work comparable to sample
period 30.2. SSA with sample period 17 is 2.4x faster
than GCSA2, mostly because it has to do less work.
BWA with sample period 32 is closer to GCSA2 in locate
performance. RLCSA is slower than the other indexes,
but the difference is smaller than with find queries due
to the optimizations for retrieving suffix array ranges.
Assuming that locate queries are 3x slower in GCSA than
in RLCSA, as GCSA does not use the optimizations, we
can estimate that GCSA2 is 4x faster than GCSA.
The remaining queries, parent and count, take a
fraction of a microsecond. As a parent query takes
comparable time to a single step of backward searching,
it will not be a bottleneck in finding maximal exact
matches. Counting the number of distinct occurrences
with a count query is faster than retrieving even a single
occurrence.
6

Discussion

GCSA2 is a path index for variation graphs. It uses a
de Bruijn graph as a k-mer index of the variation graph,
prunes it by merging redundant subgraphs, and encodes
the result with a generalization of the FM-index. The
index supports queries of length up to k exactly, and
longer queries with false positives. GCSA2 also includes

Patterns Index
351584 GCSA2
SSA
BWA
RLCSA
32
351555 GCSA2
SSA
BWA
RLCSA
64
351567 GCSA2
SSA
BWA
RLCSA
128 351596 GCSA2
SSA
BWA
RLCSA
k
16

Found Nodes Occs find() parent()
347453 2477M 872M 4.75 μs 0.42 μs
301538
–
782M 6.00 μs
–
320764
–
1564M 3.64 μs
–
301538
–
782M 23.7 μs
–
333258 112M 34.3M 10.8 μs 0.28 μs
153957
–
26.6M 10.9 μs
–
156080
–
52.9M 6.57 μs
–
153957
–
26.6M 47.6 μs
–
326101 2.63M 1.35M 22.5 μs 0.26 μs
88184
–
0.84M 17.1 μs
–
88786
–
1.60M 10.3 μs
–
88184
–
0.84M 74.3 μs
–
316500 0.32M 0.37M 45.3 μs 0.26 μs
35678
–
0.08M 23.5 μs
–
35741
–
0.12M 14.0 μs
–
35678
–
0.08M 91.7 μs
–

count()
0.87 μs
–
–
–
0.38 μs
–
–
–
0.29 μs
–
–
–
0.26 μs
–
–
–

locate()
5.85 μs
2.43 μs
4.65 μs
8.12 μs
5.44 μs
2.16 μs
3.19 μs
5.87 μs
2.92 μs
1.89 μs
2.34 μs
5.97 μs
3.13 μs
3.47 μs
3.46 μs
12.9 μs

Table 3: Query benchmarks using an order-128 GCSA2 and various FM-indexes. Pattern length; number of
patterns; index type; number matching patterns, matching nodes, and distinct occurrences; average time for
[sp, ep] = find(X), parent(sp, ep), and count(sp, ep) queries in microseconds; and average time per value for
locate(sp, ep) queries in microseconds.

extensions based on suffix trees; other extensions have
been considered but not implemented (see Appendix F).
The index is used in the variation graph toolkit vg for
e.g. read alignment based on maximal exact matches.
We can build a whole-genome index overnight on a
system with 96 gigabytes of memory and a few hundred
gigabytes of fast disk space. The resulting index takes
less than 15 gigabytes, or 1.08 bits per k-mer for the
order-128 index with extensions. Query performance is
comparable to that of fast FM-indexes for sequences.
The primary design goals for GCSA2 were query
performance and index size. The index works with
arbitrary graphs, supports queries that are long enough
to map short reads in one piece without false positives,
and provides several options for dealing with complex
regions. Other path indexes work with a more restricted
class of graphs [43, 17, 19, 26], are at least an order
of magnitude slower [17, 26], require much more space
[3, 38, 5, 28], or are theoretical proposals that have never
been implemented [44].
We may want to determine whether a pattern
matches known haplotypes or only their recombinations.
As GCSA2 does not support this directly, vg must
determine it afterwards using a separate structure
[34]. The FM-index of alignment [32, 33] embeds the
haplotypes directly in a GCSA-like index and reports the
haplotypes matching the find query. While the solution
depends on specific properties of the graph, it could be
possible to extend it to work with any GCSA.

Acknowledgements. The author thanks Erik Garrison, Richard Durbin, and Adam Novak for the fruitful
discussions while developing the GCSA2 index.
A

Proofs of Lemmas

Proof. [Lemma 3.1: No false negatives] Let S be the
k-collection used for building the path graph, and let P
be a path starting from v0 ∈ V with G.label(P ) = X.
The collection contains a sequence S ∈ S such that
S[i, i + |X| − 1] = X and S.node(S, i) = v0 .
For all positions j with i ≤ j ≤ i + |X| − 1,
there is a node vj0 ∈ V 0 with G0 .label(vj0 ) = S[j] and
G0 .key(vj0 ) = S[j, j +|G0 .key(vj0 )|−1]. By definition, path
0
graph G0 has an edge (vj0 , vj+1
) ∈ E 0 for all such positions
0
0
0
j. Hence P = vi · · · vi+|X|−1 is a path in G0 with
G0 .label(P 0 ) = X. As path P 0 starts from node vi0 ∈ V 0 ,
node vi0 is included in the set G0 .find(X). Furthermore,
v0 = S.node(S, i) ∈ G0 .value(vi0 ) ⊆ G0 .locate(X).
Proof. [Lemma 3.2: Context length] If |X| ≤ 1, the
statement is true by definition for m = |X|. Now let
Mi+1 = G0 .find(X[i + 1, |X| − 1]) be the set of all nodes
v 0 ∈ V 0 such that substring X[i + 1, i + mi+1 ] is a prefix
of key G0 .key(v 0 ), and assume that the set is nonempty.
Consider the set
[
Mi =
G0 .pred(v 0 , X[i]) = G0 .find(X[i, |X| − 1]).
v 0 ∈Mi+1

There is an edge (u0 , v 0 ) ∈ E 0 if and only if key G0 .key(u0 )

prefix-matches string G0 .label(u0 ) · G0 .key(v 0 ). Hence key
G0 .key(u0 ) prefix-matches string X[i, i + mi+1 ] for all
nodes u0 ∈ Mi .
Now let u0 ∈ V 0 be a node with key G0 .key(u0 ) prefixmatching string X[i, i + mi+1 ]. If the key is a prefix of
string X[i, i + mi+1 ], there is an edge (u0 , v 0 ) ∈ E 0 to
all nodes v 0 ∈ Mi+1 , and hence u0 ∈ Mi . Otherwise let
S[j, j + |G0 .key(u0 )| − 1] = G0 .key(u0 ) be a substring of
the k-collection used for building the path graph. As
S[j + 1, j + mi+1 ] = X[i + 1, i + mi+1 ], the substring
starting at S[j + 1] is represented by a node v 0 ∈ Mi+1 ,
and hence u0 ∈ Mi .
Set Mi is the set of all nodes u0 ∈ V 0 such that
substring X[i, i + mi+1 ] prefix-matches key G0 .key(u0 ).
If |Mi | > 1, string X[i, i + mi+1 ] is a proper prefix of key
G0 .key(u0 ) for all nodes u0 ∈ Mi due to the prefix-free
property, and we can set mi = mi+1 + 1. Otherwise we
set mi = min(mi+1 + 1, |G0 .key(u0 )|) for the only node
u0 ∈ Mi .
Proof. [Lemma 3.3: Short keys] (a) If node v has
multiple predecessors with label c, the keys of the
predecessors must be longer than string c · G0 .key(v 0 ), as
the key set is prefix-free.
(b) By the construction in the proof of Lemma 3.2,
the context length for pattern X in graph G0 is
min(|X|, |G0 .key(v 0 )|).
0

00
or (b) u00 = vK
and u0 ∈ VK0 , and the same holds for
00
0
nodes v and v . We can hence transform any path
in graph G0 into an equivalent path in graph G00 by
00
replacing nodes v 0 ∈ VK0 with node vK
.
00
00
00
Let P = v0 · · · v|P 00 |−1 be a in graph G00 . We
0
transform it into an equivalent path P 0 = v00 · · · v|P
00 |−1
0
00
in graph G . There are two cases for |P | = 1. If
00
v000 = vK
, we can replace it with any v00 ∈ VK0 . Otherwise
we use the node v00 ∈ V 0 with G0 .key(v00 ) = G00 .key(v000 ).
In the general case |P 00 | > 1, assume that we have
00
00
transformed the suffix v100 · · · v|P
into
00 |−1 of path P
0
0
0
an equivalent path v1 · · · v|P 00 |−1 in graph G . Because
G0 .value(v10 ) = G00 .value(v100 ), node v10 must have a
predecessor with label c = G00 .label(v000 ). We can choose
any such predecessor as node v00 .
Consider the predecessors v 0 ∈ G0 .pred(v10 , c) and
00
v ∈ G00 .pred(v100 , c). Their keys prefix-match string
Y = c · K1 . There are three cases:

1. If Y is a prefix of G0 .key(v 0 ), the key G0 .key(w0 ) of
every successor w0 of node v 0 prefix-matches K1 and
hence has K1 as a prefix. Therefore G0 .value(v 0 ) is
the union of sets G.pred(v, c) over v ∈ G0 .value(v10 ).
If Y is also a prefix of G00 .key(v 00 ), nodes v 0 and v 00
have identical value sets by the same reasoning.
2. If Y is a prefix of G0 .key(v 0 ) and G00 .key(v 00 ) is
a proper prefix of Y , there is only one possible
00
predecessor v 00 = v000 . Hence v 0 ∈ VK0 and v000 = vK
.

Proof. [Lemma 3.4: No short false positives] We may
assume without loss of generality that graph G0 is a
3. If key G0 .key(v 0 ) is a proper prefix of string Y , there
de Bruijn graph. Let v 0 ∈ G0 .find(X) be a node. By
is only one possible predecessor v 0 = v00 . Because
Lemma 3.3, pattern X prefix-matches key G0 .key(v 0 ).
key G00 .key(v000 ) prefix-matches string Y , it must
For every node v ∈ G0 .value(v 0 ), there is a substring
also be a proper prefix of the string. Hence either
S[i, i + k − 1] = G0 .key(v00 ) in the k-collection S used
00
.
G0 .key(v00 ) = G00 .key(v000 ) or v00 ∈ VK0 and v000 = vK
0
for building graph G , with S.node(S, i) = v. Prefix
0
0
00
00
S[i, i + |X| − 1] = X of the substring corresponds to a In every case, G .value(v0 ) = G .value(v0 ), and we can
path with label X starting from node v in graph G.
use shorter of strings Y and G00 .key(v000 ) as K0 .
We can transform any path P 0 in graph G0 into
Definition A.1. (Equivalent paths)
an equivalent path P 00 in graph G00 , and the other
0
00
Let G and G be path graphs of the same graph, and way around. Because the labels of equivalent paths
0
00
let P 0 = v00 · · · vn−1
and P 00 = v00 · · · vn−1
be paths in and the value sets of their start nodes are identical,
0
00
graphs G and G , respectively. We say that paths P 0 we have G0 .locate(X) = G00 .locate(X) for all patterns
and P 00 are equivalent, if for 0 ≤ i < n, keys G0 .key(vi0 ) X ∈ (Σ \ {#, $})∗ with |X| > 0.
and G00 .key(vi00 ) have a common prefix Ki such that
G0 .value(v 0 ) = G00 .value(v 00 ) for all nodes v 0 ∈ V 0 and Proof. [Lemma 3.6: Maximal pruning] Let Gd = (Vd , Ed )
the de Bruijn graph of graph G with the same order
v 00 ∈ V 00 having Ki as a prefix of their keys.
k as graph G0 , and let S be the k-collection used for
00
00
00
Proof. [Lemma 3.5: Pruning] Let G = (V , E ) be building the path graphs. If v 0 ∈ V 0 is a node, then
the path graph corresponding to the new key set, let Gd .value(v) = G0 .value(v 0 ) for all nodes v ∈ Vd (v 0 ).
00
vK
∈ V 00 be the node with key K, and let S be the
Assume that |G0 .key(u0 )| > |G0 .key(v 0 )| + 1 for an
k-collection used to define the path graphs.
edge (u0 , v 0 ) ∈ E 0 . String G0 .label(u0 ) · G0 .key(v 0 ) must
0 0
0
Consider the edge (u , v ) ∈ E defined by substring then be a prefix of key G0 .key(u0 ). There cannot be edges
0
S of S. The same substring also defines an edge (u0 , w0 ) to other nodes w0 =
6 v 0 , as keys G0 .key(v 0 ) and
00 00
00
00
00
0
0
0
0
(u , v ) ∈ E , where either (a) G .key(u ) = G .key(u ) G .key(w ) would prefix-match.

Let S[i, i + |G0 .key(u0 )| − 1] = G0 .key(u0 ) be a the range [spi , epi ] = G0 .find(X[i, |X| − 1]):
substring of S. Because G0 .key(v 0 ) is a substring of
[spin , epin ] = [IN.select(spi+1 , 1) + 1,
G0 .key(u0 ) and G0 is a pruned de Bruijn graph, the
IN.select(epi+1 + 1, 1)];
set of nodes S.node(S, i + 1) over all occurrences of
0
0
0
0
0
substring G .key(u ) in S is G .value(v ). As node u has
[spout , epout ] = [LF(spin , X[i]),
no other successors, set G0 .value(u0 ) is the union of sets
LF(epin + 1, X[i]) − 1];
G.pred(v, G0 .label(u0 )) over all nodes v ∈ G0 .value(v 0 ).
[spi , epi ] = [OUT.rank(spout , 1),
The above is true for all x0 ∈ G0 .pred(v 0 , G0 .label(u0 )).
0
OUT.rank(epout , 1)].
Hence we can prune graph G further using string
0
0
0
0
G .label(u ) · G .key(v ) as the new key in Lemma 3.5.
We can think this as a generalization of LF-mapping:
[spi , epi ] = G0 .LF([spi+1 , epi+1 ], X[i]).
B GCSA for Path Graphs
Query G0 .locate(X) retrieves the values G0 .value(v 0 )
0
0
0
Let G = (V , E ) be a path graph. We sort the nodes for nodes v 0 ∈ G0 .find(X) and filters out duplicates.
V 0 by their keys in lexicographic order and generate the Instead of storing the values explicitly for all nodes,
sequences BWT, IN, and OUT from the nodes in that GCSA uses a sampling scheme to save space. We
order. For each node v 0 ∈ V 0 , we append BWT with the assume that the nodes of the input graph G = (V, E) are
predecessor labels G.label(u0 ) for all edges (u0 , v 0 ) ∈ E 0 ; conveniently chosen integers. If (u, v) ∈ E is the only
0
0
IN with the indegree encoded as 0G .in(v )−1 1; and OUT outgoing edge from node u and the only incoming edge
0
0
with the outdegree as 0G .out(v )−1 1.
to node v, it should be that v = u + 1.
If node v 0 ∈ V 0 has lexicographic rank i, the range of
We sample the values G0 .value(v 0 ) for a node v 0 ∈ V 0 ,
0 0
0
incoming edges (u , v ) ∈ E to that node is [spin , epin ] = (a) if there are multiple incoming edges to node v 0 ;
[IN.select(i, 1) + 1, IN.select(i + 1, 1)]. The labels of (b) if v 0 is the source node s; or (c) if G0 .value(v 0 ) 6=
the predecessor nodes are encoded in BWT[spin , epin ]. {u + 1 | u ∈ G0 .value(u0 )} for the only incoming edge
Sorting the incoming edges by pairs (BWT[j], i), where (u0 , v 0 ) ∈ E 0 . We may also sample the values for some
BWT[j] corresponds to edge (u0 , v 0 ) ∈ E 0 , is equivalent nodes on long unary paths for performance reasons. If
to sorting them by strings G0 .label(u0 ) · G0 .key(v 0 ). As the set G0 .value(v 0 ) has not been sampled, we can derive
multiple edges may have the same sort key, our sorting it from sampled values by following the incoming edges.
algorithm must be stable. We get the desired sorting
If node v 0 ∈ V 0 with lexicographic rank i has only
order by using LF-mapping: j 7→ LF(j).
one predecessor, the lexicographic rank of the predecessor
The range of outgoing edges (u0 , v 0 ) ∈ E 0 from node is G0 .LF(i) = OUT.rank(LF(IN.select(i, 1) + 1), 1). If
u0 ∈ V 0 with lexicographic rank i0 is [spout , epout ] = the lexicographic rank G0 .LFk (i) corresponding to node
[OUT.select(i0 , 1)+1, OUT.select(i0 +1, 1)]. The edges are w0 ∈ V 0 is the first sampled node we encounter, we know
already sorted by keys G0 .key(u0 ). Because graph G0 is a that G0 .value(v 0 ) = {w + k | w ∈ G0 .value(w0 )}.
path graph, we know that key G0 .key(u0 ) prefix-matches
Let BS [0, |V 0 |−1] be a bitvector. If we have sampled
0
0
0
0
string G .label(u ) · G .key(v ). The sorting orders are the values for the node v 0 ∈ V 0 with lexicographic
therefore compatible. For every j ∈ [spin , epin ] for a rank i, we mark that as BS [i] = 1. We can then
node v 0 ∈ V 0 , having LF(j) ∈ [spout , epout ] for a node determine the rank of node v 0 among the sampled nodes
u0 ∈ V 0 implies an edge (u0 , v 0 ) ∈ E 0 .
as j = BS .rank(i, 1). For each sampled node v 0 ∈ V , we
0
We use backward searching for query G .find(X). store the size of the value set |G0 .value(v 0 )| in another
Let X ∈ (Σ \ {#, $})∗ be a pattern. If |X| = 0, query bitvector BV , using the same encoding as for bitvectors
G0 .find(X) returns the lexicographic range [0, |V 0 | − 1] IN and OUT. We store the samples in array VS in the
containing all nodes. Now assume that |X| ≥ 1 and same order, using log|V | bits each. The sampled values
that G0 .find(X[i + 1, |X| − 1]) = [spi+1 , epi+1 ]. We want for node v 0 with rank j among the sampled nodes can
to find the lexicographic range G0 .find(X[i, |X| − 1]), be found at VS [BV .select(j, 1) + 1, BV .select(j + 1, 1)].
which is the union of sets G0 .pred(v 0 , X[i]) over nodes
v 0 ∈ G0 .find(X[i + 1, |X| − 1]). We map the node C Index Construction
range [spi+1 , epi+1 ] to the range [spin , epin ] of incoming GCSA construction [43] is based on the prefix-doubling
edges; the incoming edges to the corresponding range of algorithm for suffix array construction [27]. The original
outgoing edges [spout , epout ]; and the outgoing edges to GCSA started from paths of length 1 in the input graph,
and then repeatedly joined paths of length k into paths
of length 2k, until each path had a distinct label. The
resulting path graph was essentially an order-∞ pruned
de Bruijn graph and supported queries of any length.

We use a variant of that algorithm with GCSA2.
Let G = (V, E) be the input graph. We extract all paths
of length k (typically with k = 16) from graph G. For
each path P = v0 · · · v|P |−1 , we store several fields. Key
P.key encodes G.label(P ) as a sequence of lexicographic
ranks of k-mers. If |P | is not an integer multiple of k, the
key consists of the k-mer ranks for the lexicographically
smallest (d|P |/ke · k)-mer having G.label(P ) as a prefix,
followed by the rank of the last k-mer in the largest
such (d|P |/ke · k)-mer. Value P.value is the start node
v0 of the path. We store the set of predecessor labels
{c ∈ Σ | |G.pred(v0 , c)| > 0} as P.pred. For each possible
extension node v ∈ {v ∈ V | (v|P |−1 , v) ∈ E}, we
create a separate copy of the path and store the node as
P.ext = v.
The construction uses several supporting structures.
We build an order-k de Bruijn graph Gd = (Vd , Ed ) of
the path labels and encode it as a GCSA, using the
predecessor labels P.pred for determining the edges. Let
v0 , . . . , v|Vd |−1 be the nodes of the de Bruijn graph in
lexicographic order by their keys. We use two additional
arrays: the LCP array LCP[0, |Vd | − 1], where LCP[i]
is the length of the longest common prefix of keys
Gd .key(vi−1 ) and Gd .key(vi ) (with LCP[0] = 0), and
the last character array L[0, |Vd | − 1], where L[i] =
Gd .key(vi )[k − 1]. The LCP array is stored as a wavelet
tree for fast range minimum queries [13].
Because we store path labels explicitly, we only do a
limited number of doubling steps, typically 2 or 3. After
d doubling steps, the length of the paths is 2d k, and we
can use them to build a maximally pruned order-(2d k)
de Bruijn graph. Each doubling step consists of a pruning
step, followed by an extension step. The pruning step
applies a limited form of Lemma 3.5 to lexicographic
ranges of paths. Given two paths P and P 0 , we can
determine the length of the longest common prefix of
the path labels by using the keys P.key and P 0 .key and
the LCP array. If all paths sharing a prefix start from
the same node, we merge them into a single path Q with
Q.key based on the shared prefix and Q.ext = −1.
The extension step transforms the current set of
paths of length (up to) k 0 into a set of paths of length
(up to) 2k 0 . If P is a path with P.ext = −1, we use
it as such. If P.ext = P 0 .value for paths P and P 0 , we
create a new path P P 0 . We set (P P 0 ).key according to
the concatenation of the path labels, take value and pred
from path P , and take ext from path P 0 . If we have
another path Q with Q.key = P.key such that QP 0 is a
path, and if P 0 .ext = −1, all possible 2k 0 -mer extensions
of label G.label(P P 0 ) are also labels of paths starting
from node Q.value, and the other way around. Hence
paths P P 0 and QP 0 can be represented by a single node
in a pruned de Bruijn graph.

The doubling steps are followed by the merging
step, which transforms the paths into the nodes of a
maximally pruned de Bruijn graph G00 = (V 00 , E 00 ).
We merge the paths with identical keys into the nodes
of a pruned de Bruijn graph G0 = (V 0 , E 0 ). If paths
P0 , · · · , Pm−1 all have the same key, we create a node
v 0 ∈ V 0 with the shared key as G0 .key(v 0 ) and with
Sm−1
G0 .value(v 0 ) = i=0 Pi .value. We also store the union
Sm−1
of predecessor labels as v 0 .pred = i=0 Pi .pred. We
then apply Lemma 3.5 maximally, transforming graph
G0 into graph G00 .
Storing the paths and the graphs may require
hundreds of gigabytes of memory when indexing wholegenome variation graphs. To avoid that, we keep them
on disk when possible. The subgraph corresponding to
each chromosome is stored in a separate file, with paths
sorted by their labels in lexicographic order. Extension
steps are done separately for each of the chromosomes.
The pruning step merges the sorted files. It keeps reading
paths into a buffer, until it has found a maximal range of
paths that can be merged. The merged path is written
into the new file for that chromosome, and the original
range of paths is removed from the buffer. The merging
step works in a similar way.
After creating the nodes V 00 of the maximally pruned
order-k 0 de Bruijn graph G00 = (V 00 , E 00 ), we build
the index. Sequences BWT and IN can be generated
from the predecessor sets v 00 .pred, while the outdegree
sequence OUT requires further processing. There is an
edge (u00 , v 00 ) ∈ E 00 if and only if string c · G00 .key(v 00 )
prefix-matches key G00 .key(u00 ) and c ∈ v 00 .pred. We
determine the edges and produce the outdegree sequence
by scanning the node file sequentially with σ + 1 pointers.
The pointer corresponding to node v 00 scans the entire
file, while each of the remaining σ pointers scans only the
range of nodes u00 with G00 .label(u00 ) = c for a character
c ∈ Σ. We also sample the nodes for locate queries
during the same scans.
Checking whether key G00 .key(u00 ) prefix-matches
string c · G00 .key(v 00 ) can be done by using the GCSA for
de Bruijn graph Gd = (Vd , Ed ) and the last character
array L. If the lexicographic rank of k-mer X is i, the
lexicographic rank of k-mer (cX)[0, k − 1] is Gd .LF(i, c).
If 2k-mer X is encoded with k-mer ranks (i, j), we can encode string cX as a lexicographic range of k-mer rank sequences, with (Gd .LF(i, c), Gd .LF(j, L[i]), Gd .LF(0, L[j]))
as the lower bound and Gd .LF(|Vd | − 1, L[j]) as the last
rank in the upper bound. String c · G00 .key(v 00 ) prefixmatches key G00 .key(u00 ) if and only if the lexicographic
ranges of the k-mer rank sequences overlap.
Complex regions of the variation graph must be
pruned before indexing. While this happens before index
construction begins, the construction algorithm has fea-

tures that can make the pruning less destructive. Pruning heuristics often create regions that are completely
missing from the index. The same mechanism that saves
memory by having each chromosome in a separate file
can be used to index overlapping subgraphs without indexing any paths between them. By having the pruned
graph in one file and the reference path in another file,
we can guarantee that no region is completely missing
from the index. We can also index selected paths in complex regions by duplicating nodes for prefix-doubling and
mapping the duplicates back to the original nodes during
the merging step. If we index the paths corresponding to
known haplotypes in complex regions, we can guarantee
that the index contains all observed variation.
D

Simplified GCSA Encoding

Let G0 = (V 0 , E 0 ) be a path graph with G0 .pred(v 0 , c) ≤ 1
for all nodes v 0 ∈ V 0 and characters c ∈ Σ. This is true
for de Bruijn graphs, and also for maximally pruned
de Bruijn graphs (Lemmas 3.6 and 3.3). We can use the
simplified encoding of the original GCSA [43] with such
path graphs.
We replace the sequences BWT and IN with indicator
bitvectors Bc [0, |V 0 | − 1] for all c ∈ Σ. If node v 0 ∈ V
with lexicographic rank i has a predecessor with label
c ∈ Σ, we set Bc [i] = 1. The backward step becomes:
spout = C[X[i]] + BX[i] .rank(spi+1 , 1);
epout = C[X[i]] + BX[i] .rank(epi+1 + 1, 1) − 1;
[spi , epi ] = [OUT.rank(spout , 1), OUT.rank(epout , 1)].
Two expensive queries (IN.select() and BWT.rank()) are
replaced with a cheap Bc .rank().
Computing G0 .LF(i) can expensive, as we have to
look at Bc [i] for all c ∈ Σ to determine the character
used in the backward step. If the alphabet is small,
this is still faster than the select queries in the general
encoding. We can further improve the time/space tradeoff by compressing the bitvectors Bc for rare characters
(e.g. N, #, and $) and checking Bc [i] first for the frequent
characters (e.g. bases) when computing G0 .LF(i).
E

Suffix Tree of a Path Graph

resented by the lexicographic range [sp, ep] matching the
path label `(v) from the root to the node. These ranges
can be determined from the LCP array. If we use the
FM-index with an LCP array supporting next/previous
smaller value queries and range minimum queries, we
can support the full functionality of the suffix tree [11].
If we build an x-ary tree over the LCP array, with
each internal node storing the minimum LCP value in
the corresponding range, we can support the required
queries in O(x logx |V 0 |) time with O(logx |V 0 |) random
x
memory accesses, while using x−1
|V 0 | log k bits of space.
Because the nodes of a path graph may match patterns
that do not prefix-match their keys, we have to be careful
with the suffix tree operations we use. By Lemma 3.2,
the ranges returned by find queries always correspond to
prefixes of the pattern. Hence we can safely use parent
queries with such ranges.4
BWA-MEM [21] aligns reads to a reference genome
by finding maximal exact matches between the read
and the reference. It indexes both the reference and its
reverse complement, allowing it to extend the pattern in
both directions. GCSA cannot use similar techniques, as
we cannot guarantee that the length of the lexicographic
range matching pattern X is the same as the length of
←
−
the range matching its reverse complement X (that the
←
−
key set contains key K if and only if it contains key K ).
However, we can search for maximal exact matches by
using LF-mapping and the parent operation [35].
In an ordinary FM-index, the length of a lexicographic range tells the number of distinct pointers in the
range. In GCSA, each node v 0 ∈ V 0 may have multiple
values (pointers), and a value may occur in multiple
nodes. As locate queries can be slow, we need another
way to support efficient counting queries. The problem
is similar to determining the document frequency of a
pattern. Given a collection of documents, the document
frequency of pattern X is the number of distinct documents that contain occurrences of the pattern. We can
use a bitvector of length 2n − d − 1, where n is the total
size of the document collection and d is the number of
documents, to compute the frequencies efficiently [41].
Let v be a suffix tree node corresponding to lexicographic range [sp, ep], and let count(v) be the document
frequency of the label `(v). If nodes v0 , . . . , vm−1 are the
children of node v,Pthe number of redundant documents
m−1
in them is R(v) = i=0 count(vi )−count(v). We create
an array R[0, n − 2] based on the inorder traversal of the
suffix tree. If the ith internal node we encounter is our
first visit to node v, we set R[i] = R(v). We set R[j] = 0
for any subsequent visits to the same node. Range

Let G0 = (V 0 , E 0 ) be an order-k path graph, and let
0
v00 , . . . , v|V
0 |−1 be its nodes in lexicographic order. The
LCP array of graph G0 is an array LCP[0, |V 0 |−1], where
LCP[i] is the length of the longest common prefix of keys
0
G0 .key(vi−1
) and G0 .key(vi0 ), with LCP[0] = 0. If we
build a trie of keys G0 .key(v 0 ) for v 0 ∈ V 0 and compact
the unary paths into single edges, the resulting tree can
be considered the suffix tree of graph G0 .
4 The shorter query in the variable-order de Bruijn graph [2] is
The LCP interval tree [1] is an alternate representation of the suffix tree. Each node v of the suffix tree is rep- essentially a parent query.

R[sp, ep − 1] covers the internal nodes in the subtree
with node v as the root. We can determine
P document
frequencies as count(v) = (ep + 1 − sp) − R[sp, ep − 1].
If we encode array R in unary, with value x becoming
x
0P
1, we get a bitvector BR , where we can compute sums
b
i=a R[i] with select queries as
(BR .select(b + 1, 1) − b) − (BR .select(a, 1) + 1 − a).
For value counting, we use array R to store the
number of redundant values as above. We use another
array A[0, |V 0 | − 1] to store the number of additional
values in each node vi0 ∈ V 0 as A[i] = |G0 .value(vi0 )| − 1,
and encode it as a bitvector BA in the same way as
array R above. The number of distinct values in range
[sp, ep] = G0 .find(X) for a pattern X is
G0 .count(X) = G0 .count(sp, ep)
=

ep
X

(A[i] + 1) −

i=sp

ep−1
X

R[i].

While graphs are a natural formalism for representing genetic variation, they cannot adequately represent
certain types of rearrangements. For example, if sequence S can occur in different positions of the genome
(e.g. ASBC and ABSC), we can either have a single
copy or multiple copies of S in the graph. Neither option
is good in a reference genome. With a single copy, we
create paths in the reference that do not correspond to
any valid genome. With the second option, we lose the
information that both copies of S are the same sequence.
One solution is to use context-free grammars. As long
as the grammar is non-nested, we can handle it with
the hypertext index. We build a hypertext index for
a high level graph, where each node is labeled with a
nonterminal symbol, while each nonterminal expands
into a subgraph indexed in GCSA.
References

i=sp

The bitvectors are often highly compressible [12], but
GCSA already uses one of the compression schemes
implicitly when it prunes the de Bruijn graph.
F Using the Hypertext Index
The hypertext index [44] is based on graphs G = (V, E),
where the label G.label(v) of a node v ∈ V is a string
over alphabet Σ. The labels are indexed in FM-index
F , while the reverse labels are indexed in FM-index
R. When we search for a pattern X ∈ Σ∗ of length
|X| ≥ 2, some of the matches may cross edges. In
order to find matches crossing one edge, we search for
suffixes X[i, |X| − 1] in the forward index F and the
reverses of prefixes X[0, i − 1] in the reverse index R,
for all 1 ≤ i ≤ |X| − 1. For each value of i, we combine
the partial matches into complete matches with a twodimensional range query in the edge matrix E, using the
lexicographic ranges for the reverse of X[0, i − 1] and for
X[i, |X| − 1] as the query ranges.
We often have to prune complex regions of the input
graph before indexing it. This causes false negatives:
paths that exist in the input graph but not in the index.
We can avoid the false negatives with a generalization of
the hypertext index. Instead of pruning the input graph
heuristically, we create a primary graph based on known
haplotypes and build a GCSA index for both strands
of the graph. We then create a matrix of additional
edges corresponding to potential recombinations in the
path graph, always crossing from the reverse complement
strand to the forward strand. We search for pattern X
←
−
and its reverse complement X in the index, and combine
←
−
the results find(X[i, |X|−1]) and find( X [|X|−i, |X|−1])
with a range query.

[1] M. I. Abouelhoda, S. Kurtz, and E. Ohlebusch. Replacing suffix trees with enhanced suffix arrays. Journal of
Discrete Algorithms, 2(1):53–86, 2004.
[2] C. Boucher et al. Variable-order de Bruijn graphs. In
Proc. 2015 IEEE Data Compression Conference (DCC
2015), pages 383–392. IEEE, 2015.
[3] A. Bowe, T. Onodera, K. Sadakane, and T. Shibuya.
Succinct de Bruijn graphs. In Proc. 12th Workshop
on Algorithms in Bioinformatics (WABI 2012), volume
7534 of LNCS, pages 225–235. Springer, 2012.
[4] M. Burrows and D. J. Wheeler. A block sorting lossless
data compression algorithm. Technical Report 124,
Digital Equipment Corporation, 1994.
[5] B. Cazaux, T. Lecroq, and E. Rivals. From indexing
data structures to de bruijn graphs. In Proc. 25th
Annual Symposium on Combinatorial Pattern Matching
(CPM 2014), volume 8486 of LNCS, pages 89–99.
Springer, 2014.
[6] D. M. Church et al. Extending reference assembly
models. Genome Biology, 16:13, 2015.
[7] A. Danek, S. Deorowicz, and S. Grabowski. Indexing
large genome collections on a PC. PLOS ONE,
9(10):e109384, 2014.
[8] A. Dilthey et al. Improved genome inference in the
MHC using a population reference graph. Nature
Genetics, 47(6):682–688, 2015.
[9] P. Ferragina, F. Luccio, G. Manzini, and S. Muthukrishnan. Compressing and indexing labeled trees, with
applications. Journal of the ACM, 57(1):4, 2009.
[10] P. Ferragina and G. Manzini. Indexing compressed text.
Journal of the ACM, 52(4):552–581, 2005.
[11] J. Fischer, V. Mäkinen, and G. Navarro. Faster entropybounded compressed suffix trees. Theoretical Computer
Science, 410(51):5354–5364, 2009.
[12] T. Gagie et al. Document counting in compressed space.
In Proc. 2015 IEEE Data Compression Conference
(DCC 2015), pages 103–112. IEEE, 2015.

[13] T. Gagie, G. Navarro, and S. J. Puglisi. New algorithms on wavelet trees and applications to information
retrieval. Theoretical Computer Science, 426-427:25–41,
2012.
[14] E. Garrison et al. vg: the variation graph toolkit, 2014–
2016. https://github.com/vgteam/vg.
[15] S. Gog, T. Beller, A. Moffat, and M. Petri. From theory
to practice: Plug and play with succinct data structures.
In Proc. 13th International Symposium on Experimental
Algorithms (SEA 2014), volume 8504 of LNCS, pages
326–337. Springer, 2014.
[16] R. Grossi, A. Gupta, and J. S. Vitter. High-order
entropy-compressed text indexes. In Proc. fourteenth
annual ACM-SIAM symposium on Discrete algorithms
(SODA 2003), pages 841–850. SIAM, 2003.
[17] L. Huang, V. Popic, and S. Batzoglou. Short read
alignment with populations of genomes. Bioinformatics,
29(13):i361–i370, 2013.
[18] S. Huang et al. Indexing similar DNA sequences. In
Proc. The Sixth International Conference on Algorithmic Aspects in Information and Management (AAIM
2010), volume 6124 of LNCS, pages 180–190. Springer,
2010.
[19] D. Kim et al. HISAT2, 2015–2016. https://ccb.jhu.
edu/software/hisat2/index.shtml.
[20] D. Kim, B. Langmead, and S. L. Salzberg. HISAT:
a fast spliced aligner with low memory requirements.
Nature Methods, 12(4):357–360, 2015.
[21] H. Li. Aligning sequence reads, clone sequences and
assembly contigs with BWA-MEM. arXiv:1303.3997,
2013.
[22] H. Li and R. Durbin. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics,
25(14):1754–1760, 2009.
[23] A. Limasset, B. Cazaux, E. Rivals, and P. Peterlongo.
Read mapping on de bruijn graphs. BMC Bioinformatics, 17:237, 2016.
[24] Y. Lin and P. A. Pevzner. Manifold de Bruijn graphs. In
Proc. 14th Workshop on Algorithms in Bioinformatics
(WABI 2014), volume 8701 of LNCS, pages 296–310.
Springer, 2014.
[25] B. Liu, H. Guo, M. Brudno, and Y. Wang. deBGA:
read alignment with de Bruijn graph-based seed and
extension. Accepted to Bioinformatics, 2016.
[26] S. Maciuca, C. del Ojo Elias, G. McVean, and Z. Iqbal.
A natural encoding of genetic variation in a BurrowsWheeler transform to enable mapping and genome
inference. bioRxiv, 2016.
[27] U. Manber and G. Myers. Suffix arrays: A new method
for on-line string searches. SIAM Journal on Computing,
22(5):935–948, 1993.
[28] S. Marcus, H. Lee, and M. C. Schatz. SplitMEM: a
graphical algorithm for pan-genome analysis with suffix
skips. Bioinformatics, 30(24):3476–3483, 2014.
[29] T. Marschall et al. Computational pan-genomics:
Status, promises and challenges. BioRxiv, 2016.
[30] V. Mäkinen and G. Navarro. Succinct suffix arrays
based on run-length encoding. Nordic Journal of

Computing, 12(1):40–66, 2005.
[31] V. Mäkinen, G. Navarro, J. Sirén, and N. Välimäki.
Storage and retrieval of highly repetitive sequence collections. Journal of Computational Biology, 17(3):281–308,
2010.
[32] J. C. Na et al. FM-index of alignment: A compressed
index for similar strings. Theoretical Computer Science,
638:159–170, 2016.
[33] J. C. Na et al. FM-index of alignment with gaps.
arXiv:1606.03897, 2016.
[34] A. M. Novak, E. Garrison, and B. Paten. A graph
extension of the positional Burrows-Wheeler transform
and its applications. bioRxiv, 2016.
[35] E. Ohlebusch, S. Gog, and A. Kügel. Computing
matching statistics and maximal exact matches on
compressed full-text indexes. In Proc. 17th Symposium
on String Processing and Information Retrieval (SPIRE
2010), volume 6393 of LNCS, pages 347–358. Springer,
2010.
[36] D. Okanohara and K. Sadakane. Practical entropycompressed rank/select dictionary. In Proc. Ninth
Workshop on Algorithm Engineering and Experiments
(ALENEX 2007), pages 60–70. SIAM, 2007.
[37] B. Paten, A. Novak, and D. Haussler. Mapping to a
reference genome structure. arXiv:1404.5010, 2014.
[38] J. Pell et al. Scaling metagenome sequence assembly
with probabilistic de Bruijn graphs. Proc. National
Academy of Sciences of the United States of America,
109(33):13272–13277, 2012.
[39] E. A. Rødland. Compact representation of k-mer
de Bruijn graphs for genome read assembly. BMC
Bioinformatics, 14:313, 2013.
[40] K. Sadakane. Compressed suffix trees with full functionality. Theory of Computing Systems, 41(4):589–607,
2007.
[41] K. Sadakane. Succinct data structures for flexible
text retrieval systems. Journal of Discrete Algorithms,
5(1):12–22, 2007.
[42] K. Schneeberger et al. Simultaneous alignment of
short reads against multiple genomes. Genome Biology,
10(9):R98, 2009.
[43] J. Sirén, N. Välimäki, and V. Mäkinen. Indexing graphs
for path queries with applications in genome research.
IEEE/ACM Transactions on Computational Biology
and Bioinformatics, 11(2):375–388, 2014.
[44] C. Thachuk. Indexing hypertext. Journal of Discrete
Algorithms, 18:113–122, 2013.
[45] The 1000 Genomes Project Consortium. A global
reference for human genetic variation. Nature, 526:68–
64, 2015.
[46] S. Wandelt, J. Starlinger, M. Bux, and U. Leser. RCSI:
Scalable similarity search in thousand(s) of genomes.
Proc. VLDB Endowment, 6(13):1534–1545, 2013.
[47] P. Weiner. Linear pattern matching algorithms. In
Proc. 14th Annual IEEE Symposium on Switching and
Automata Theory (FOCS 1973), pages 1–11. IEEE,
1973.

