<p>For much faster and stable<sup>$</sup> implementation, sambamba can handle this using following one-liner:</p>

<pre><code>sambamba view -t 12 -h -f bam -F "mapping_quality &gt;= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam
</code></pre>

<p>See <a href="http://lomereiter.github.io/sambamba/docs/sambamba-view.html" rel="nofollow noreferrer">manpage</a> and <a href="https://github.com/biod/sambamba/wiki/%5Bsambamba-view%5D-Filter-expression-syntax" rel="nofollow noreferrer">sambamba synatx for filters</a> for details on further usage.</p>

<p>$: Besides <code>grep</code> based search being slow, it may return false <del>positive</del> negative hits (@gringer's reply above). I also do not find <code>awk</code> based approach reliable because fields like XA, SA, etc. are optional fields in <a href="https://samtools.github.io/hts-specs/SAMv1.pdf" rel="nofollow noreferrer">SAM format</a>, and not <em>positionally</em> fixed fields.</p>

CC BY-SA 4.0
2018-06-22T06:05:36.037Z
In what way were the hits false positive?
gringer
CC BY-SA 4.0
2018-06-22T15:45:16.470Z
Hmm.. I guess it should be false negative because you were excluding those search patterns with `grep -v`, right?
Samir
CC BY-SA 4.0
2018-06-22T19:35:07.953Z
The only likely way grep could produce the wrong results is if someone modified the read names or a read group to contain one of those strings, which seems quite unlikely.
Devon Ryan
CC BY-SA 4.0
2018-06-22T20:38:06.180Z
agree, it is quite unlikely but nonetheless, grep based approach is inherently slow.
Samir
CC BY-SA 4.0
2018-06-22T20:58:35.220Z
The slow part is in the conversion from BAM to SAM by `samtools view`. The `grep` program is very quick at what it does.
gringer
CC BY-SA 4.0
2018-06-22T02:09:02.160Z
2018-06-22T15:48:08.723Z
2018-06-22T15:48:08.723Z
Samir
<p><strong>Update</strong> - as of <a href="https://github.com/samtools/samtools/pull/1346" rel="nofollow noreferrer">January 2021</a>, samtools can now do filtering based on an expression that includes tag variables. In this case, this expression can be used to exclude any reads that have either an XA or SA tag:</p>
<pre><code>samtools view -b mapped.bam -e '!([XA] | [SA])' &gt; unique_mapped.bam
</code></pre>
<p>For more details on the samtools expression parameter, see the samtools documentation:</p>
<p><a href="https://www.htslib.org/doc/samtools.html#FILTER_EXPRESSIONS" rel="nofollow noreferrer">https://www.htslib.org/doc/samtools.html#FILTER_EXPRESSIONS</a></p>
<hr />
<p>Original answer follows....</p>
<p>To exclude <em>all possible</em> multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use <code>grep</code> on the uncompressed SAM fields:</p>
<pre><code>samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b &gt; unique_mapped.bam
</code></pre>
<p>Explanation follows...</p>
<p>I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.</p>
<p>... but what if that were not possible?</p>
<p>My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [<code>samtools view -bq 1</code>] probably won't work.</p>
<p>Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a <em>lot</em> of repetitive sequence:</p>
<pre><code>ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    &lt;AABBEGABFJKKKIM7GHKKJK&gt;JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG&gt;D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;
</code></pre>
<p>BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.</p>
<p>However, the alternative locations are shown by the presence of the <code>XA</code> tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the <code>XA</code> tag will be able to exclude multiply-mapped reads. The <code>samtools view</code> man page suggests that <code>-x</code> will filter out a particular tag:</p>
<pre><code>$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    &lt;AABBEGABFJKKKIM7GHKKJK&gt;JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG&gt;D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67
</code></pre>
<p>... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the <em>read</em>. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.</p>
<p>In this particular case, I can use <code>grep -v</code> on the uncompressed SAM output to exclude alignment lines that have the <code>XA</code> tag (and re-compress to BAM afterwards, just to be tidy):</p>
<pre><code><span class="math-container">$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b &gt; output_filtered.bam
$</span> samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
&lt;no output&gt;
</code></pre>
<p>Hurray! reads filtered. As a little aside, this <code>grep</code> search has a fairly substantial computational load: it's looking for some string with the text <code>XA:Z:</code> somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads <code>HAXXA:Z:AWESOME!:&lt;readNumber&gt;</code>, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the <code>XA:Z:</code>.</p>
<p>Now I do a check for <em>any</em> duplicated read names, just to be sure:</p>
<pre><code><span class="math-container">$ samtools view output_filtered.bam | awk '{print $</span>1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986
</code></pre>
<p>Oh... damn. I wonder what they are:</p>
<pre><code>$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;
</code></pre>
<p>Aha! Supplemental alignments, which use the <em>official</em> <code>SA</code> tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, <em><strong>both of the alignments still have MAPQ scores of over 3</strong></em>. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because <em>both</em> read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.</p>
<p>Luckily, BWA appears to put the <code>SA</code> tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the <code>SA</code> search as an additional <code>grep</code> filter:</p>
<pre><code><span class="math-container">$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b &gt; output_filtered.bam
$</span> samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
&lt;no output&gt;
</code></pre>
<p>Done. Easy peasy! <code>&lt;/sarcasm&gt;</code></p>
<p>... that original &quot;high-effort&quot; solution of using a different aligner doesn't look so bad now.</p>

CC BY-SA 3.0
2017-06-07T06:15:14.663Z
that's a neat solution, and well argued. I don't think it's that high effort if you condense the whole thing to a few lines of shell script either. :)
roblanf
CC BY-SA 3.0
2017-06-07T06:43:25.843Z
The thing with using a MAPQ to filter out multimappers is that it uses the presumption of a multimapper being "equally well mapping", which for even the examples presented here won't be the case. I've set your reply as the accepted answer on biostars but asked op there to clarify a bit on what exactly is wanted.
Devon Ryan
CC BY-SA 3.0
2017-06-07T06:57:33.667Z
Yes, it's not hard if you know what to expect, which is why I've tried to run through all the problems that might happen.
gringer
CC BY-SA 3.0
2017-06-07T08:34:51.717Z
About the read names possibly containing "XA:Z:", a way to achieve robustness would be to use a proper SAM format parser, such as pysam, if one programs in python, rather than grep. I don't know how much slower this would be.
bli
CC BY-SA 3.0
2017-06-07T10:03:15.300Z
Worst case, you can write a couple lines of python or likely use Pierre's jvarkit to filter things out.
Devon Ryan
CC BY-SA 4.0
2017-06-07T04:21:28.520Z
2022-08-19T08:45:41.163Z
2022-08-19T08:45:41.163Z
gringer
<p><em>This is based on a question from <a href="https://www.biostars.org/u/39680/" rel="nofollow noreferrer">betsy.s.collins</a> on BioStars. The original post can be found <a href="https://www.biostars.org/p/256448/" rel="nofollow noreferrer">here</a>.</em></p>
<blockquote>
<p>Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the <a href="https://github.com/Generade-nl/TULIP" rel="nofollow noreferrer">TULIP</a> assembler/scaffolder, which works best for reads that map to unique genomic locations.</p>
</blockquote>
<p>The questioner is referring to genomically-unique alignments, which are something different from multiple possible alignments at a single location. If the only multiple alignment is at the same genomic location, those alignments are still unique in the genomic sense. The impossibility of finding the exact, correct, alignment is a well-known problem, and there are a few downstream methods (e.g. left normalisation) to make sure that multiple local alignments and/or sequencing errors have reduced impact.</p>
<p>The concept of &quot;uniquely mapped reads&quot; is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.</p>
<p>Uniquely mapped reads (i.e. a read that maps to a single location in the genome) are sometimes preferred when running analyses that depend on the quantification of reads, rather than just the coverage (e.g. RNASeq). Duplicated reads require additional filtering or normalisation on the analysis side, and most downstream programs won't consider concepts like &quot;a fifth of the read maps somewhere here on chromosome 5, and four-fifths of the read maps somewhere here on chromosome 14.&quot;</p>

<sam><phylogenetics><bwa>
Obtaining uniquely mapped reads from BWA mem alignment
# Obtaining uniquely mapped reads from BWA mem alignment

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

*This is based on a question from [betsy.s.collins](https://www.biostars.org/u/39680/) on BioStars. The original post can be found [here](https://www.biostars.org/p/256448/).*

> Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the [TULIP](https://github.com/Generade-nl/TULIP) assembler/scaffolder, which works best for reads that map to unique genomic locations.

The questioner is referring to genomically-unique alignments, which are something different from multiple possible alignments at a single location. If the only multiple alignment is at the same genomic location, those alignments are still unique in the genomic sense. The impossibility of finding the exact, correct, alignment is a well-known problem, and there are a few downstream methods (e.g. left normalisation) to make sure that multiple local alignments and/or sequencing errors have reduced impact.

The concept of "uniquely mapped reads" is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.

Uniquely mapped reads (i.e. a read that maps to a single location in the genome) are sometimes preferred when running analyses that depend on the quantification of reads, rather than just the coverage (e.g. RNASeq). Duplicated reads require additional filtering or normalisation on the analysis side, and most downstream programs won't consider concepts like "a fifth of the read maps somewhere here on chromosome 5, and four-fifths of the read maps somewhere here on chromosome 14."

**Tags:** sam, phylogenetics, bwa

---

## 2 answers

### Answer by username_0 on 2017-06-07. Score: 14 | Accepted answer

Update - as of January 2021, samtools can now do filtering based on an expression that includes tag variables. In this case, this expression can be used to exclude any reads that have either an XA or SA tag:
samtools view -b mapped.bam -e '!([XA] | [SA])' > unique_mapped.bam

For more details on the samtools expression parameter, see the samtools documentation:
https://www.htslib.org/doc/samtools.html#FILTER_EXPRESSIONS

Original answer follows....
To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields:
samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

Explanation follows...
I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.
... but what if that were not possible?
My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [samtools view -bq 1] probably won't work.
Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a lot of repetitive sequence:
ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.
However, the alternative locations are shown by the presence of the XA tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the XA tag will be able to exclude multiply-mapped reads. The samtools view man page suggests that -x will filter out a particular tag:
$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the read. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.
In this particular case, I can use grep -v on the uncompressed SAM output to exclude alignment lines that have the XA tag (and re-compress to BAM afterwards, just to be tidy):
$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'


Hurray! reads filtered. As a little aside, this grep search has a fairly substantial computational load: it's looking for some string with the text XA:Z: somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads HAXXA:Z:AWESOME!:, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the XA:Z:.
Now I do a check for any duplicated read names, just to be sure:
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

Oh... damn. I wonder what they are:
$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

Aha! Supplemental alignments, which use the official SA tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, both of the alignments still have MAPQ scores of over 3. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because both read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.
Luckily, BWA appears to put the SA tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the SA search as an additional grep filter:
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d


Done. Easy peasy! 
... that original "high-effort" solution of using a different aligner doesn't look so bad now.

> **Comments:**
> **username_3**: that's a neat solution, and well argued. I don't think it's that high effort if you condense the whole thing to a few lines of shell script either. :)  
> **username_2**: The thing with using a MAPQ to filter out multimappers is that it uses the presumption of a multimapper being "equally well mapping", which for even the examples presented here won't be the case. I've set your reply as the accepted answer on biostars but asked op there to clarify a bit on what exactly is wanted.  
> **username_0**: Yes, it's not hard if you know what to expect, which is why I've tried to run through all the problems that might happen.  
> **username_4**: About the read names possibly containing "XA:Z:", a way to achieve robustness would be to use a proper SAM format parser, such as pysam, if one programs in python, rather than grep. I don't know how much slower this would be.  
> **username_2**: Worst case, you can write a couple lines of python or likely use Pierre's jvarkit to filter things out.  

---

### Answer by username_1 on 2018-06-22. Score: 4

For much faster and stable$ implementation, sambamba can handle this using following one-liner:
sambamba view -t 12 -h -f bam -F "mapping_quality >= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam

See manpage and sambamba synatx for filters for details on further usage.
$: Besides grep based search being slow, it may return false positive negative hits (@gringer's reply above). I also do not find awk based approach reliable because fields like XA, SA, etc. are optional fields in SAM format, and not positionally fixed fields.

> **Comments:**
> **username_0**: In what way were the hits false positive?  
> **username_1**: Hmm.. I guess it should be false negative because you were excluding those search patterns with `grep -v`, right?  
> **username_2**: The only likely way grep could produce the wrong results is if someone modified the read names or a read group to contain one of those strings, which seems quite unlikely.  
> **username_1**: agree, it is quite unlikely but nonetheless, grep based approach is inherently slow.  
> **username_0**: The slow part is in the conversion from BAM to SAM by `samtools view`. The `grep` program is very quick at what it does.  

---

