Graph Indexing Beyond Wheeler Graphs

Let $S=S[1]\mathrm{\cdots }S[n]$ be a string of length $\left|S\right|$ over an alphabet $\mathrm{\Sigma }$ of size $\left|\mathrm{\Sigma }\right|$. If the string and the alphabet are clear from the context, we often denote their sizes as $n=\left|S\right|$ and $\sigma =\left|\mathrm{\Sigma }\right|$. We write substrings as $S[i:j]=S[i]\mathrm{\cdots }S[j]$ and call substrings of the form $S[:j]=S[1:j]$ and $S[i:]=S[i:n]$ prefixes and suffixes, respectively. Let $\$\in \mathrm{\Sigma }$ be the smallest symbol in the alphabet. We often use it as a sentinel $S[n]=\$$ that occurs only at the end of the string.

Let $S$ be a string of length $n$, let $i\in \mathbb{Z}$ be an offset in it, and let $c\in \mathrm{\Sigma }$ be a symbol. We define $S.\mathrm{rank}(i,c)$ as the number of occurrences of symbol $c$ in the prefix $S[1:i]$. For convenience, we define $S.\mathrm{rank}(i,c)=0$ when and $S.\mathrm{rank}(i,c)=S.\mathrm{rank}(n,c)$ when $i>n$.

We define the inverse operation $S.\mathrm{select}(i,c)$ as the position of the occurrence of $c$ of rank $i$. Let $n_c=S.\mathrm{rank}(n,c)$ be the total number of occurrences of symbol $c$ in the string. For $1\le i\le n_c$, we have $j=S.\mathrm{select}(i,c)$ such that $S.\mathrm{rank}(j,c)=i$ and $S.\mathrm{rank}(j-1,c)=i-1$. We also define $S.\mathrm{select}(0,c)=0$ and $S.\mathrm{select}(n_c+1,c)=n+1$ for convenience and leave the function undefined for the remaining values of $i$.

A bitvector is a data structure that stores a binary string and supports efficient rank/select queries. If we store the string as such, we can support both queries in $\mathrm{O}\left(1\right)$ time and $n+\mathrm{o}\left(n\right)$ bits of space [14, 37]. When the bitvector is sparse, we can achieve better space usage by using Elias–Fano encoding to store the positions with ones [31, 42]. If there are $n_1$ ones in a bitvector of length $n$, we can store it in $n_1\log (n/n_1)+\mathrm{O}\left(n_1\right)$ bits of space and support select queries in $\mathrm{O}\left(1\right)$ time and rank queries in $\mathrm{O}\left(\log (n/n_1)\right)$ time.

Let $𝒮=(S_1,\mathrm{\dots },S_m)$ be an ordered collection of $\left|𝒮\right|$ strings of total length $\Vert 𝒮\Vert$. We assume that each string $S_i$ ends with a sentinel ${\$}_i$ and that iff to make the suffixes of the strings all distinct in lexicographic comparisons. However, we represent the sentinels using the same symbol $\$$ in the actual strings. When the collection is clear from the context, we use $m=\left|𝒮\right|$ and $N=\Vert 𝒮\Vert$ to denote the number of strings and the total length of the strings. When the collection consists of a single string $S$, we may consider the collection and the string to be the same.

Let $(i,j)$ and $(i^{\prime },j)$ be two positions in the same column. If we have $\mathrm{𝖬𝖲𝖠}[i,j]\ne -$ (if $\mathrm{𝖬𝖲𝖠}[i^{\prime },j]\ne -$), let $S_i[x]$ ($S_{i^{\prime }}[y]$) be the corresponding symbol. We interpret the alignment in the following way:

• $\mathrm{■}$

  If $\mathrm{𝖬𝖲𝖠}[i,j]=\mathrm{𝖬𝖲𝖠}[i^{\prime },j]\ne -$, positions $S_i[x]$ and $S_{i^{\prime }}[y]$ are a match.

• $\mathrm{■}$

  If $-\ne \mathrm{𝖬𝖲𝖠}[i,j]\ne \mathrm{𝖬𝖲𝖠}[i^{\prime },j]\ne -$, positions $S_i[x]$ and $S_{i^{\prime }}[y]$ are a mismatch.

• $\mathrm{■}$

  If $\mathrm{𝖬𝖲𝖠}[i,j]\ne -$ and $\mathrm{𝖬𝖲𝖠}[i^{\prime },j]=-$, string $S_i$ has an insertion relative to string $S_{i^{\prime }}$ and string $S_{i^{\prime }}$ has a deletion relative to string $S_i$ at the position.

Recombinations let us build new strings by breaking existing strings in two at aligned positions and concatenating the prefix of one string with the suffix of another. Using the above notation, if positions $S_i[x]$ and $S_{i^{\prime }}[y]$ are a match, we can recombine strings $S_i$ and $S_{i^{\prime }}$ at those positions as $S_i[:x-1]S_{i^{\prime }}[y:]$ and $S_{i^{\prime }}[:y-1]S_i[x:]$. This corresponds to switching between rows $i$ and $i^{\prime }$ in column $j$ of the multiple sequence alignment.

In later sections, we discuss data structures for indexing all possible recombinations of a string collection compatible with a given alignment. We formalize this as follows.

According to rule 3, each suffix $S^{\prime \prime }[j^{\prime \prime }:]$ of the recombined string $S^{\prime \prime }$ inherits the alignment of the suffix starting from the corresponding position in the source strings $S$ and $S^{\prime }$.

