Fast Pseudoalignment Queries on Compressed Colored de Bruijn Graphs

We study the problem of accelerating pseudoalignment queries directly over the compressed, repetition-aware representations developed for the Fulgor index. More precisely, we devise new query algorithms that exploit the repetitiveness of the color sets to actually perform less work, for both intersection and union of color sets. Our experimental results across large pangenomes show that such repetition-aware pseudoalignment queries can be generally performed $2-3\times$ faster. To make a concrete example, the smallest Fulgor index mentioned above now processes the same workload ($6.6\times {10}^6$ pseudoalignment queries) in 34 minutes, rather than the previous 1 hour and 50 minutes, and without consuming more space than 6.3 GB of RAM.

This representation compresses a color set using one of three different strategies, each tailored to the density of the set. Among the three density classes – sparse, dense, and very dense – the last one presents the most significant opportunities for efficiency gains. Remember that very dense color sets are compressed by storing their complement set, denoted with $\overline{C}$, which is sparse by definition. So, rather than iterating through a large set, we can directly leverage the complement set for faster processing: given $z$ very dense sets, represented as ${\overline{C}}_1,\mathrm{\dots },{\overline{C}}_z$, we can use de Morgan’s laws to compute the intersection of their original sets, that is ${\bigcap }_{i=1}^z{C}_i=\overline{\left({\bigcup }_{i=1}^z{\overline{C}}_i\right)}.$ In other words, the intersection of $z$ sets is equivalent to the complement of the union of the complement sets. The refined intersection algorithm is therefore:

1. 1.

   Compute the union of the complements of all very dense color sets.

2. 2.

   Compute the intersection of all other color sets.

3. 3.

   Return the intersection between the complement of the set computed at step 1. and the set computed at step 2.

We recall that in this variant the color sets are stored as the symmetric differences $D_i=C_i\mathrm{\Delta }{A}_j$ between the original color set $C_i$ and the representative set $A_j$ of the group ${𝒩}_j$, such that $C_i\in {𝒩}_j$ and $|D_i|\ll |A_j|$. Note that simply using Algorithm 1 over this representation is inefficient. In fact, to decode $C_i$, it is necessary to read its whole differential set $D_i$ and representative set $A_j$, for which we have that $|D_i|+|A_j|\ge |C_i|$. Furthermore, Algorithm 1 does not exploit the fact that multiple color sets can be part of the same group ${𝒩}_j$, and thus have the same representative set $A_j$. Consequently, each $A_j$ could be processed redundantly, leading to a significant performance bottleneck.

Consider the case where all color sets belong to the same group ${𝒩}_j$. For ease of notation, throughout the rest of this section, we will simply write $A$ and $𝒩$, instead of $A_j$ and ${𝒩}_j$. Then, the time complexity of the intersection is $\text{𝒪}\left(z|A|+{\sum }_{i=1}^z(|D_i|+|C_i|)\right).$ This complexity is not ideal because we would like to scan each $D_i$ only rather than $D_i$ plus $C_i$, and $A$ only once rather than $z$ times. To do so, we give the following lemma (the proof can be found in Appendix A).

In other words, Lemma 2 says that a color is part of the intersection if it appears only in the representative set or, exclusively, in all differential sets. Note that this alternative formulation does not involve the original $C_i$ but only their symmetric differences $D_i$ and the representative set $A$.

The intersection of $z$ differential color sets, belonging to the same group, can be implemented as in Algorithm 2. The algorithm has a time complexity of $\mathrm{\Theta }\left(N+|A|+{\sum }_{i=1}^z|D_i|\right)$, since we have to read only once the representative and differential sets, plus the array $counts[1..N]$. For $c=1,\mathrm{\dots },N$, we define $counts[c]:=|\{1\le i\le z:c\in D_i\}|$, i.e., the number of differential sets that include $c$. We use this array to compute the intersection between the $A$ and every ${\overline{D}}_i$, and the intersection between $\overline{A}$ and every $D_i$, as follows. For any $c=1,\mathrm{\dots },N$, if $c\in A$ and $counts[c]=0$ ($c$ appears in all ${\overline{D}}_i$), or $c\notin A$ and $counts[c]=z$ ($c$ appears in all $D_i$), then $c$ is part of the result. Note that while one can simply rely on Algorithm 1 to intersect the color sets without the additional $\mathrm{\Theta }(N)$ cost, this algorithm should be executed twice: once for each operand in the union of Lemma 2. Given that $|A|+|\overline{A}|=|D_i|+|{\overline{D}}_i|=N$, this results in a time complexity of $\text{𝒪}(zN)$, which is worse than that of Algorithm 2, especially considering that $|D_i|\ll |A|\le N$.

