Pangenome Graph Indexing via the Multidollar-BWT

Let $\mathrm{\Sigma }=\{{\text{c}}_1,\mathrm{\dots },{\text{c}}_{\sigma }\}$ be a finite, ordered alphabet of cardinality $\sigma =|\mathrm{\Sigma }|$, and let $\text{\$},{\text{\$}}_1,{\text{\$}}_2,\mathrm{\dots }$ be an infinite number of terminal symbols. Let $\prec$ denote the lexicographic order; we assume that $\text{\$}\prec {\text{\$}}_1\prec {\text{\$}}_2\prec \mathrm{\cdots }\prec {\text{c}}_1\prec \mathrm{\dots }\prec {\text{c}}_{|\mathrm{\Sigma }|}$. We define a string $S$ over $\mathrm{\Sigma }$ to be a concatenation of $n$ symbols $S=S[1..n]$. We denote the $i$-th prefix of $S$ as $S[1..i]$, the $i$-th suffix as $S[i..n]$, and the substring spanning positions $i$ through $j$ as $S[i..j]$. The set of all non-empty strings made of symbols from $\mathrm{\Sigma }$ is denoted by ${\mathrm{\Sigma }}^{+}$. If $P$ and $T$ are two strings, we say that $P$ occurs in $T$ at position $1\le i\le |T|$ (or, equivalently, that $i$ is an occurrence of $P$ in $T$) if $T[i..i+|P|-1]=P$.

Given an input string $S$ ending with symbol $\text{\$}\notin \mathrm{\Sigma }$, we consider all rotations of $S$ in lexicographical order. We denote the matrix of these rotations, which is commonly referred to as the Burrows-Wheeler Transform Matrix [3], as $ℳ$. Since the rows of $ℳ$ are stored lexicographically all occurrences of any pattern $P$ occur together in a contiguous range, implying that finding all occurrences of $P$ in $S$ correspond to finding the contiguous range of rows beginning with $P$. Although $ℳ$ requires $𝒪(n^2)$-space to be stored, it is sufficient to store and use only the first column (denoted as F) and the last column (denoted as L) of $ℳ$: this crucial property is known in the literature as the LF-mapping [8], which maps the $i$-th occurrence of a symbol c in L to the $i$-th occurrence of a symbol c in F. In fact, such two occurrences of c in columns $F$ and $L$, correspond to the same symbol in the original text $T$.

The Burrows-Wheeler Transform of $S$, denoted as ${\text{BWT}}_S$, is in fact column L. We denote ${\text{BWT}}_S$ as simply BWT if there is no ambiguity. The starting positions of all sorted rotations is a permutation of $\{1,\mathrm{\dots },n\}$, and is referred to as the suffix array [17] of $S$ since it is equivalent to the lexicographical sorting of all the suffixes of $S$. We denote this by ${\text{SA}}_S$, or just SA if $S$ is clear from the context. It is not difficult to see that $\text{BWT}[i]$ is the symbol preceding $\text{F}[i]$ in $S$ and thus BWT can be computed from the suffix array: $\text{BWT}[i]=S[\text{SA}[i]-1]$, where, due to the rotational nature of the transform, we set $\text{BWT}[i]=\text{\$}$ in case of $\text{SA}[i]=1$.

A classic problem in computer science is the string matching problem: given strings $T$ (the text) and $P$ (the pattern), the goal is to count and/or locate all occurrences of $P$ in $T$. Arguably one of the most successful ideas to solve this problem, especially in the presence of large texts, is the celebrated FM-index [8] which heavily relies on the LF-mapping. Consider the interval $[i..j]$ of rows of $ℳ$ that start with $P$. As already mentioned above, ${\text{BWT}}_T[i..j]$ corresponds to the symbols preceding $\text{F}[i..j]$ in $T$, hence, for any $\text{c}\in \mathrm{\Sigma }$, the interval $[i^{\prime }..j^{\prime }]$ of rows that start with $\text{c}\cdot P$ can be retrieved via the LF-mapping. This operation, which is referred to c-backward step or c-backward extension (and for which we will omit the c-prefix when there is no ambiguity), can be implemented to run in either constant time [8] or, in case one would like to lessen the memory burden, in $𝒪(\log \sigma )$ using suitable data-structures [10]. In the following, any such interval $[i..j]$ in BWT will also be called c-backward interval, and the complete sequence of backward steps to retrieve the entire $P$ as the backward search.

There has been extensive work in the literature to generalize the Burrows-Wheeler Transform from single strings to collections of strings. We will use one such generalization, namely the Multidollar Burrows-Wheeler Transform [18, 4].

The following two simple properties, which are direct consequences of Definition 1 and which are discussed in Example 9, will be used in our algorithm.

In this paper, we are concerned with a more general version of the pattern matching problem where a pattern is searched inside a labelled graph.

