BAMQL: a query language for extracting reads from BAM files

Abstract

Background

It is extremely common to need to select a subset of reads from a BAM file based on their specific properties. Typically, a user unpacks the BAM file to a text stream using SAMtools, parses and filters the lines using AWK, then repacks them using SAMtools. This process is tedious and error-prone. In particular, when working with many columns of data, mix-ups are common and the bit field containing the flags is unintuitive. There are several libraries for reading BAM files, such as Bio-SamTools for Perl and pysam for Python. Both allow access to the BAM’s read information and can filter reads, but require substantial boilerplate code; this is high overhead for mostly ad hoc filtering.

Results

We have created a query language that gathers reads using a collection of predicates and common logical connectives. Queries run faster than equivalents and can be compiled to native code for embedding in larger programs.

Conclusions

BAMQL provides a user-friendly, powerful and performant way to extract subsets of BAM files for ad hoc analyses or integration into applications. The query language provides a collection of predicates beyond those in SAMtools, and more flexible connectives.

Background

Binary Alignment/Map (BAM) provides a common format to hold large quantities of genomic read data after alignment to a reference genome. Reads are annotated with supplementary information not present in FASTA or FASTQ files, such as target position and chromosome, and information about the mate pair.

Often, only a subset of reads are needed for further analysis. SAMtools, the standard software for manipulation of BAM files, provides options to perform limited filtering of BAM files and collect matching reads [1]. The interface is neither expressive nor user-friendly: the user must indicate which bit flags they require not using names, but the numeric values of those flags. Further, the selection condition is restricted to a filter that describes which bit flags must match, which must not, a minimum quality score and a list of matching read groups. This limits query expressiveness; for example, selecting reads with poor quality is impossible.

General purpose programming languages can be used with BAM file manipulation libraries to filter a file. Since most subsetting operations are ad hoc, this is cumbersome and requires substantial boilerplate code. More concerning, there is a disparity between the checking done by the language and what is needed by the BAM APIs. For instance, it is possible to read only a subset of a file using an index, but this requires correctly matching chromosome names between the index request and the filtering expression. A mismatch may result in a silent omission of data that is difficult to detect due to the complex, mixed-language query.

To simplify the subsetting process, but retain the ability to have powerful queries, we developed BAM Query Language (BAMQL), a domain-specific language for matching BAM reads.

Implementation

BAMQL queries are specified by predicates and logical connectives. The predicates include ones which match BAM flags, mapping information, auxiliary read information, position and the sequence. For example, to output a BAM file with all reads which are on chromosome 1 and have a mate read, one could do the following: bamql -f input.bam -o output.bam ’chr(1) & paired?’. An outline of all available connectives and predicates can be found in Table 1.

Table 1 Descriptions of BAMQL connectives and predicates

BAMQL has a library which compiles a query into native code, then checks the reads in the input BAM file, saving the matched reads to a user-specified output BAM file. A recursive descent parser reads the query and the resulting structure is converted into native code using LLVM 3.4.0 or later [2]. A driver program reads the input BAM file, checking each read against the query and tracking the number of accepted and rejected reads, which it then writes to BAM files. All BAM file operations are done using HTSlib 1.3 [1]. The command line interface uses filename specifiers for input and output so, when utilizing pipes in the bash environment, /dev/stdin and /dev/stdout will allow for writing across pipes.

If an indexed BAM is provided, the query automatically determines if only a subset of the file needs to be analysed to improve performance. SAMtools has the same performance gain by specifying regions, but, there is the possibility of mismatching the regions specified and those in filter.

Results

To test the efficacy, we compared several equivalent queries written in BAMQL, SAMtools (v1.3) + GNU AWK (v4.0.1), Sambamba (v0.5.9), Python using pysam (v0.8.2) [3], Perl using Bio-SamTools (v1.41) [4], and C using HTSlib (v1.1). We compared execution times and the quantity of code written. Execution time is an imperfect metric, as it is dependent on various properties of the system at runtime, such as CPU load, disk operations, and the hardware specifications of the system. By keeping these factors equivalent when performing benchmarks, we can see relative differences between the tools. All tests are provided in the additional files [see Additional file 1] and examples of queries and a list of predicates can be found on the github wiki page [see Availability of data and materials]. Six tests were defined:

• all – copies the input to establish a baseline I/O cost.

• paired – keeps only reads which are paired.

• chrY – keeps reads on the Y chromosome, whether it is named “Y” or “24”, optionally prefixed with “chr”, irrespective of case.

• mt – keeps all reads which they or their mate pair are unmapped or mapped to the mitochondrial genome (named “M”, “MT” or “25”, optionally prefixed with “chr”, irrespective of case).

• nt – keeps all reads which have base C at position 13353 of the reference genome after alignment.

• chain – filters chromosomes 1, 2 and 3 into separate BAM files.

The tests were written as simply as possible. No profiling or intentional optimisations were done as this test code is meant to be a single-use ad hoc query.

Method-comparison