It is also important to consider that accessing the sets sequentially to build the $counts$ array leads to better cache locality and further improves performance in practice.

In the general case, when not all color sets belong to the same group, we first apply Algorithm 2 to every group and then apply Algorithm 1 to the intermediate results.

This representation organizes color sets in two layers: the first layer is composed of meta colors, which are pointers to groups of partial color sets that form the second layer. Meta colors pairs of integers $(j,p)$, with $j$ being the group identifier, and $p$ the identifier of the partial color set inside the group.

Since groups form a partition, we can perform the intersection independently on each group and then join the results. Consequently, if group $j$ does not appear in a meta color set, it can be discarded directly, as this will lead to an empty group intersection. This means we can devise a two-layer intersection algorithm, where the first pass filters out groups that do not appear in all meta color sets, and the second one performs the intersection on the partial colors within each group. Moreover, instead of computing the intersection of all partial color sets, we can compute the intersection between the distinct partial color sets.

These ideas lead to Algorithm 3. If $ℱ$ represents the first-layer intersection (i.e., among the group identifiers, only), its time complexity is $\text{𝒪}\left({\sum }_{i=1}^z|M_i|+{\sum }_{(j,p)\in M_i:j\in ℱ}|P_{j,p}|\right)$. Although this complexity is higher than that of Intersect (Algorithm 1) in the worst case – due to the additional need to intersect meta color sets alongside the partial color sets that constitute the original color sets – it is always beneficial in practice on repetitive collections as we are going to see in our experimental analysis in Section 5.

The main drawback of Algorithm 4 is that, besides not being cache-friendly, it has to perform a membership query to every color set for each color encountered. In this way, all color sets are queried the same number of times, with potentially many negative results in the case of short color sets.

A better approach, which we call counting union, is to scan the color sets to build the score of each color incrementally. To do so, the scores are maintained in an array $scores[1..N]$, so that $scores[c]={\mu }_c$ at the end of the scan. Then, it is sufficient to scan this array and return all colors whose score is $\ge T$. This algorithm has a time complexity of $\mathrm{\Theta }\left(N+{\sum }_{i=1}^z|C_i|\right)$. It is worth noting that this complexity, despite being better then $\mathrm{\Theta }\left(z\left|{\bigcup }_{i=1}^z{C}_i\right|\right)$ in the worst case (all sets are disjoint), may still be worse when $z$ is small and the sets $C_i$ share many integers (i.e., their union is small). Nevertheless, our experiments in Section 5 show that this new algorithm is more efficient overall than Algorithm 4, especially thanks to the improved cache locality, as color sets are scanned linearly Yet, we can improve even more as follows.

As explained in Section 2.2, color sets with a density higher than $3/4$ are stored using their complement set. Therefore, reading the complement set is at least $3\times$ faster than reading the set directly. We can exploit this fact by noting that we are not interested in the actual score of a color, but in the difference between the score and the threshold. In fact, either adding ${\gamma }_i$ to the score of every color $c\in C_i$, or subtracting ${\gamma }_i$ from the threshold $T$ and from ${\mu }_c$ for all $c\in {\overline{C}}_i$, yields the same result. We call this technique dynamic thresholding, since we adjust the threshold in the presence of very dense sets, making $T$ a dynamic value instead of fixed. The following lemma formalizes this technique, and Algorithm 5 shows how to use it for threshold-union, together with the counting-union technique described above. (In the statement of the lemma and its proof in Appendix A, “v.d.” stands for “very dense”.)

This lemma implies that we can directly process the complement sets of the very dense color sets, which are sparse by definition, thus avoiding iterating over very large sets.

As similarly done for intersection, consider the case where the $z$ differential color sets belong to the same group and hence are relative to the same representative set $A$. Let the score of the group be $\gamma ={\sum }_{i=1}^z{\gamma }_i$. To perform the union, in this case, we leverage a variant of Lemma 2, as formalized in the following lemma (and proven in Appendix A).

In this way, it is sufficient to read the differential and representative sets once, plus the array of scores to determine the colors $c$ for which ${\mu }_c\ge T$, obtaining a total time complexity of $\mathrm{\Theta }\left(N+{\sum }_{i=1}^z|D_i|+{\sum }_{j=1}^r|A_j|\right)$, with $r$ being the number of groups in the index.

We can adapt the two-layer intersection algorithm described in Section 3 to also work for the threshold-union strategy, by considering the scores of the colors. In the first layer, we compute the score of each group and filter out all groups whose score does not reach the threshold, as it is impossible for their partial sets to gain a score $\ge T$. Then, in the second layer, for each group that “survives” the initial filtering, we deduplicate its meta colors before performing the union of the distinct partial color sets.