Let $𝒮$ be an ordered collection of strings. We define the (generalized) suffix array of collection $𝒮$ as an array of pointers $\mathrm{SA}[1:N]$ to the suffixes of the strings in lexicographic order. Each pointer $\mathrm{SA}[x]=(i,j)$ refers to the suffix $S_i[j:]$ of string $S_i$. Given two pointers $\mathrm{SA}[x]=(i,j)$ and $\mathrm{SA}[x^{\prime }]=(i^{\prime },j^{\prime })$, we have iff in lexicographic order. We therefore call $x$ the lexicographic rank of suffix $S_i[j:]$ among the suffixes of collection $𝒮$. Note that for $1\le x\le m$, we have $\mathrm{SA}[x]=(x,\left|S_x\right|)$.

The multi-string Burrows–Wheeler transform (BWT) [11] of collection $𝒮$ is a string $\mathrm{BWT}[1:N]$ representing a permutation of the symbols in strings $S\in 𝒮$. We define it using the suffix array. If $\mathrm{SA}[x]=(i,j)$ and $j>1$, the BWT stores the preceding symbol $\mathrm{BWT}[x]=S_i[j-1]$. When $j=1$, we set $\mathrm{BWT}[x]=\$$.

In addition to the suffix array and the Burrows–Wheeler transform, we often define other arrays over the lexicographically sorted suffixes. The longest common prefix (LCP) array $\mathrm{LCP}[1:N]$ stores the length of the longest common prefix of two lexicographically adjacent suffixes. If $\mathrm{SA}[x]=(i,j)$ and $\mathrm{SA}[x-1]=(i^{\prime },j^{\prime })$, we store the LCP of suffixes $S_i[j:]$ and $S_{i^{\prime }}[j^{\prime }:]$ in $\mathrm{LCP}[x]$, with $\mathrm{LCP}[1]=0$.

The Burrows–Wheeler transform can be inverted using a function called LF-mapping that maps the lexicographic rank of a suffix to the lexicographic rank of its cyclic predecessor. We define it in two parts. Let $\mathrm{SA}[x]=(i,j)$. If $j>1$, we define $\mathrm{LF}(x)=y$ such that $\mathrm{SA}[y]=(i,j-1)$. We sometimes write this as $\mathrm{SA}[\mathrm{LF}(x)]=\mathrm{SA}[x]-1$. For $j=1$, we define $\mathrm{LF}(x)=y$ such that $\mathrm{SA}[y]=(i,\left|S_i\right|)$.

We can compute the LF-mapping for most values of $x$ using the BWT. Let $\mathrm{C}[1:\sigma ]$ be an array such that $\mathrm{C}[c]$ is the total number of occurrences of all symbols in collection $𝒮$. Let $\mathrm{SA}[x]=(i,j)$. If $j>1$, which is equivalent to $\mathrm{BWT}[x]\ne \$$, we have $\mathrm{LF}(x)=\mathrm{C}[\mathrm{BWT}[x]]+\mathrm{BWT}.\mathrm{rank}(x,\mathrm{BWT}[x])$. Here $\mathrm{C}[\mathrm{BWT}[x]]$ is the number of suffixes starting with a symbol smaller than $\mathrm{BWT}[x]=S_i[j-1]$ and $\mathrm{BWT}.\mathrm{rank}(x,\mathrm{BWT}[x])$ is the lexicographic rank of suffix $S_i[j-1:]$ among the suffixes starting with symbol $\mathrm{BWT}[x]$.

Inverse LF-mapping, denoted as $\mathrm{\Psi }$, can be computed as $\mathrm{\Psi }(x)={\mathrm{LF}}^{-1}(x)$ by identifying a symbol $c$ such that and then applying $\mathrm{\Psi }(x)=\mathrm{BWT}.\mathrm{select}(x-\mathrm{C}[c],c)$.

We can generalize LF-mapping for hypothetical suffixes. Let $S$ be the lexicographic rank of string $X$ among the suffixes of collection $𝒮$. That is, $x$ is the number of suffixes $Y$ of the collection such that $Y\le X$. Then, given a symbol $c\in \mathrm{\Sigma }$ with $c\ne \$$, we define the generalized LF-mapping as the lexicographic rank of string $cX$. We can compute it as $\mathrm{LF}(x,c)=\mathrm{C}[c]+\mathrm{BWT}.\mathrm{rank}(x,c)$.

Let $\mathrm{find}(X)=[sp:ep]$ be the lexicographic range of suffixes starting with string $X$. That is, for each $x\in [sp:ep]$, we have $\mathrm{SA}[x]=(i,j)$ such that $S_i[j:j+|X|-1]=X$. Then for any symbol $c$, the range of suffixes starting with $cX$ is $\mathrm{find}(cX)=[\mathrm{LF}(sp-1,c)+1:\mathrm{LF}(ep,c)]$. This backward search step forms the core of the FM-index [22]. Given a pattern $P$, we can determine the range of suffixes $\mathrm{find}(P)$ starting with $P$ with $\mathrm{O}\left(\left|P\right|\right)$ rank operations on the BWT.

