Software | Open | Published:
BAMQL: a query language for extracting reads from BAM files
BMC Bioinformaticsvolume 17, Article number: 305 (2016)
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.
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.
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.
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 . 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.
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.
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 . 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 . 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.
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) , Perl using Bio-SamTools (v1.41) , 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.
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 . 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.
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.
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].
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.
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
BAM, binary alignment/map; BAMQL, BAM query language
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools. Bioinformatics. 2009; 25(16):2078–9. doi:10.1093/bioinformatics/btp352.
Lattner C, Adve V. LLVM: a compilation framework for lifelong program analysis transformation. In: Proceedings of the International Symposium on Code Generation and Optimization, 2004. Washington: IEEE Computer Society: 2004. p. 75–86. doi:10.1109/CGO.2004.1281665.
Heger A, Belgard TG, Finkernagel F, Goodstadt L, Goodson M, Jacobs KB, Lunter G, Martin M, Schiller B. Pysam: Python Interface for the SAM/BAM Sequence Alignment and Mapping Format. via PyPI. https://github.com/pysam-developers/pysam. Accessed 2016-02-01.
Stein L. Bio-SamTools. via CPAN. http://search.cpan.org/~lds/Bio-SamTools. Accessed 2016-02-12.
Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015; 31(12):2032–4. doi:10.1093/bioinformatics/btv098.
Special thanks to Dr. Julia F. Hopkins for testing.
This study was conducted with the support of the Ontario Institute for Cancer Research to PCB through funding provided by the Government of Ontario. This work was supported by Prostate Cancer Canada and is proudly funded by the Movember Foundation – Grant #RS2014-01. PCB was supported by a Terry Fox Research Institute New Investigator Award and a CIHR New Investigator Award. This project was supported by Genome Canada through a Large-Scale Applied Project contract to PCB, Dr. Sohrab Shah and Dr. Ryan Morin.
APM and CML designed and implemented the software. APM did benchmarking and first draft of the paper. PS wrote the second draft, performed additional benchmarks, and made bug fixing updates to the software. MF, RGB and PCB supervised research. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
This directory contains tests written in BAMQL, C, Perl, Python, SAMtools and Sambamba languages. A binary file with the appropriate version is provided for Sambamba. (ZIP 1495 kb)
This file contains the complete results of the test-suite in a csv format. (XLS 32 kb)