We present bafiq , a command-line rust tool for repeated filtering and counting of sequencing reads stored in the BAM file format. This software can be easily incorporated into NGS computational workflows, especially in situations where one needs to repeatedly process mapped reads in large files based on their FLAG values. Commands in bafiq follow query language familiar from samtools and similar software where, however, indexing is only based on the location of reads in a reference genome. Our tool allows significant computational time savings in situations where a single BAM file is queried repeatedly based on diferent FLAG combinations, leveraging a lfag-based index. We made an efort to match or outperform the commonly used tasks as well as in sequence counting tasks.
wrapping these in Java [
CEUR Workshop
ISSN1613-0073 an SQL-like language for querying BAM files. The speed improvements are mostly based on parallel execution and code optimizations.
A well-known speed-centered addition to BAM/CRAM files is an index based on their
coordinatebased sorting. This functionality made its way into samtools and many other BAM processors as well.
The index divides the mapped and sorted reads in a BAM file into bins, which are then referenced in
the index for fast data retrieval from a given range. Genome browsers such as IGV [
In our work with BAM files we often needed to process them in various ways that, unlike the example above, depended on other components of the read data, not the coordinates. BAM format allows for reads to be marked using FLAGs and TAGs. For example, we may want to eliminate read pairs that do not map concordantly, or reads where their mate in the pair did not map. Or, we may want to only work with reads mapping with a certain minimum mapping quality (MAPQ value). While existing tools ofer arguments to select specific flags or tags, these operations do not take advantage of indexing, as if the queries and filtering operations were not supposed to be done more than once within a given workflow. We noticed that this is not always the case. One may want to count certain types of reads in the alignment files, repeat the counting procedure for diferent settings of a variable, then filter the file for only a subset of reads to create a ”cleaned” version for further work. Also in certain educational settings, an instructor may require students to query the same BAM file as part of a structured exercise.
To support higher speeds in these scenarios without compromising on the real-world data volume, we propose to create a new type of index file that can be stored along with the original alignment file, similarly to the coordinate-based (.bai) index that can be made using samtools. The bafiq tool described here is a command-line program written in rust that takes a BAM file as an input along with a number of filtering arguments. It then either counts or extracts the reads in SAM format that pass the query criteria. When queried, bafiq first looks for an existing flag index (.bfi) of the input file. If it does not exist, it is created upon first query. Any subsequent calls to bafiq on the same input BAM file will take advantage of the index and return results faster. Counts are returned immediately as they are part of the index, sequence retrieval based on the query follows standard I/O limitations with the added benefit of the index allowing selective compressed blocks skipping.
The SAM format [
Leveraging this observation makes pre-computing of all the flag combinations and the number of the corresponding reads upfront and embedding that into the .bfi index feasible. This way all count queries are resolved immediately as bafiq retrieves that count directly from the index. That also means the very first bafiq count query time is almost exclusively spent by building the index.
Learning the sequence number of given flag combinations in a BAM file is useful in itself however the downstream analysis often also requires the sequence extraction. For that reason the .bfi index stores the gzip file byte ofsets of the blocks which contain at least one sequence of the given flag combination. This allows for faster extraction of the actual sequences due to skipping of entire blocks during the retrieval process without processing the whole file start to finish. paired, mapped, read rev. strand, first in pair paired, mapped, mate rev. strand, second in pair paired, mapped, mate rev. strand, first in pair paired, mapped, read rev. strand, second in pair paired, mapped, mate rev. strand, second in pair, PCR/optical dup. paired, mapped, read rev. strand, first in pair, PCR/optical dup. paired, mapped, read rev. strand, second in pair, PCR/optical dup. paired, mapped, mate rev. strand, first in pair, PCR/optical dup. paired, read rev. strand, first in pair paired, mate rev. strand, second in pair
The BAM file is compressed using gzip which stores sequences in blocks. For reference a BAM ifle containing human chromosome 1 sequences obtained from a WGS (whole-genome sequencing) experiment of 30x read depth (common scenario) contains close to 70M sequences. A typical BAM block contains 180 sequences on average (see Section 6.1). That is roughly 380,000 blocks potentially (worst-case scenario) distributed among the flag combinations. The lists of blocks containing sequences of given flag combinations are referred to as bins. As each block can host sequences with diferent flag combinations it can be part of multiple bins. As a result the total number of block ofsets stored across all bins can grow quickly. Therefore we have briefly explored possible means of compression of the index itself to keep the size reasonable ideally without compromising bafiq’s speed.
As established earlier the list of bins is quite sparse (100 out of 4,096). Stripping empty bins is thus the ifrst obvious step. That in itself however does not provide large enough space-saving benefit so we explored further.
Each bin contains an ordered sequence of byte ofsets which are growing in nature. Thus storing the ofsets as deltas from the first saves a lot of space and is trivial to decode at runtime.
For the referenced example of chromosome 1 the .bfi index compression gains amount to ∼30% using just the 2 methods: Original size: Compressed size: Compression ratio: Space saved: 6,496,336 bytes (6.2 MB) 4,471,326 bytes (4.3 MB) 1.45x 31.2%
If we break it down by technique the ratio is as expected: Bin Sparsity: 129,280 bytes saved ( 6.0%) Delta encoding: 2,025,458 bytes saved (94.0%)
We have also explored the dictionary compression which would replace repeated ofset subsequences in the index with a literal to be expanded upon decompression. However the initial implementation did not scale well ( ( 2) in complexity) and was not the primary focus so we have decided not to include it. It might be a subject of future optimizations should the index size becomes an issue.
Since our main goal was to enable speedy exploratory analysis using flags we have explored and compared multiple strategies (Table 2).1 of reading the input BAM file and building the index.
The index build process boils down to (1) read, (2) decompress, (3) index and optionally (4) compress.
Each of the steps can be solved sequentially, in parallelized fashion or some combination of the two
with various benefits and tradeofs. Considering our goals we have decided to implement following 3
strategies: Channel Producer-Consumer (CPC) (utilizing crossbeam-channels[
The Work Stealing strategy leverages ryon[
The first computational challenge is the fast gzip block discovery in the BAM file with respect to the underlying OS I/O limitations. Two approaches were explored: (1) Discover-All-First and (2) Streaming. The Discover-All-First approach uses single thread to read the whole file (mapped into memory) organized into gzip blocks for downstream processing. While the Streaming approach is feeding the downstream gzip block processors continuously as they are being discovered.
By nature of the memory mapping (leveraging memmap2[
The next phase takes in discovered gzip blocks and decompresses them in memory for the downstream BAM records processing. Here the challenge is how to eficiently utilize existing decompression libraries (as the decompression itself is outside of the scope of this work).
We are leveraging libdeflater [
In this phase the just decompressed block becomes addressable for the contents. With the simple
operation of reading just the flag bit we opted-out of using the htslib-rust[
Both WS and CPC strategies build up local indices within scoped threads and merge at the end of decompression and extraction phases. The CM uses the same merge-tree algorithm although merges to the main index after every micro-batch is processed to minimize memory impact.
To define flag combination for querying the BAM sequences we opted for -f (required flags) and -F (forbidden flags) to keep familiarity.
As the SAM[
Hex
Read Flag
We have measured the CPU and memory utilization as well as overall time to finish a task of all 3 selected bafiq strategies alongside samtools as a reference. As bafiq’s advantage comes primarily after the index is built we also added our version of quick counting and filtering ( bafiq fast-count ) more resembling to what samtools view -c does to demonstrate the underlying performance.
Apart from the actual index building we bench-marked the sequence retrievals using the .bfi index.
Two short-read BAM files originating from a typical NGS WGS 30x pipeline were selected for benchmarking - chr22 (1.3GB) and chr1 (8.2GB) aligned to hg38 reference. As many of bafiq features build on parallelized processing all benchmarking runs had explicit –threads setting aligned with samtools for fair comparison.
All benchmark runs were performed either to produce or consume uncompressed indices as the final index size was not an issue (∼0.5% of the input file).
As seen in the index building benchmarks (Table 3) the best performing strategy was Work Stealing with its index building time close to samtools view -c across the diferent number of threads available, even slightly faster when both constrained to 2 threads. As the subsequent query for sequence count is fetched pre-computed from the index and resolve under 20ms we consider the index building a time-equivalent task to bafiq index + bafiq query .
We can clearly see that the cost of channel management overhead and work organization required by Channel Producer-Consumer strategy is very costly especially with limited threads but becomes less significant as the number of available threads grows.
Due to the memory-mapping feature both CPC and WS strategies can exhaust memory up to the original BAM file size which is not pracical for full-genome alignments. Because of that we have included more memory-conservative strategy (Constant Memory (CM)) which, although not as performant as WS, keeps a constant memory footprint (∼100 MB) independent of the input file size.
The sequence retrieval (Table 4) of bafiq view leverages stored gzip block byte ofsets (indexed for each of the flag combination) so that reading from the original BAM file does not depend on a pre-scanning step anymore and can limit its read to focused blocks. As each block stores 180 sequences on average traversing it is relatively cheap operation.
Where .bfi index performs quite well are use-cases of retrieving rare flag combinations relative to the bulk of records (chr22 0x4 ∼3,922). However its utility as a performance helper decreases with the abundance of matching reads (as seen in chr22 0x2 and 0x10 flags).
Despite that observed performance decrease in most measured scenarios the index-based retrieval was still faster than without it.
As the .bfi index file is stored in a binary form to save space we have developed a tool to quickly glance over the contents including basic index statistics. For example a .bfi index of chromosome 1 can look as follows: Total records: Total bins: Non-empty bins: Total blocks: Reads per block:
69095506 56 56 379169 min=35, max=310, avg=182.2 We have shown that an exploratory analysis of aligned reads based on FLAG combinations can vastly benefit from a dedicated index to obtain sequence counts (from s to ms) and for sequence retrieval scenario where the flag combination matches smaller subset. However in some cases where the flag combination matches large proportion of sequences of the BAM file the benefits of stored ofsets in current implementation diminish.
Nevertheless the code base and insights from individual presented strategies can serve as a demonstrator of utility of building narrowly focused indices. The present FLAG index scope could be also extended to include other fields from the header (such as TAG or MAPQ).
During the preparation of this work, the authors used Claude Sonnet 4 in order to: Assist in implementing parts of the rust code base. After using the tool, the authors reviewed and edited the content as needed and take full responsibility for the code content.
The authors have not employed any Generative AI tools for the publication content.
The source code for the bafiq tool is available via GitHub (https://github.com/honzakotrs/bafiq) .