Another core operation is $\mathrm{locate}(x)=\mathrm{SA}[x]=(i,j)$, which converts the lexicographic rank of a suffix to its starting position. To compute it, the FM-index needs suffix array samples in addition to the BWT. Let $B[1:N]$ be a bitvector that marks the sampled suffix array positions with $1$, and let $S$ be the subsequence of the suffix array corresponding to those positions. If $B[x]=1$, we can compute $\mathrm{locate}(x)$ as $S[B.\mathrm{rank}(x,1)]$. Otherwise we rely on iterated LF-mapping and the fact that $\mathrm{SA}[\mathrm{LF}(x)]=\mathrm{SA}[x]-1$ in the general case. If we need $k$ steps to find the nearest sampled position, we have $\mathrm{locate}(x)=\mathrm{locate}({\mathrm{LF}}^k(x))+k$. For this scheme to function, we need to sample $\mathrm{SA}[x]=(i,j)$ whenever $j=1$. If we sample $\mathrm{SA}[x]$ every time $j\equiv 1(mods)$ for a parameter $s>0$, we can compute $\mathrm{locate}(x)$ with $\mathrm{O}\left(s\right)$ rank operations.

Let $G=(V,E)$ be a graph with a finite set of vertices $V$ and a set of edges $E\subseteq V\times V$. We assume that the graph is directed. For any two vertices $u,v\in V$, with $u\ne v$, edges $(u,v)$ and $(v,u)$ are distinct. A path is a string $P[1:n]$ of vertices such that $(P[i],P[i+1])\in E$ for . We call the graph acyclic if there is no path that visits a vertex $v\in V$ more than once.

A vertex-labeled graph is a tuple $G=(V,E,\mathrm{\Sigma },\mathrm{\ell })$, where $(V,E)$ is a graph, $\mathrm{\Sigma }$ is an alphabet, and $\mathrm{\ell }:V\to \mathrm{\Sigma }$ is a function that assigns a label $\mathrm{\ell }(v)\in \mathrm{\Sigma }$ to each vertex $v\in V$. If $P$ is a path, we write $\mathrm{\ell }(P)=\mathrm{\ell }(P[1])\mathrm{\cdots }\mathrm{\ell }(P[n])$ to refer to the string formed by concatenating the vertex labels. We can similarly define edge-labeled graphs as tuples $G=(V,E,\mathrm{\Sigma },\mathrm{\ell })$, where $\mathrm{\ell }:E\to \mathrm{\Sigma }$ is a function that assigns labels to edges. Later, we extend the notions to vertices and edges labeled with non-empty strings.

Note that the set of labels $\mathrm{\ell }(P)$ for paths $P$ in the alignment graph is the set of substrings of strings $S\in {𝒮}^{\sim }$. If $\sim$ is an equivalence relation based on a multiple sequence alignment such that $S[j:]\sim S^{\prime }[j^{\prime }:]$ when the corresponding positions match, the resulting alignment graph is always acyclic.

A vertex-labeled tree with root $v_r\in V$ is a vertex-labeled graph $T=(V,E,\mathrm{\Sigma },\mathrm{\ell })$ such that the indegree of the root is $0$ and the indegree of every other vertex $v\in V$ is $1$, with a unique path $P_v$ from the root to every vertex $v\in V$. Given an edge $(u,v)\in E$, we call vertex $u$ the parent of vertex $v$ and $v$ a child of $u$. We call vertices with outdegree $0$ the leaves of the tree.

XBWT uses vertices in place of suffixes. The vertices are sorted by the label $\mathrm{\ell }(P_v)$ in colexicographic order (in lexicographic order by reverse path labels). For a given vertex $u\in V$, the BWT stores the concatenation of the labels $\mathrm{\ell }(v)$ for each outgoing edge $(u,v)\in E$. We encode the outdegree $d$ in unary as a one followed by $d-1$ zeros and concatenate the outdegrees in the same order to form bitvector $B_{out}$. In order to handle the leaves, we add a technical vertex $v_{\$}\notin V$ such that $\mathrm{\ell }(v_{\$})=\$$. For each leaf $v$, we encode an edge $(v,v_{\$})\notin E$. As with strings, we assume that $\mathrm{\ell }(v)\ne \$$ for every vertex $v\in V$. If the colexicographic rank of vertex $v\in V$ is $x$, the labels of its children are stored in $\mathrm{BWT}[B_{out}.\mathrm{select}(x,1):B_{out}.\mathrm{select}(x+1,1)-1]$. See Figure 1 for an example. With an uncompressed wavelet tree [29] for the BWT and an uncompressed bitvector for $B_{out}$, the XBWT data structure requires $\mathrm{O}\left(\left|V\right|\log \sigma \right)$ bits of space, not including the $\mathrm{C}$ array.

Let $\mathrm{find}(X)=[sp:ep]$ be the colexicographic range of vertices $v\in V$ with path labels $\mathrm{\ell }(P_v)$ ending with string $X$. Given a symbol $c\in \mathrm{\Sigma }$, we want to find the range $\mathrm{find}(Xc)$ of vertices with path labels ending with $Xc$. We first expand the range into a colexicographic range of outgoing edges $[s{p}_{out}:e{p}_{out}]=[B_{out}.\mathrm{select}(sp,1):B_{out}.\mathrm{select}(ep+1,1)-1]$. Every edge $(u,v)\in E$ in that range is encoded with the label $\mathrm{\ell }(v)$ of the child vertex, and the colexicographic rank of that vertex is based on path label $\mathrm{\ell }(P_v)=\mathrm{\ell }(P_u)\mathrm{\ell }(v)$. Hence we can use a backward search step (Section 2.4) to get the range of children with label $c$: $\mathrm{find}(Xc)=[\mathrm{LF}(s{p}_{out}-1,c)+1:\mathrm{LF}(e{p}_{out},c)]$ Note that the backward search step becomes a forward search step, as the XBWT is based on colexicographic order. Given a string $X$, we can therefore compute $\mathrm{find}(X)$ using $\mathrm{O}\left(\left|X\right|\right)$ rank operations on $\mathrm{BWT}$ and $\mathrm{O}\left(\left|X\right|\right)$ select operations on $B_{out}$.