The efficiency of this process, unlike the ones described previously, critically depends on the value of $T$ (remember that this value is controlled via a parameter ): if $T$ is small, fewer groups are filtered away, increasing the computational load on the second layer. If ${ℱ}_T$ indicates the set of groups whose score is at least $T$, then the complexity of the algorithm is $\text{𝒪}\left(N+{\sum }_{i=1}^z|M_i|+{\sum }_{(j,p)\in M_i:j\in {ℱ}_T}|P_{j,p}|\right)$.

For all experiments, we use a $k$-mer length of $k=31$, which is common throughout prior work. For the threshold-union queries, we set the threshold to the customary value $\tau =0.8$. Reported times are relative to a second run of each experiment, where the indexes are loaded from the disk cache, benefiting mostly the larger indexes than the smaller ones. To avoid recording I/O overhead, the pseudoalignment output of all the experiments was written to /dev/null.

To better measure the impact of our algorithms, we tested them under a high-hit query workload, i.e., with a mapping rate (percentage of mapped reads over the total number of queried reads) of about $90\%$ or more, as it is more informative of the efficiency of the color set processing step. A low-hit workload, on the other hand, stresses the speed of negative $k$-mer lookups on the dictionary, and is therefore not relevant for the scope of this work. Figure 7 in Appendix B shows how query time is distributed among the different steps of pseudoalignment: as evident, processing the color sets is the most resource-intensive step for these high-hit workload queries (with the exception of GB, where processing the color sets is not the most time-consuming step because the sets are very short on average; see below). Although not reported herein, similar results were observed for the threshold-union.

For the rest of the paper, we use the following naming convention to distinguish between the different representations used: h-Fulgor employs the DA representation, while d-Fulgor and m-Fulgor employ repetition-aware compression schemes, the DR and MP representations, respectively. The md-Fulgor variant instead – our most succinct index to date – combines the two repetition-aware representations: it first partitions the color sets using the MP approach, then applies the DR scheme on the partial color sets of each group.

The evaluations were performed using the following datasets, as also used in prior work [13, 8]: 3,682 Escherichia Coli (EC) genomes from NCBI [1]; different collections of Salmonella Enterica (SE), ranging from 5,000 to 150,000 genomes from the “661K” collection [6]; 30,691 genomes assembled from human gut-bacteria samples (GB), from [17]. More details on these pangenomes can be found in [8, Table 1]. We remark that the GB dataset is much more diverse than the others: it features very short color sets (44 integers on average) but contains an order of magnitude more unitigs than the SE and EC datasets.

The queried reads consist of all FASTQ records in the first read file of the following nucleotide runs: SRR1928200 for EC, SRR801268 for SE, and ERR321482 for GB. These files contain around $6.5-7$ million reads each. The lengths of each read in the FASTQ file are: 100 bases (70 $k$-mers) for EC, 53 bases (23 $k$-mers) for SE, and 90 bases (60 $k$-mers) for GB.

Table 1 reports the results of the full-intersection algorithm, while Figure 4 provides a graphical comparison between the most significant improvements observed. Similarly, Table 2 and Figure 5 report the results of the threshold-union algorithm. In both tables, for each Fulgor variant, we show the query time before applying any optimization and after using the optimizations described in the previous sections, together with the RAM used while performing the queries. Importantly, the RAM usage does not change when using the new algorithms, so we report it once.

The clear result of these experiments is that all indexes are now substantially faster, by $2-3\times$ on average. More precisely:

• $\mathrm{■}$

  h-Fulgor, while being the fastest (but also largest) variant without optimizations, is now even faster thanks to the density-aware optimization, with speedup factors of $1.2\times$ for full-intersection, and $1.8\times$ for threshold-union.

• $\mathrm{■}$

  d-Fulgor is twice as fast for both pseudoalignment algorithms after the optimizations. As a result, its performance is now comparable to that of h-Fulgor – while it remains generally $20\%$ slower, it also consumes $3\times$ less space.

• $\mathrm{■}$

  The tailored m-Fulgor intersection algorithm leads to a consistent $20\%$ improvement over h-Fulgor. This result is remarkable considering it comes with a large reduction in memory usage, e.g., $8\times$ less space on the largest SE-150K collection. Table 3 from Appendix C motivates why this is the case: the number of meta colors is nearly three orders of magnitude smaller than the number of colors. Furthermore, by applying filtering and deduplication of the meta color sets, the number of colors intersected is reduced by $36-46\%$, significantly reducing the time spent in intersecting the partial color sets.

• $\mathrm{■}$

  md-Fulgor, which combines the benefits of both the DR and MP representations, is the representation that improves the most with the optimizations, reaching speedups of up to $3.2\times$ for both algorithms.