<p>BWA-MEM can be used as a library. File <a href="https://github.com/lh3/bwa/blob/master/example.c" rel="noreferrer">bwa/example.c</a> shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header <a href="https://github.com/lh3/bwa/blob/master/bwamem.h" rel="noreferrer">bwa/bwamem.h</a> 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.</p>

CC BY-SA 3.0
2017-06-07T13:22:27.007Z
That is extremely cool.
Konrad Rudolph
CC BY-SA 3.0
2017-06-06T23:07:51.320Z
2017-06-06T23:07:51.320Z
user172818
<p><a href="http://seqan.de/" rel="nofollow noreferrer">SeqAn</a> supports BWT tables for use with their parametrisable alignment algorithms.</p>

<p>To use it, <a href="http://seqan.readthedocs.io/en/seqan-v1.4.2/Tutorial/SimpleReadMapping.html" rel="nofollow noreferrer">follow the general outline for building a SeqAn short-read aligner</a>, and use the <a href="http://docs.seqan.de/seqan/2.2.0/index.html?p=FMIndex" rel="nofollow noreferrer"><code>FMIndex</code></a> specialisation instead of — as in the example — the <code>IndexQGram</code>.</p>

CC BY-SA 3.0
2017-06-06T19:44:56.890Z
2017-06-06T19:50:47.123Z
2017-06-06T19:50:47.123Z
Konrad Rudolph
<p>If you're looking for a python API, nanopore have created one (although I am not sure how "maintained" it is) <a href="https://github.com/nanoporetech/bwapy" rel="nofollow noreferrer">https://github.com/nanoporetech/bwapy</a></p>

CC BY-SA 4.0
2019-04-29T00:39:30.740Z
2019-04-29T00:39:30.740Z
Michael Hall
<p>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 <a href="http://brinda.cz/publications/brinda_phd.pdf" rel="noreferrer">this thesis</a>). 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.</p>

<p>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 <a href="https://docs.python.org/3/library/subprocess.html#subprocess.Popen" rel="noreferrer">subprocess.Popen</a> for running the mapper and <a href="http://pysam.readthedocs.io/" rel="noreferrer">PySam</a> for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using <a href="https://github.com/mengyao/complete-striped-smith-waterman-library" rel="noreferrer">SSW</a> with its Python binding).</p>

<p>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).</p>

<p>There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:</p>

<ul>
<li><p><a href="http://seqan.rtfd.io" rel="noreferrer">SeqAn</a>. See the <a href="http://seqan.readthedocs.io/en/master/Tutorial/DataStructures/Indices/FMIndex.html?highlight=bwt" rel="noreferrer">FMIndex tutorial</a> and the <a href="http://seqan.readthedocs.io/en/master/Tutorial/Algorithms/Alignment/PairwiseSequenceAlignment.html" rel="noreferrer">Pairwise Sequence Alignment tutorial</a> for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a <a href="http://seqan.readthedocs.io/en/seqan-v1.4.2/Tutorial/SimpleReadMapping.html" rel="noreferrer">quick development of a read mapper</a>, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in <a href="https://github.com/seqan/seqan/tree/master/apps/yara" rel="noreferrer">YARA</a>.</p></li>
<li><p><a href="https://github.com/simongog/sdsl-lite" rel="noreferrer">SDSL-Lite</a>. This is a general library for succinct data structures. See the <a href="http://simongog.github.io/assets/data/sdsl-slides/tutorial#1" rel="noreferrer">tutorial slides</a> for an idea of how it works. For instance, <a href="https://github.com/iqbal-lab/gramtools" rel="noreferrer">GramTools</a> are based on SDSL.</p></li>
</ul>