If we know the colexicographic rank $x_r$ of the root $v_r$, we can match entire path labels instead of suffixes by starting the above search from range $[x_r:x_r]$. We can also use the XBWT for tree navigation. Given a vertex $v\in V$ with colexicographic rank $x$, the range of its children is $[B_{out}.\mathrm{select}(x,1):B_{out}.\mathrm{select}(x+1,1)-1]$. The colexicographic rank of its $i$-th child is $\mathrm{LF}(B_{out}.\mathrm{select}(x,1)+i-1)$ and the rank of its parent is $B_{out}.\mathrm{rank}(\mathrm{\Psi }(x),1)$.

The initial XBWT paper proposed using the structure as a fast and space-efficient representation for hierarchical documents, such as XML [21]. In such applications, it is important that a vertex can have multiple children with the same label (such as list items). When storing string collections (often called string dictionaries in this context), we can collapse the tree into a trie, improving space usage. Compared to other representations for string collections, the XBWT is smaller but slower [36]. The slowness comes from symbol-to-symbol navigation and the poor memory locality of iterated LF-mapping. However, the same FM-index functionality enables substring searches not supported by faster representations based on hashing or front coding.

The Aho–Corasick algorithm can find the occurrences of multiple patterns $𝒫$ in a string $S$ in linear time $\mathrm{O}\left(\left|S\right|+\Vert 𝒫\Vert +occ\right)$, where $occ$ is the number of occurrences. It is based on a trie $T=(V,E,\mathrm{\Sigma },\mathrm{\ell })$ of the patterns augmented with failure links (suffix links). A failure link from vertex $v\in V$ points to a vertex $u\in V$ such that $\mathrm{\ell }(P_u)$ is the longest proper suffix of $\mathrm{\ell }(P_v)$ among all vertices $u\in V$. The XBWT can be augmented with a $2\left|V\right|(1+\mathrm{o}\left(1\right))$-bit balanced parentheses representation to support failure links, making it a space-efficient Aho–Corasick automaton [35]. Here the relative slowness matters less, as the algorithm is based on symbol-to-symbol navigation anyway.

Because the XBWT represents a collection of strings, it can be built from other representations, such as the suffix array and the Burrows–Wheeler transform [35]. The link also works in the other direction: we can go from the XBWT directly to a run-length encoded multi-string BWT [12]. As a byproduct, we can efficiently find the order of the strings that minimizes the number of runs in the BWT [12]. Because the order of the strings determines how well the suffixes close to the end of the string can be run-length encoded, this can be valuable with highly repetitive collections of short strings.

In this section we define the Spectral Burrows-Wheeler transform, SBWT. We begin with few basic definitions:

padded $k$-spectrum adds a minimal set of $-padded dummy $k$-mers to ensure that in the de Bruijn graph (to be defined shortly), every non-dummy $k$-mer has an incoming path of length at least $k$:

The number of dummy vertices can be calculated with a simple formula. It is indeed equal to the number of internal vertices of $𝒯({\mathrm{R}}^{\prime })$, where $𝒯({\mathrm{R}}^{\prime })$ is the trie of the source vertex set of a set of strings.

For example, if we consider two strings $S_1=$ ACAGTG and $S_2=$ ATCAGA, and $k=3$, then $\mathrm{R}={\mathrm{S}}_3(S_1,S_2)=\{$ACA, CAG, AGT, GTG, ATC, TCA, AGA$\}$, ${\mathrm{R}}^{\prime }=\{$ACA, ATC$\}$, and the padded $k$-spectrum ${\mathrm{S}}_3^{+}(S_1,S_2)={\mathrm{S}}_3(S_1,S_2)\cup \{$$$$, $$A, $AC, $AT$\}$.

See figure 2 for an illustration of a de Bruijn graph.

We can now define the Spectral BWT.

Continuing the example above and adding a new sequence:
$𝒮=\{\text{ACAGTG},\text{ATCAGA},\text{TTGTCAGTGT}\}$. The colexicographically ordered list of ${\mathrm{S}}_3^{+}(𝒮)$ is: $$$, $$A, ACA, TCA, AGA, $AC, ATC, GTC, CAG, GTG, TTG, $$T, $AT, AGT, TGT, $TT, and the SBWT of $𝒮$ is the sequence of sets: {A,T}, {C,T}, {G}, $\mathrm{\varnothing }$, $\mathrm{\varnothing }$, {A}, {A}, $\mathrm{\varnothing }$, {A,T}, {T}, $\mathrm{\varnothing }$, {T}, {C}, {C,G}, $\mathrm{\varnothing }$, {G}.

The sets in the SBWT represent the labels of outgoing edges in the vertex-centric de Bruijn graph, such that we only include outgoing edges from $k$-mers that have a different suffix of length $k-1$ than the preceeding $k$-mer in the colexicographically sorted list. The padding of dollar-symbols is a technical detail required to make the SBWT work. It ensures that every $k$-mer of the $k$-spectrum has an incoming path of length at least $k$ in the de Bruijn graph. Alternatively, the sequences of the set $𝒮$ on which the SBWT is built could be made cyclic, thus removing the need for the padded $k$-spectrum ${\mathrm{S}}_k^{+}(𝒮)$.

