BAMQL: a query language for extracting reads from BAM files
© The Author(s) 2016
Received: 8 April 2016
Accepted: 21 July 2016
Published: 11 August 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.
Descriptions of BAMQL connectives and predicates
Is satisfied if expr is not satisfied.
expr | expr
Is satisfied if at least one operand is satisfied.
expr ^ expr
Is satisfied if at exactly one operand is satisfied.
expr & expr
Is satisfied only if both operands are satisfied.
cond then then_expr else else_expr
If cond is satisfied, this expression will only be satisfied if then_expr is satisfied. If cond is not satisfied, then this expression will only be satisfied if else_expression is satisfied.
The read is paired in sequencing.
The read is mapped in a proper pair
The read is unmapped.
The mate is unmapped.
The read is mapped to the reverse strand.
The mate is mapped to the reverse strand.
The read matches the specified SAM flag.
The read is the first read in a pair.
The read is the second read in a pair.
The alignment is not primary.
The read failed platform/vendor quality control checks.
The read is either a PCR or optical duplicate
The alignment is supplementary.
The read is mapped to the specified chromosome.
Matches the read if the proability of error is less than the specified probability.
The mate is mapped to the specified chromosome.
Matches if both the read and it’s mate pair are mapped, but only if to different chromosomes.
Matches all reads that cover the specified position or any higher position.
Matches all reads that cover the specified position or any lower position.
Matches all sequences that cover the range of position from specified start to end, i.e. N matches any base.
Matches all reads that are classified as the IUPAC-style nucleotide n at the specified position, i.e. N matches N.
Matches all reads that exactly match the IUPAC-style nucleotide n at the specified position.
Matches the read group.
header ~ /regex/
Matches a Perl-compatible regular expression against the read’s header.
The read is chosen based on a uniform pseudo-random variable.
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.
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.
Comparison of the filtering capabilities of the various tools
Supports multiple filters
Supports logical connectives
Supports chaining filters
Filters reads by mapping
Filters reads by flags
Filters reads by sequence
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 %.
Comparison of performance for selecting mitochondrial sequences (mt)
Comparison of performance for selecting reads with base at a sequence position (nt)
Comparison of performance for selecting paired-end sequences (paired)
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].
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
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
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- 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.View ArticlePubMed CentralGoogle Scholar
- 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.Google Scholar
- 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.