Crate sbwt 

Introduction

This crate contains an implementation of the Bit Matrix SBWT data structure, as described in Small Searchable k-Spectra via Subset Rank Queries on the Spectral Burrows-Wheeler Transform, for the DNA alphabet ACGT. A CLI for the main features of the library can be found at sbwt-rs-cli. The data structure uses a variant of the Burrows-Wheeler transform to compress a set of k-mers in way that allows fast lookup queries. If the input k-mers are consecutive k-mers from longer underlying sequences, the index takes typically around 5 bits per distinct k-mer, supporting k-mer lookup queries at a speed of around 1 μs / k-mer on modern hardware.

Queries can be further sped up by using the Longest common suffix array, (see here) taking roughly log(k) bits of space per k-mer. The LCS array also enables the computation of k-bounded matching statistics and shortest frequency-bounded suffixes (both described here). Two SBWTs can be combined with set operations: the union, the intersection and the difference (described here). Finally, the crate provides an interface for traversing the node-centric de Bruijn graph of the k-mers, and an interface for traversing the variable-order de Bruijn graph.

API Quick start

use sbwt::*;
use std::io::BufReader;
use std::io::BufWriter;
use std::fs::File;
use std::path::Path;
 
// Build the sbwt
let seqs: Vec<&[u8]> = vec![b"AACTGACTGATCGTCTTGACTCGTTTATCTACGGT", b"ACTGACAGCTCTGCGATGCGA"];
let seq_stream = sbwt::SliceSeqStream::new(seqs.as_slice());
let k = 6;
let (sbwt, lcs) = BitPackedKmerSortingDisk::new_from_slices(&seqs, k)
    .mem_gb(2)
    .dedup_batches(false)
    .temp_dir(Path::new("./temp"))
    .n_threads(4).build_lcs(true).add_rev_comp(true)
    .run();

// Query a k-mer
let query_kmer = b"GACTCG";
if sbwt.search(query_kmer).is_some() {
    println!("{} is found", &String::from_utf8_lossy(query_kmer));
}

// Query all k-mers of a longer query using k-bounded matching statistics
let lcs = lcs.unwrap(); // Ok because we used build_lcs(true)
let streaming_index = StreamingIndex::new(&sbwt, &lcs);
let long_query = b"TGATACGTCTTAGTGACTCGTTT";
for (i, (len, range)) in streaming_index.matching_statistics_iter(long_query).enumerate() {
    // Kmer ending at long_query[i] exists iff len == k
    println!("Longest match ending at {} has length {} and colex range {:?}", i, len, range);
}

// Write index to disk for later use
sbwt.serialize(&mut BufWriter::new(File::create("index.sbwt").unwrap())).unwrap();
lcs.serialize(&mut BufWriter::new(File::create("index.lcs").unwrap())).unwrap();
 

Reverse complements and input preprocessing

If your input is very large, you may want to preprocess it to remove duplicate k-mers reduce construction time, memory and disk space. We recommend using a specialized tool such as GGCAT, Cuttlefish or BCALM2 for this purpose.

This crate considers a k-mer distinct from its reverse complement. In the use case where a k-mer is considered equal to its reverse complement, you need to either feed both directions to the index, or feed just one direction but query each k-mer both ways. The former approach can be easily implemented by enabling reverse complements in the construction algorithm. For the latter approach, we recommend using one of the unitig construction tools mentioned above to turn the data into canonical unitigs containing each k-mer in only one orientation. After this, it is very important to reorient the unitigs to maintain a consistent orientation for neighboring k-mers since as explained here, each k-mer without a direct predecessor increases the index size. The payoff of this approach is up to two times smaller index size compared to indexing both orientations explicitly, with the drawback that now queries need to be run in both orientations.

De Bruijn graph operations

With the addition of two extra bitvectors, the data structure supports traversal on the node-centric de-Bruijn graph of the input k-mers. The set of nodes is the set of distinct k-mers, and there is an edge from x to y iff x[1..k) = y[0..k-1). The label of the edge is the last character of y. The graph is not aware of reverse complements, so if you want to traverse the bi-directed de Bruijn graph that includes edges where one or both of the endpoints are reverse complemented, you need to take care of that logic yourself. See Dbg for the API.

Variable-order de Bruijn graphs