In the following, whenever we mention the term graph we are assuming that it is node-labelled. We extend the labelling function $\mathrm{\ell }(\cdot )$ to any path $\pi =(v_1,\mathrm{\dots },v_k)$ by simply representing $\mathrm{\ell }(\pi )=\mathrm{\ell }(v_1)\mathrm{\cdots }\mathrm{\ell }(v_k)$.

In this context, in order to locate where a pattern $P$ occurs in the graph $G$ we define two notions: (i) the coarse occurrence, in which we are interested only in the starting node (say $v$) of a path where $P$ occurs; (ii) the fine occurrence, where we are also interested in the full path where $P$ occurs. Formally:

Notice that, in contrast to the string-against-string matching, both a fine and a coarse occurrence can be testified by more than one path.

The goal of this work is to solve the following:

GCSA (Generalized Compressed Suffix Array) [23] and GCSA2 [21, 22] are BWT and FM-index generalizations developed initially for directed acyclic graphs (GCSA) and later extended to support also cyclic graphs (GCSA2). These approaches were inspired by similar BWT-based indexes for de Bruijn graphs [2]. The latter is currently integrated as a part of the VG toolkit [9].

Previous experimental evaluations proposed in GCSA2’s publication proved the efficiency of this index in terms of performance; according to the authors, GCSA2 performs pattern matching on graph path label concatenations almost as well as classical FM-index approaches on texts of similar sizes [21].

In this paper, we propose an alternative graph indexing method with the goal of overcoming two of the main limitations of GCSA2: (i) the high memory requirements for indexing the graph, especially in the presence of complex regions and (ii) the possibility of obtaining false positive results for long queries. Regarding (i), this behavior is strongly related to the GCSA2 construction algorithm which relies on various iterations of the prefix-doubling algorithm to build the underlying de-Bruijn graph. This approach consists of repetitively joining $k$-length paths into paths of length $2k$ in the input graph; even with two or three iterations, on large graphs, the algorithm requires hundreds of gigabytes of RAM even if it stores all the possible temporary data on disk, as stated in [21]. With respect to issue (ii), the length of the queries is bounded by $k$, where $k$ is the k-mer size used in the underlying de-Bruijn graph, currently the maximum size possible is $k=256$, due to technical details in the implementation¹¹1https://github.com/jltsiren/gcsa2.

Consider now one step in Algorithm 1. The complexity of this step scale in the size of the backward-intervals set, i.e. we have to perform one backward-extension for each interval. If we have merged all the consecutive intervals, solving coarse Problem 6, in the worst case we need to consider $|\text{BWT}|/2$ for the next step. On the other hand, if we consider the backward-intervals multiset in the case of fine occurrences we can not bind its cardinality, making it difficult to establish a time and space complexity bound for finding all fine occurrences.

Focusing on the first case, having to iterate the pattern $P$ and having to perform at most $|\text{BWT}|/2$ extension for each symbol of $P$, we have in total $|P|(|\text{BWT}|/2)$ backward steps. In addition, the time required to retrieve the ending nodes in each backward-interval and their predecessors is bounded by the number of possible backward intervals multiplied by the number of nodes, i.e. $(|\text{BWT}|/2)\cdot |V|$. Having $|\text{BWT}|\gg |V|$, we can conclude that our algorithm runs in $𝒪(|P||\text{BWT}|B)$ time, where $B$ is the time of a single backward step. Regarding space, we are bounded by the size of the BWT and the size of the graph. We need to store the multidollar-BWT index, which requires $𝒪(|\text{BWT}|)$ space, the graph adjacent list, which requires $𝒪(|V|+|E|)$ space, and the set of backward-intervals, which requires $𝒪(|\text{BWT}|)$ space.

Recalling how we designed our approach, after the first iteration we will add to the backward-intervals set one interval for each predecessor of every node that starts with the last pattern symbol. In fact, without considering the merging step or any further optimization, we have to check all the dollar symbols in the backward-interval we get after the first extension. In Example 9, suppose that the pattern ends with the symbol $C$. With the first backward-extension, we fall into the interval $[9,17)$, which contains three dollars: ${\$}_3$, ${\$}_5$ and ${\$}_4$. By Property 7, $C$ is the first symbol of labels $L_4$, $L_5$ and $L_1$ and we need to add one interval for each predecessor of these three nodes.