It must be noted that nt, mt and chrY tests could not be written with SAMtools alone and required output to be piped through AWK scripts for the actual filtering process. Similarly, Sambamba required an AWK script to implement filtering for the nt test.

The difference in user-friendliness between BAMQL and SAMtools is evident in the relatively simple chrY test. The BAMQL query for a chrY test (listing ??) is a one-line query since BAMQL recognizes the differing names for the human Y chromosome and checks them automatically.

The same query implemented using SAMtools (listing ??) however requires the query to be run using an AWK filter in order to manually check for the multiple ways Y chromosome could be labelled in the BAM file.

Sambamba was more expressive than SAMtools in this case. However, Sambamba requires chromosomes to be specified using reference ID rather than chromosome name [5]. Sambamba’s expressiveness also do not extend to complex queries such as the chain test that requires chaining filters. While BAMQL had a dedicated function for handling a chain of filters, each with an optional output file, Sambamba and SAMtools both required separate queries for each of the filters.

A summary of the differences in capabilities of the tools used for testing is shown in Table 2.

Table 2 Comparison of the filtering capabilities of the various tools

Performance

All tests were run three times on an Intel Core i5 3.4 GHz Ubuntu 14.04 workstation with a 10 GB indexed, sorted BAM file containing male human whole exome sequence as input. The reads from each implementation were compared, and were identical in all cases. The coefficient of variation of user time was < 4 % while that of system and wall times were noisier, up to 15 %.

The mt and nt tests, shown in Tables 3 and 4 respectively, are more processing-intensive while the paired and chain tests, shown in Tables 5 and 6 respectively, are less processing-intensive and keeps more sequences, causing it to be more output-intensive. Two notable outliers are Perl and Sambamba. Perl consumes much more CPU time than the other tools; this is evidenced by high user times in the mt test. Tests implemented with Sambamba alone, paired and mt, have a low wall time due to parallel processing capabilities of Sambamba. However, Sambamba+AWK and SAMtools+AWK have high system times, due to the overhead to cope with pipe operations. C, Python and Perl seem to incur a high overhead in the chain test due to multiple output files opened for depositing data. Overall, BAMQL, C and Python performed well. In the output-dominated tests, paired and chain, the differences in real processing become clouded by the cost of writing the output. Notice that the user and wall times show a smaller percentage difference in paired compared to mt. In the more logic-intensive test, mt, the languages separate out more.

Table 3 Comparison of performance for selecting mitochondrial sequences (mt)

Table 4 Comparison of performance for selecting reads with base at a sequence position (nt)

Table 5 Comparison of performance for selecting paired-end sequences (paired)

Table 6 Comparison of performance for applying chain of filters to separate by chromosomes (chain)

With the exception of the nt test, the BAMQL implementation is the most performant with 24 % to 90 % faster system time than the C implementation. In terms of code volume, BAMQL required much less code; typically half of the SAMtools+AWK and a tenth of the C implementation. An additional csv file shows the compete results for all tests [see Additional file 2].

Discussion

Since most BAM filtering operations are ad hoc, the runtime performance of BAMQL is not especially important, although it outperforms most standard approaches. We believe that BAMQL’s major impact is the simplicity of the queries and the automatic guards against common mistakes. In particular, BAMQL generates index requests from the BAM and does the filtering using the same query. This guards against accidental mismatches in chromosome names caused by independent index requests and filtering expressions (e.g. SAMtools with AWK filter). It also automatically corrects for common chromosomal name irregularities. These features make it easy to write correct queries and, by reducing boilerplate, makes them simple to read and understand.

Since the performance of BAMQL queries is good, it is also possible to use them in other contexts. Queries can be compiled to native machine code, making it trivial to integrate into a larger C or C++ program. There are queries that are outside the domain of BAMQL: BAMQL queries are stateless, while general purpose languages are not. Given the low overhead of integration, it is reasonable to split a stateful query between C or C++ and BAMQL.

Conclusions

BAMQL is designed to improve the efficacy and user-friendliness of the current tools like SAMtools and Sambamba. It provides a larger collection of predicates than SAMtools or Sambamba. These predicates can be joined together with a wider set of logical connectives, allowing for more expressive queries. By compiling queries into native code using LLVM, BAMQL is also able to perform as well as C but with significantly simpler code. BAMQL provides a powerful yet intuitive approach to ad hoc querying BAM formatted data.

Availability of data and materials

The complete test-suite as well as results supporting the conclusions of this article are included as additional files. The test-suite is also available on the homepage. Ubuntu packages for the latest version of BAMQL are available.Project Name: BAMQL Project home page: https://labs.oicr.on.ca/boutros-lab/software/BAMQL Version: 1.1 Operating System: Platform independent Programming language: C++, C Other requirements: LLVM 3.4.0 or higher, CLANG, HTSlib and libuuid are required to compileLicense: MIT Any restrictions to use by non-academics: None

Abbreviations

BAM, binary alignment/map; BAMQL, BAM query language