An SBWT of order k also contains the de Bruijn graphs of all orders k’ ≤ k. The VoDbg struct exposes them as a single variable-order de Bruijn graph, where the order can be changed on the fly during a traversal in addition to the usual node-centric operations. This requires the LCS array, select support, and a previous/next smaller value structure. See the vodbg module for the API.

Set operations

Two indexes of the same order k can be combined into an index of the union, the intersection or the difference of their k-mer sets, All three start from a MergeInterleaving, which can also be queried for the sizes of the three sets without materializing the result index. These operations are described in Fast Set Operations for Compact k-mer Sets.

Construction algorithms

The crate provides four construction algorithms:

• BitPackedKmerSortingDisk: sorting of bit-packed k-mers using lots of temporary disk space. This is the lowest RAM option, but can be slow.
• BitPackedKmerSortingMem: Sort all k-mers in memory. This is the fastest algorithm, but uses a lot of RAM.
• BuildByBoundedSuffixSort: Runs a k-bounded suffix sort with SIMD suffix comparison. Good for short and medium values of k (k < 1000).
• BuildByLibsais: builds from a suffix array of a concatenation of the input strings, computed with the libsais C-library. Only available when the optional libsais feature is enabled. The running time is O(n), independent of k. Good for very large values of k (k > 1000).

The two bit-packed k-mer sorting algorithms stream the input and support k up to 256. The two suffix-sorting algorithms instead hold the concatenation of the input sequences in memory, and have no upper bound on k.

Details on the space usage of the index

The index exploits overlaps between k-mers to encode them in small space. We say that a k-mer x is a source k-mer if it has no incoming edges in the node-centric de Bruijn graph of the input k-mers S, that is, there does not exist a k-mer y ∈ S such that y[1..k) = x[0..k-1). The number of bits in SbwtIndex<SubsetMatrix> is 5(n + n’) plus a small constant, where n is the number of distinct k-mers in the dataset, and n’ is the number of nodes in the trie of all prefixes of length k-1 of all source k-mers (See here for more details on the inner workings of the SBWT to understand what is going on). When the k-mers are from biological sequences, and unitigs are oriented consistently, the term n’ is typically negligible, but if the k-mers are for example a randomly sampled subset, then the benefit of overlaps is lost, and the term n’ dominates. In the worst case, n’ can be up to n(k-1) + 1.

Limitations

The implementation only supports the DNA alphabet ACGT. For best compression, the input k-mers should originate from a longer underlying sequence so that sbwt is able to exploit the overlaps for better compression. For non-overlapping k-mer sets, a simple hash table is likely a better choice.

Citation

If you use the SBWT in your research, please cite as follows:

@inproceedings{alanko2023small,
  title={Small searchable k-spectra via subset rank queries on the spectral Burrows-Wheeler transform},
  author={Alanko, Jarno N and Puglisi, Simon J and Vuohtoniemi, Jaakko},
  booktitle={SIAM Conference on Applied and Computational Discrete Algorithms (ACDA23)},
  pages={225--236},
  year={2023},
  organization={SIAM}
}

Modules

benchmark

Benchmarking queries on an existing SbwtIndex.

build_by_suffix_sorting

Construction algorithms for an SbwtIndex via suffix sorting. Good for large k.

dbg

de Bruijn graph operations using SbwtIndex.

sbwt_index_variant

A wrapper for SbwtIndex<SS> that can hold an index with any SubsetSeq implementation that is defined in the crate.

vodbg

Variable-order de Bruijn graph operations using SbwtIndex.

Structs

BitPackedKmerSortingDisk

A construction algorithm based on sorting of bit-packed k-mers using temporary disk space.

BitPackedKmerSortingMem

A construction algorithm based on sorting of bit-packed k-mers in entirely in RAM. Faster and scales better with parallelism than BitPackedKmerSortingDisk, but takes more RAM.

BuildByBoundedSuffixSort

A construction algorithm based on sorting the k-bounded contexts of the input. Unlike BitPackedKmerSortingMem and BitPackedKmerSortingDisk, this algorithm holds a concatenation of the input sequences in memory, and is not limited to k <= 256. Runs in time O(nk), where n is the length of the input.

BuildByLibsais

A construction algorithm that uses the libsais crate to build a suffix array of the concatenation of the (reversed) input sequences, which crate::build_by_suffix_sorting::build then turns into the SBWT. Like BuildByBoundedSuffixSort, this algorithm holds a concatenation of the input sequences in memory, and is not limited to k <= 256.

FastXReader

A crate::SeqStream over the sequences of a FASTA or FASTQ file.