We now describe how to implement efficient $k$-mer membership queries on the spectrum of the input strings, using only the information encoded in the SBWT. It is helpful to think of the SBWT as encoding the set ${\mathrm{R}}^{+}$ of $k$-mers as a list arranged in colexicographical order (as in the example above). In the context of $k$-mer lookup, the SBWT allows navigation of this ordered list via a search routine based on subset rank queries.

The search routine, introduced below, works by searching the $k$-mer character by character from left to right. It does so maintaining the interval of $k$-mers in the colexicographically-sorted list that are suffixed by the prefix which has been processed so far.

Since the $k$-mers in ${\mathrm{R}}^{+}$ are colexicographically sorted in the SBWT, the subset of $k$-mers in ${\mathrm{R}}^{+}$ that share a string $\alpha$ as a suffix are next to each other, and they are called the colexicographic interval of $\alpha$. The colexicographic interval of a string $\alpha$ longer than $k$ is defined as the range of $k$-mers in the sorted list of ${\mathrm{R}}^{+}$ suffixed by the last $k$ characters of $\alpha$. Now let ${[s:e]}_{\alpha }$ denote the colexicographic interval of the string $\alpha$, where $s$ and $e$ are respectively the colexicographic ranks of the smallest and largest $k$-mer in the SBWT having substring $\alpha$ as a suffix. Finding the interval $[s^{\prime }:e^{\prime }]$ of $\alpha c$, where $c$ is a character from $\mathrm{\Sigma }-\{\$\}$, is equivalent to following all edges labeled with $c$ from the vertices in $[s:e]$. Due to the way the edges are defined, the end points of these edges are a contiguous range $[s^{\prime }:e^{\prime }]$ such that $s^{\prime }$ is the destination of the first outgoing edge labeled $c$ from $[s:e]$, and $e^{\prime }$ is the destination of the last one. Let $C[c]$ be the number of edge labels with label smaller than $c$ in the graph. We now have the following formulas:

where $subsetran{k}_c$ is a subset rank query on the SBWT sequence. The “+1” at the start of the formulas is to skip over the vertex of ${\text{ \$}}^k$. The values $C[c]$ can be precomputed for all characters $c\in \mathrm{\Sigma }$. By iterating these formulas $k$ times, we have the $k$-mer search routine.

The simplest data structures for subset rank queries is the plain matrix. See Figure 3 for an example.

The rows of the matrix $M$ are indexed for constant time rank queries [37] (implementation from [28]). A single rank query on row $i$ at index $j$ in $M$ would answer in constant time $subsetran{k}_c(j)$, where $c$ is the $i^{th}$ character of the alphabet. There exist alternative succinct data structures for subset rank, such as the subset wavelet tree, details of which can be found in [4, 5]. A concatenation of the sets indexed as a wavelet tree leads to something resembling the BOSS data structure. All $k$-mers of a longer string can be searched efficiently by adding the $\mathrm{LCS}$ array of the SBWT [3]. At the end of each single $k$-mer seach, the $\mathrm{LCS}$ array allows to extend the colexicographic interval to the last $k-1$ characters observed in the previous $k$-mer.

Consider a collection that consists of $m$ identical copies of the same string $S$. Then for any position $j$, the suffixes starting at that position are consecutive in lexicographic order. Because the preceding symbol $S[j-1]$ (if $j>1$; otherwise $S[n]=\$$) is also the same for every string, the suffixes are in the same run of symbol $S[j-1]$ in the BWT of the collection. Hence the BWT of a collection of $m$ identical strings always has the same number of runs, regardless of the number of copies. We now consider what happens when we apply edits to a collection of identical strings.

Let $𝒮=(S_1,\mathrm{\dots },S_m)$ be a highly repetitive collection of similar strings. We assume that there is a sequence of $s$ edit operations that transforms a collection of $m$ copies of string $S$ into $𝒮$, with $s\ll mn$. We further assume that there is a multiple sequence alignment corresponding to the edit operations, where suffixes that have been derived from the same suffix of the original string $S$ start with a match or a mismatch in the same column. Note that different sequences of edit operations leading to the same collection may correspond to different alignments.

Initially, there are $n$ structural runs, one for each position in string $S$. Each edit operation that changes a substring $S_i[j:j^{\prime }]$ can affect the suffixes of $S_i$ in three ways:

1. 1.

   The symbol in the BWT can change for the suffixes starting from $S_i[j+1]$ to $S_i[j^{\prime }+1]$.

2. 2.

   Suffixes starting from $S_i[j-k]$ to $S_i[j^{\prime }]$, for some $k\ge 0$, may leave their current structural runs due to changes in lexicographic ranks.

3. 3.

   New unaligned suffixes can be created and existing suffixes may be removed if substring $S_i[j:j^{\prime }]$ is replaced with a string of different length.

Every affected suffix can increase the number of structural runs by at most two: by becoming a new run and by splitting an existing run in two.

Mäkinen et al. [38] showed that if the original string $S$ is random and we apply random edits, $k$ is at most $\mathrm{O}\left({\log }_{\sigma }N\right)$ in the expected case. If we assume single-symbol edits, $s$ edit operations increase the number of structural runs by at most $\mathrm{O}\left(s{\log }_{\sigma }N\right)$ in the expected case. Because structural runs are (possibly non-maximal) BWT runs, there are at most $r+\mathrm{O}\left(s{\log }_{\sigma }N\right)$ runs in the BWT of collection $𝒮$ in the expected case. Here $r$ is the number of runs in the BWT of string $S$.

