Scalable sequence database search using partitioned aggregated Bloom comb trees

Abstract

Motivation: The Sequence Read Archive public database has reached 45 petabytes of raw sequences and doubles its nucleotide content every 2 years. Although BLAST-like methods can routinely search for a sequence in a small collection of genomes, making searchable immense public resources accessible is beyond the reach of alignment-based strategies. In recent years, abundant literature tackled the task of finding a sequence in extensive sequence collections using k-mer-based strategies. At present, the most scalable methods are approximate membership query data structures that combine the ability to query small signatures or variants while being scalable to collections up to 10 000 eukaryotic samples. Results. Here, we present PAC, a novel approximate membership query data structure for querying collections of sequence datasets. PAC index construction works in a streaming fashion without any disk footprint besides the index itself. It shows a 3-6 fold improvement in construction time compared to other compressed methods for comparable index size. A PAC query can need single random access and be performed in constant time in favorable instances. Using limited computation resources, we built PAC for very large collections. They include 32 000 human RNA-seq samples in 5 days, the entire GenBank bacterial genome collection in a single day for an index size of 3.5 TB. The latter is, to our knowledge, the largest sequence collection ever indexed using an approximate membership query structure. We also showed that PAC's ability to query 500 000 transcript sequences in less than an hour.

Availability and implementation: PAC's open-source software is available at https://github.com/Malfoy/PAC.

Figure 1.

Left: A SBT structure. Right: A Bloom left-comb tree for the same database. Four datasets of a database are represented using grey shapes. Toy example Bloom filters are represented as bit vectors associated with shapes, and the content of each node is recalled. In the case of the Bloom comb tree, the inner nodes’ Bloom filters (highlighted zone in grey) are not explicitly represented but are encoded using an Aggregated Bloom filter instead (integer vector). Runs of 1’s corresponding to the integers are colored. For instance, the leftmost 1 is found at Levels 1 and 2, therefore a run of length two in the vector.

Figure 2.

(A) Two Bloom comb trees and their Aggregated Bloom filters. The leaves Bloom filters are in the middle, colored in grey. They are the Bloom filters representing the input datasets. The left-comb tree (top) and right-comb tree (bottom) are built by aggregating the Bloom filter’s k-mers. The Aggregated Bloom filter V_l represents the leftmost path of the top comb. Respectively, V_r represents the rightmost path of the bottom comb. We colored in green the second bit of leaves in order to show its encoding in the Aggregated Bloom filters. In the left comb, this bit is set to 1 only in the root for the longest branch; therefore, the value is 1 (green) in V_l. Conversely, the longest branch of the right tree has a run of five 1’s on the second position, hence, a 5 in V_r. (B) Using V_l and V_r from (A), we show how Aggregated Bloom filters define restrained search spaces in the leaves. Again, the leaves’ Bloom filters are grey, and the values of the vectors are printed vertically. Arrows represent the intervals given by each V, and the final search space at each position is the dotted grey area.

Figure 3.

(A) The content of one partition in PAC. The Aggregated Bloom comb trees constructed for a given partition are represented using Bloom filter leaves (in gray) and two Aggregated Bloom filters V_l and V_r. (B) A single query in the same partition as (A). The index has been inverted. A k-mer enters a PAC by finding the partition corresponding to its minimizer. V_l and V_r are loaded, and a search space (here $[3]$ only) is defined given the position obtained by hashing the k-mer (here the k-mer’s hash value is 2). Then, a slice is extracted from the inverted index at the given position, and bits are checked to find 1’s. Here, there is a single position, so necessarily an 1.

Figure 4.

Results on increasing human dataset sizes. We report total CPU hours for constructing the different indexes, including the Bloom filter construction and index constructions. X-axis is in $\text{log}\hspace{0.17em}2$ scale and Y-axis is in $\text{log}\hspace{0.17em}10$ scale. All methods were run with 12 threads.

Figure 5.

Results on increasing human dataset sizes. We report disk footprints, both temporary and final, for constructing the different indexes. Note that PAC does not use temporary disk. X-axis is in $\text{log}\hspace{0.17em}2$ scale and Y-axis is in $\text{log}\hspace{0.17em}10$ scale. All methods were run with 12 threads.

Figure 6.

Results on query batches. We present query CPU times computed on batches of 1 to 10k transcript sequences. The Y-axis is on a $\text{log}\hspace{0.17em}10$ scale. The queries were computed using 12 threads.