I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. As far as I can tell, though, both of these tools and their derivatives are invoked strictly via a command-line interface. Are there any libraries that implement BWT-based read alignment with a C/C++ API?
Python bindings would also be great, but probably too much to expect.
First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.
It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).
Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).
There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:
SeqAn. See the FMIndex tutorial and the Pairwise Sequence Alignment tutorial for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a quick development of a read mapper, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in YARA.
SDSL-Lite. This is a general library for succinct data structures. See the tutorial slides for an idea of how it works. For instance, GramTools are based on SDSL.
proc.communicate gives me BWA's SAM output as a string. It looks like none of the usual string trickery will work, and somehow I'll have to go through the file system to get PySAM to read it.
$\endgroup$
Commented
Jun 12, 2017 at 20:17
BWA-MEM can be used as a library. File bwa/example.c shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header bwa/bwamem.h contains basic documentation. Paired-end mapping is doable, but not well exposed. Several teams, including GATK and MS genomics, have been using bwa-mem as a library.
SeqAn supports BWT tables for use with their parametrisable alignment algorithms.
To use it, follow the general outline for building a SeqAn short-read aligner, and use the FMIndex specialisation instead of — as in the example — the IndexQGram.
If you're looking for a python API, nanopore have created one (although I am not sure how "maintained" it is) https://github.com/nanoporetech/bwapy