Let $𝒮$ be a string collection and ${\sim }_0^b$ an equivalence relation that denotes that two suffixes are in the same structural BWT run. Let $G^{{\sim }_0^b}=(V,E,\mathrm{\Sigma },\mathrm{\ell })$ be the alignment graph of the collection with respect to ${\sim }_0^b$.

Graphs like this have been used as space-efficient representations of the collection. Both SimDNA [33] and FM-index of alignment [40, 39] index the collection using a graph based on a specific type of alignment. The alignment alternates between common regions, where all strings are aligned, and non-common regions, where they diverge. BWT tunneling [7] also uses a similar graph to improve BWT-based data compression by not storing run lengths within unary paths. Our goal is instead to extend the FM-index of the collection to index all strings in the closure ${𝒮}^{{\sim }_0^b}$ using the graph.

We could augment the FM-index with a bitvector $B_r[1:N]$ that marks the start (or equivalently the end) of each structural run with ones. After each backward search step, we use some rank/select operations on $B_r$ to extend the interval $[sp:ep]$ to the smallest interval $[s{p}^{\prime }:e{p}^{\prime }]$ that contains it and consists of entire structural runs. If we continue the search from interval $[s{p}^{\prime }:e{p}^{\prime }]$, we are effectively accepting recombinations within each vertex in the range. No matter which suffixes we used to reach a vertex, we can continue the search to the predecessor of any suffix in the vertex. See Figure 4 for an example.

Because we discard the information on how we arrived to the current vertex, we no longer need to keep track of the exact BWT interval. Hence we choose the GCSA representation [48] (Figure 5) that corresponds more closely to the structure of the graph. We keep the vertices in the same order as the structural runs. Let ${\mathrm{rank}}_V:V\to [1:|V|]$ be a function that maps a vertex $v\in V$ to its lexicographic rank ${\mathrm{rank}}_V(v)$. Lexicographic range $[sp:ep]$ now corresponds to vertices $v\in V$ such that ${\mathrm{rank}}_V(v)\in [sp:ep]$.

The data structure consists of three components $B_{in}$, $\mathrm{BWT}$, and $B_{out}$. We use them all for computing a backward search step from $\mathrm{find}(X)=[sp:ep]$ to $\mathrm{find}(cX)=[s{p}^{\prime }:e{p}^{\prime }]$ for a string $X$ and a symbol $c$.

1. 1.

   Bitvector $B_{in}$ encodes the indegree of each vertex in unary. Similar to using $B_{out}$ in XBWT (Section 3.1), we map the range of vertices to a range of incoming edges as $[s{p}_{in}:e{p}_{in}]=[B_{in}.\mathrm{select}(sp,1):B_{in}.\mathrm{select}(ep+1,1)-1]$.

2. 2.

   The BWT range for vertex $v\in V$ is the concatenation of $\mathrm{\ell }(u)$ for each vertex $u\in V$ with $(u,v)\in E$ in an arbitrary order. A backward search step with this BWT maps the range $[s{p}_{in}:e{p}_{in}]$ of incoming edges to a range $[s{p}_{out}:e{p}_{out}]=[\mathrm{LF}(s{p}_{in}-1,c)+1:\mathrm{LF}(e{p}_{in},c)]$ of outgoing edges, ordered by the predecessor vertex.

3. 3.

   Bitvector $B_{out}$ encodes the outdegree of each vertex in unary. With rank queries, we get the interval of the predecessors of the vertices in $[sp:ep]$ with label $c$ as $[s{p}^{\prime }:e{p}^{\prime }]=[B_{out}.\mathrm{rank}(s{p}_{out},1):B_{out}.\mathrm{rank}(e{p}_{out},1)]$.

Together these take $\left|E\right|(\log \sigma +2)(1+\mathrm{o}\left(1\right))$ bits, assuming an uncompressed wavelet tree for the BWT and uncompressed bitvectors and excluding the $\mathrm{C}$ array. To handle the edge cases with source/sink vertices, we add technical edges connecting vertices with outdegree $0$ to vertices with indegree $0$.

This data structure also facilitates graph navigation. For instance, if $v\in V$ is a vertex and $x={\mathrm{rank}}_V(v)$, we can navigate backward along a unary path from $v$. This is effectively the same as using $\mathrm{LF}(x)$ to navigate from a suffix to its predecessor. The rank of the first incoming edge to $v$ is $x_{in}=B_{in}.\mathrm{select}(x,1)$. If $B_{in}[x_{in}+1]=0$, vertex $v$ has multiple predecessors and the path is no longer unary. With $x_{out}=\mathrm{LF}(x_{in})$, we get the rank of the corresponding outgoing edge, and $B_{out}.\mathrm{rank}(x_{out},1)$ is the rank of the predecessor vertex. Navigation like this is used in locate queries (Section 5.5).

The alignment graph based on structural runs allowed us to index the original strings as well as their recombinations within the runs. Our ultimate goal is allowing recombinations at any aligned position. In principle, we already have the necessary tools for it. Because ${\sim }_0^b$ is a refinement of ${\sim }_0$, we have ${𝒮}^{{\sim }_0}={({𝒮}^{{\sim }_0})}^{{\sim }_0^b}$. In other words, recombinations at aligned positions include recombinations at aligned positions within the same BWT run. We could therefore build the GCSA structure for the closure ${𝒮}^{{\sim }_0}$. The challenge is doing that without representing the closure explicitly.