We can expect the cardinality of this new set of backward-intervals to be $\approx |V|/|\mathrm{\Sigma }|$, assuming a uniform distribution of the first symbols on the labels $L$. Iteration by iteration the size of the set decreases exponentially, having discarded all the mismatching paths. In Figure 1 we show this behaviour. One of the main bottlenecks of this algorithm is the extension of the backward-intervals set in the first iterations, as the cardinality of the set is exponentially larger when we are searching short suffixes of the pattern. To alleviate this issue and allow to scale on large pangenome graphs, such as the ones from HPRC, we implemented a preprocessing step that computes the $ℒ$ sets for any possible $x$-length string over the alphabet $\mathrm{\Sigma }$, for fixed $x$ during the indexing phase. In detail, we precompute all the matches for all ${|\mathrm{\Sigma }|}^x$ strings. We refer to this precomputed set as cache of length $x$. This cache allows us to skip all the first $x$ iterations during the querying phase and access directly the backward-interval sets for $P[|P|-x]$, thus skipping the exponential search. Therefore, trade-off the initial computation time and disk usage to gain a considerable speedup in time during the querying. The entire cache needs to be stored on disk and it needs to be entirely loaded into memory when we query the index. Increasing the value of $x$ allows us to skip more iterations, but at the cost of disk and RAM usage, therefore we utilized the results from Figure 1 to select a threshold that gives us a compromise.

By computing the cache the indexing time is increased, thus to optimize this task, we developed a parallel algorithm that starts from the backward-intervals of each $\sigma \in \mathrm{\Sigma }$ which is then extended with each ${\sigma }^{\prime }\in \mathrm{\Sigma }$; such extensions can be computed in parallel. The process is repeated $x$ times until the matches of all the $x$-length strings are computed. With a naïve approach we would have many repeated extensions to compute that are not necessary; for example two patterns $CAAA$ and $TAAA$ all have the common $AAA$ suffix, that does not require to be recomputed for each pattern, since the backward-interval set is the same for all common suffixes. Therefore we implemented a tree-based approach in which each iteration extends a suffix with the all the $\sigma$ without the need to recompute already processed patterns. With this approach, we perform only the necessary ${\sum }_{i=1}^x{|\mathrm{\Sigma }|}^i$ extensions, instead of the all the $x\cdot {|\mathrm{\Sigma }|}^x$; furthermore, each of these extensions can be executed in parallel, since they are independent. Thanks to this optimization we can heavily reduce the impact of the cache computation time.

Figure 1 clearly shows the trend of the cardinality of the backward-interval sets $ℒ$ - further motivating the use of the caching step - that exponentially decreases at each iteration, reaching a small enough value at around step 5. The vertical lines in the plot indicate the average position where the cardinalities reach a value of 1. In the experimentation, we tested various cache lengths to evaluate its impact, suggesting to use a value of around 7 to obtain a good trade-off between computational time, memory usage and querying speed. Increasing the length too much would also result in computing patterns that would rarely occur when performing pattern matching (due to their higher granularity) that results in wasted resources.

While technically the cardinality of $ℒ$ is not strictly monotonic decreasing, as in the presence of multiple predecessors with the same ending symbol dimension may increase for a few iterations, it is also likely to drop shortly after or to stabilize to the number of actual occurrences of the pattern.

In our experiments, we tested three different pangenome graphs:

• $\mathrm{■}$

  Drosophila pangenome graph with 10M nodes, 12M edges, a maximum max degree of 15, an average degree of 1.18 and 144M bp. It was built by VG using 205 haplotypes [15]. This graph has at most 32 bp per node and has been pruned to delete complex regions with many variants that prevent indexing using GCSA2³³3https://github.com/vgteam/vg/issues/4404;

• $\mathrm{■}$

  Human Pangenome Reference Consortium (HPRC) pangenome graph [14] built using MiniGraph [13]. This graph has 12M nodes, 13M edges, a maximum degree of 74, an average degree of 1.02 and 3.2B bp. We slightly modified this graph to have all nodes with at most 256 bp, trying to avoid issues with GCSA2 indexing⁴⁴4https://github.com/jltsiren/gcsa2/issues/44 and fixing the labels’ maximal length to the maximum k-mers size supported by the tool (256bp);

• $\mathrm{■}$

  HPRC pangenome graph built using MiniGraph-Cactus [11]. This graph has 80M nodes, 110M edges, a maximum degree of 39, an average degree of 1.39, and 3.1B bp. No modification has been necessary to use this graph.

Table 5 shows the file sizes on disk of all gindex components for the three input graphs. In detail, we need to store the RopeBwt2 index, the topology of the graph as the adjacency list of the predecessors, the label names and the BWT dollars order. The latter is needed because RopeBwt2 does not distinguish different dollars.

To test the efficacy of gindex we tested the method on all the pangenome graphs with three different modes of execution:

• $\mathrm{■}$

  Cache length $x$: indexing with a cache computation of length $x$ and returning the starting nodes of each match (coarse occurrence);

• $\mathrm{■}$

  No Cache: indexing without using any cache and returning the starting nodes of each match (coarse occurrence);

• $\mathrm{■}$

  Full Path: indexing without using any cache and returning the complete paths that match the query (fine occurrence).

These settings would allow us to explore in detail the properties of both the method and its implementation.

Note that we do not perform a full path match with caching mode because this would require to include in the cache also all the possible paths in the precomputed intervals, therefore greatly affecting memory usage, computation time and implementation complexity without a clear advantage, considering the fact that the only other available implementation (to the best of our knowledge) GCSA2 only reports the starting node of the match.