LcsArray

An array that stores the lengths of the longest common suffixes of consecutive k-mers in colexicographic order.

MatchingStatisticsIterator

An iterator that produces values of the matching statistics one by one. Create with StreamingIndex::matching_statistics_iter(). Useful for streaming over the matching statistics without having to keep the values in memory all at once. Using Iterator::collect on the iterator will give the same vector as StreamingIndex::matching_statistics.

MergeInterleaving

An interleaving plan for combining two SbwtIndex structures with a set operation.

PrefixLookupTable

A table storing the SBWT intervals of all 4^p possible p-mers.

SbwtIndex

The SBWT index data structure. Construct with BitPackedKmerSortingMem or BitPackedKmerSortingDisk. For the SubsetSeq trait implementation, we recommend using the bit matrix implementation SubsetMatrix.

SeqStreamWithPossiblyRevComp

Wraps a crate::SeqStream so that, if enabled, every sequence is followed by its reverse complement.

SliceSeqStream

Creates a crate::SeqStream out of a slice of ASCII sequences.

StreamingIndex

An index that uses right extensions and left contractions to find matches in a streaming fashion. The Copy trait is implemented because this struct contains only references to components of the index, so it can be copied around cheaply.

SubsetMatrix

An implementation of SubsetSeq with a matrix of sigma indicator bit vectors: the i-th bit of the j-th bit vector is 1 if and only if the i-th subset contains the j-th character. Rank and select queries are reduced to bit vector rank and select queries on the indicator bit vectors.

VecSeqStream

Creates a crate::SeqStream out of a slice of ascii vectors.

Enums

Orientation

Constants

CARGO_API_VERSION

Traits

ContractLeft

Contracting a search pattern from the left.

ExtendRight

Extending a search pattern to the right.

SeqStream

A stream of ASCII-encoded DNA-sequences. This is not necessarily a standard Rust iterator because we want to support streaming sequences from disk, which is not possible with a regular iterator due to lifetime constraints of the Iterator trait.

SubsetSeq

This trait represents a sequence of subsets from alphabet {0, 1, …, sigma-1}, where sigma is the alphabet size. The trait provides access to the subsets and rank and select queries for the elements inside the subsets.

Functions

build_from_kmers_on_disk

Build the index from the files written by sort_and_dedup_kmers_into_file. Pass first_mers_file if and only if that function was called with add_all_dummy_paths set. The input files are not deleted.

difference

Computes the set difference index1 \ index2: the SbwtIndex containing exactly the k-mers that are present in index1 but not in index2.

intersect

Computes the intersection of two SbwtIndex structures. The intersection k-mers are those present in both index1 and index2. This mirrors merge in structure but uses AND logic for incoming edges instead of OR, and restricts output positions to those shared by both SBWTs.

is_dna

Returns whether the given ASCII character is one of A,C,G,T,a,c,g,t.

load_from_cpp_plain_matrix_format

Loads an index that was previously serialized either with the C++ API https://github.com/algbio/SBWT, or the associated CLI. Supports only version v0.1.

merge

Merge index1 and index2 according to interleaving. After the merge, a PrefixLookupTable with prefix length new_prefix_lookup_table_length will be added to new index. The number of threads used in the merge is n_threads. Indexes are passed in as Arcs because if those are the only existing references, this function can free the input SBWTs early which lowers the memory peak. Passing as Arc also allows for use cases where the caller still wants to hold onto the sbwts: in that case dropping the Arcs here will not free the memory. Same goes for the interleaving.

optimize_unitig_orientation

Given an iterator of canonical unitigs of the node-centric de Bruijn graph of order k, returns a vector of orientations, one for each sequence, aiming to minimize the number of unitigs which do not have an incoming edge in the de Bruijn graph of order k. Canonical unitigs means that each k-mer occur in only one orientation. The input can be a VecSeqStream, a SliceSeqStream or any struct implementing SeqStream.

reverse_complement_in_place

Reverses the given ASCII DNA sequence and replaces each nucleotide with its complement.

sort_and_dedup_kmers_into_file

Sort and deduplicate the reverse k-mers of the input sequences into a file on disk, to be later consumed by build_from_kmers_on_disk. Returns the path to the k-mers file, and if add_all_dummy_paths is set, also the path to the file of first k-mers of each input sequence. The returned files are NOT deleted: it is the caller’s responsibility to delete them once they are no longer needed.