Before proceeding, we relax the definition of structural runs. The GCSA structure only requires that the equivalence classes of suffixes are lexicographic ranges that start with the same symbol. By not requiring that they also share the preceding symbol, we simplify the discussion and make the alignment graph $G^{{\sim }_0^s}$ smaller.

The algorithm we use for building GCSA [48, 49] is similar to the prefix-doubling approach to suffix array construction. We start from graph $G_1=G^{{\sim }_0}=(V,E,\mathrm{\Sigma },\mathrm{\ell })$ and build a series of graphs $G_l=(V_l,E_l,\mathrm{\Sigma },{\mathrm{\ell }}_l)$ for increasing path lengths $l$. If we consider the graphs finite automata, they all recognize the same language ${𝒮}^{{\sim }_0}$. There are two types of vertices $v\in V_l$:

1. 1.

   The vertex corresponds to a path $P_v$ of length $\left|P_v\right|=l$ starting from vertex $\mathrm{locate}(v)=P_v[1]\in V$ in graph $G^{{\sim }_0}$. For every maximal path $P$ starting from vertex $v$ in graph $G_l$, we require that $\mathrm{\ell }(P_v)$ is a prefix of ${\mathrm{\ell }}_l(P)$.

2. 2.

   The vertex corresponds to a half-open lexicographic range $[X_v:Y_v)$ of path labels. There is an arbitrary path $P_v$ of length $\left|P_v\right|\le l$ starting from vertex $\mathrm{locate}(v)=P_v[1]\in V$ in graph $G^{{\sim }_0}$ such that $\mathrm{\ell }(P_v)\in [X_v:Y_v)$. For every maximal path $P$ starting from vertex $v$ in graph $G_l$, we require that ${\mathrm{\ell }}_l(P)\in [X_v:Y_v)$. For every other vertex $u\in V_l$ with $u\ne v$, we require that $\mathrm{\ell }(P_u)\notin [X_v:Y_v)$ (if $u$ is of type 1) or $[X_u:Y_u)\cap [X_v:Y_v)=\mathrm{\varnothing }$ (if $u$ is of type 2).

In both cases, we set ${\mathrm{\ell }}_l(v)=\mathrm{\ell }(\mathrm{locate}(v))$.

If all vertices $v\in V_l$ are of type 2, we have a graph $G_l$ that is equivalent to $G^{{\sim }_0}$, with vertices that can be sorted unambiguously by the labels of maximal paths starting from them. For each suffix $S_i[j:]$ of the closure ${𝒮}^{{\sim }_0}$, there is a vertex $v\in V_l$ such that $S_i[j:]\in [X_v:Y_v)$ and a path $P$ starting from $v$ with label ${\mathrm{\ell }}_l(P)=S_i[j:]$. Therefore the vertices of $G_l$ correspond to (possibly non-maximal) structural runs. If we merge them to form maximal runs, we get the graph $G^{{\sim }_0^s}$. See Figure 5 for an example.

We derive graph $G_{2l}$ from graph $G_l$ with a doubling step followed by a pruning step. Let $u\in V_l$ be a type 1 vertex, and let ${\mathrm{rank}}_l(u)$ be the lexicographic rank of $\mathrm{\ell }(P_u)$ in the set of distinct path labels $\{\mathrm{\ell }(P_v)\mid v\in V_l\}$. In the doubling step, we combine it with all vertices $v\in V_l$ such that $P_u{P}_v$ is a valid path. For every combination, we create a new vertex $w$ with $P_w=P_u{P}_v$ and ${\mathrm{rank}}_{2l}(w)=({\mathrm{rank}}_l(u),{\mathrm{rank}}_l(v))$. The new vertex will be of the same type as vertex $v$. This requires that graph $G^{{\sim }_0}$ is reverse deterministic; see below for a discussion. If $u\in V_l$ is a type 2 vertex, we copy it as vertex $w\in V_{2l}$ with ${\mathrm{rank}}_{2l}(w)=({\mathrm{rank}}_l(u),0)$. The doubling step finishes by sorting the new vertices by their ranks and replacing the pairs of ranks with integer values.

In the pruning step, we merge lexicographic ranges of type 1 vertices into type 2 vertices. Let $V_{2l}^x=\{v\in V_{2l}\mid {\mathrm{rank}}_{2l}(v)=x\}$ be the set of vertices with rank $x$. If all vertices $v\in V_{2l}^x$ share the same $\mathrm{locate}(v)$ vertex, the shared path label uniquely defines the starting vertex. We can therefore merge them into a type 2 vertex. Similarly, if sets $V_{2l}^x$ and $V_{2l}^{x+1}$ also share the starting vertex, we can merge them into a single vertex.

If $S_{\max }$ is the longest string in ${𝒮}^{{\sim }_0}$, the prefix-doubling algorithm finishes after at most $\lceil \log \left|S_{\max }\right|\rceil$ doubling steps. In the worst case, the final graph $G^{{\sim }_0^s}$ can be exponentially larger than the graph $G^{{\sim }_0}$ we started from. However, if the distances between the positions of the original string $S$ with edits in any string $S^{\prime }\in 𝒮$ are sufficiently high, graph $G^{{\sim }_0^s}$ will have $\mathrm{O}\left(\left|V\right|\right)$ vertices and $\mathrm{O}\left(\left|E\right|\right)$ edges in the expected case [49]. If we assume a random string and random edits, the phase transition to exponential growth occurs at average distance $\mathrm{O}\left({\log }_{\sigma }\left|S\right|\right)$.

