text: user73: This is based on a question from betsy.s.collins on BioStars. The original post can be found here.
  
  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 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."
  
  user2978: 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.
  
  user73: 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    <AABBEGABFJKKKIM7GHKKJK>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    <AABBEGABFJKKKIM7GHKKJK>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:]]'
  <no output>
  
  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!:<readNumber>, 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
  <no output>
  
  Done. Easy peasy! </sarcasm>
  ... that original "high-effort" solution of using a different aligner doesn't look so bad now.
  
original_text: 73: <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>
   2978: <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>
   73: <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>
  