CC BY-SA 3.0
2017-06-06T21:55:16.143Z
Never heard of YARA … better-known mappers also use SeqAn, e.g. both TopHat and Bowtie (though not for the mapping itself, as far as I know).
Konrad Rudolph
CC BY-SA 3.0
2017-06-12T20:12:35.160Z
Sigh, I was hoping to avoid writing BWA's SAM output to disk, but PySAM treats input strings as file names and explicitly doesn't support StringIO objects. Guess I'll have to look into named pipes or something like that.
Daniel Standage
CC BY-SA 3.0
2017-06-12T20:15:33.590Z
PySAM supports also file objects – see http://pysam.readthedocs.io/en/latest/api.html#pysam.AlignmentFile. AlignmentFile(filepath_or_object,....: "If filepath_or_object is a python File object, the already opened file will be used."
Karel Břinda
CC BY-SA 3.0
2017-06-12T20:17:45.197Z
Yeah, but running `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.
Daniel Standage
CC BY-SA 3.0
2017-06-06T20:17:34.890Z
2017-06-06T20:46:40.420Z
2017-06-06T20:46:40.420Z
Karel Břinda
<p>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?</p>

<p>Python bindings would also be great, but probably too much to expect.</p>

CC BY-SA 3.0
2017-06-07T08:51:12.840Z
Not C/C++, but just to mention that rust-bio seems to have BWT-related capabilities: <https://rust-bio.github.io/>
bli
<phylogenetics><bwa><read-mapping><api>
Library for computing BWT-based alignments
# Library for computing BWT-based alignments

(Asked by: username_0 on 2017-06-06)

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.

**Tags:** phylogenetics, bwa, read-mapping, api
> **Comments:**
> **username_1**: Not C/C++, but just to mention that rust-bio seems to have BWT-related capabilities: <https://rust-bio.github.io/>  

---

## 4 answers

### Answer by username_5 on 2017-06-06. Score: 12 | Accepted answer

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](http://brinda.cz/publications/brinda_phd.pdf)). 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](https://docs.python.org/3/library/subprocess.html#subprocess.Popen) for running the mapper and [PySam](http://pysam.readthedocs.io/) for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using [SSW](https://github.com/mengyao/complete-striped-smith-waterman-library) 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](http://seqan.rtfd.io). See the [FMIndex tutorial](http://seqan.readthedocs.io/en/master/Tutorial/DataStructures/Indices/FMIndex.html?highlight=bwt) and the [Pairwise Sequence Alignment tutorial](http://seqan.readthedocs.io/en/master/Tutorial/Algorithms/Alignment/PairwiseSequenceAlignment.html) 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](http://seqan.readthedocs.io/en/seqan-v1.4.2/Tutorial/SimpleReadMapping.html), but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in [YARA](https://github.com/seqan/seqan/tree/master/apps/yara).
* [SDSL-Lite](https://github.com/simongog/sdsl-lite). This is a general library for succinct data structures. See the [tutorial slides](http://simongog.github.io/assets/data/sdsl-slides/tutorial#1) for an idea of how it works. For instance, [GramTools](https://github.com/iqbal-lab/gramtools) are based on SDSL.

> **Comments:**
> **username_3**: Never heard of YARA … better-known mappers also use SeqAn, e.g. both TopHat and Bowtie (though not for the mapping itself, as far as I know).  
> **username_0**: Sigh, I was hoping to avoid writing BWA's SAM output to disk, but PySAM treats input strings as file names and explicitly doesn't support StringIO objects. Guess I'll have to look into named pipes or something like that.  
> **username_5**: PySAM supports also file objects – see http://pysam.readthedocs.io/en/latest/api.html#pysam.AlignmentFile. AlignmentFile(filepath_or_object,....: "If filepath_or_object is a python File object, the already opened file will be used."  
> **username_0**: Yeah, but running `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.  

---

### Answer by username_2 on 2017-06-06. Score: 9

BWA-MEM can be used as a library. File [bwa/example.c](https://github.com/lh3/bwa/blob/master/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](https://github.com/lh3/bwa/blob/master/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.

> **Comments:**
> **username_3**: That is extremely cool.  

---

### Answer by username_3 on 2017-06-06. Score: 4

[SeqAn](http://seqan.de/) supports BWT tables for use with their parametrisable alignment algorithms.

To use it, [follow the general outline for building a SeqAn short-read aligner](http://seqan.readthedocs.io/en/seqan-v1.4.2/Tutorial/SimpleReadMapping.html), and use the [`FMIndex`](http://docs.seqan.de/seqan/2.2.0/index.html?p=FMIndex) specialisation instead of — as in the example — the `IndexQGram`.


---

### Answer by username_4 on 2019-04-29. Score: 1

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>


---