We do not need to maintain the edges of the intermediate graphs. Determining them for the final graph $G^{{\sim }_0^s}=(V^{\prime },E^{\prime },\mathrm{\Sigma },{\mathrm{\ell }}^{\prime })$ is simple enough. Given a vertex $v^{\prime }\in V^{\prime }$ of the final graph and an edge $(u,\mathrm{locate}(v^{\prime }))\in E$ of the original graph, we know that there must be a corresponding edge $(u^{\prime },v^{\prime })\in E^{\prime }$ such that $\mathrm{locate}(u^{\prime })=u$ in the final graph. We can therefore generate the edges as $(u,v^{\prime })$ and sort them by $(\mathrm{\ell }(u),{\mathrm{\ell }}^{\prime }(P_{v^{\prime }}))$. This also sorts them by ${\mathrm{\ell }}^{\prime }(P_{u^{\prime }})$. By scanning the lists of vertices and edges in sorted order, we can then match edges $(u,v^{\prime })$ to source vertices $u^{\prime }$ by $\mathrm{locate}(u^{\prime })=u$.

This construction algorithm only works if the original graph $G^{{\sim }_0}$ is reverse deterministic. For every vertex $v\in V$, we require that each predecessor $u\in V$ with $(u,v)\in E$ has a distinct label $\mathrm{\ell }(u)$. If two predecessors $u,u^{\prime }\in V$ share a label, $\mathrm{\ell }(uP)=\mathrm{\ell }(u^{\prime }P)$ for every path $P$ starting from vertex $v$. The pruning step cannot merge them, because no extension of the path label can determine the starting vertex uniquely. On the other hand, if the graph is reverse deterministic and $v\in V$ is a type 2 vertex, $\mathrm{\ell }(u{P}_v)$ always determines the unique predecessor $u$.

Determinizing a finite automaton takes exponential time in the worst case, because the states of the determinized automaton correspond to sets of states of the non-deterministic automaton. However, when the textbook algorithm was used to determinize the graph before GCSA construction, it was observed to be fast [49]. This was later explained using Wheeler graphs [16]. The exponent in the time complexity is not the total number of states but the maximal number of states that are not comparable in lexicographic order. Intuitively, the closer the graph is to a Wheeler graph, the faster the determinization is. If the graph is simple enough that the prefix-doubling algorithm only increases its size by a constant factor, it must already be close to being a Wheeler graph.

Let $G^{{\sim }_0}=(V,E,\mathrm{\Sigma },\mathrm{\ell })$ be the original alignment graph and $G^{{\sim }_0^s}=(V^{\prime },E^{\prime },\mathrm{\Sigma },{\mathrm{\ell }}^{\prime })$ the final graph in GCSA construction. During the construction, we annotated each vertex $v\in V_l$ in each intermediate graph $G_l=(V_l,E_l,\mathrm{\Sigma },{\mathrm{\ell }}_l)$ (including $G^{{\sim }_0}$ and $G^{{\sim }_0^s}$) with the corresponding vertex $\mathrm{locate}(v)\in V$ of the original graph. As the notation implies, we want GCSA locate queries return positions in the original graph. While it would be easier to return positions in $G^{{\sim }_0^s}$ (we built GCSA for that graph, after all), it is only a technical construct. Positions in the original graph (or columns in the alignment) are more likely to be meaningful to the user.

A locate query in an FM-index relies on the fact that $\mathrm{SA}[x]=\mathrm{SA}[\mathrm{LF}(x)]+1$. If a position has not been sampled, we can derive the value from the predecessor. Let $\mathrm{id}:V\to [1:|V|]$ be a function that gives each vertex $v\in V$ a unique integer identifier. To apply the same idea with GCSA, we require that $\mathrm{id}(v)=\mathrm{id}(u)+1$ for each edge $(u,v)\in E$ such that $u$ is the only predecessor of $v$ and $v$ is the only successor of $u$. If we sample $\mathrm{id}(\mathrm{locate}(v))$ for each vertex $v\in V^{\prime }$ that has multiple predecessors or $\mathrm{id}(\mathrm{locate}(v))\ne \mathrm{id}(\mathrm{locate}(u))+1$ for the only predecessor $u$, the locate algorithm works as before [48, 49]. (We also want to sample $\mathrm{id}(\mathrm{locate}(v))$ at regular intervals on unary paths to ensure query performance.)

This approach was originally designed for sampling the LCP array [45]. If $\mathrm{BWT}[x]$ is not at the start of a run, $\mathrm{LCP}[x]=\mathrm{LCP}[\mathrm{LF}(x)]-1$. We therefore need to sample the irreducible values, while the rest can be derived from the predecessor. More generally, the approach can be understood as a sampled tag array [8]. We annotate each position in the FM-index with a tag, which can store arbitrary information related to the position. If most tags can be derived easily from the only predecessor, we can store the tag array space-efficiently and still support efficient queries.

As mentioned earlier, there is a phase transition in GCSA construction. As the distances between positions with edits in the alignment get below a critical threshold, the number of vertices and edges in the final graph $G^{{\sim }_0^s}$ starts growing rapidly. This makes the GCSA impossible to build in practice for many graphs. One way to handle this is adding a context length requirement to the definition of alignment [48, 49].