<p>I generated BAM files (sorted by coordinates) by aligning human RNA reads against the human reference genome using <code>BWA MEM</code> and <code>TopHat2</code>. I would like to know how I can retain only reads that mapped exclusively to a single site in the reference genome with one or zero mismatches for downstream analysis.</p>
<p>For retaining reads mapped to only one site I tried the following two commands, using <code>grep</code></p>
<pre><code>samtools view -h file1.sorted.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' |
samtools view -b &gt; file1.filtered_unqiue.sorted.bam
</code></pre>
<p>or filtering by <code>mapping quality</code></p>
<pre><code>samtools view -q 10 -b file1.sorted.bam &gt; file1.filtered_unique.sorted.bam
</code></pre>
<p>Then I filtered <code>file1.filtered_unique.sorted.bam</code> of either two commands above for mismatches with</p>
<pre><code>samtools view -e '[NM]&lt;=0' -O BAM -o
file1.filtered_unique_mismatches.sorted.bam   
file1.filtered_unique.sorted.bam
</code></pre>
<p>I am quite unsure about the correctness of the results.</p>
<p>I also came across a one-liner combining both efforts, where the BAM is first intersected with the coordinates of interest and then filtered for mismatches:</p>
<pre><code>samtools view -h file1.sorted.bam chr:start-end | egrep '(^@|(NM:i:(0|1)))' &gt; file1.filtered_unique_mismatches.sorted.bam
</code></pre>
<p>However, here I think this would also include multiple mapped reads since I only look for reads mapped in a specified region?</p>
<p>Does anyone have an advice what the best solution is?</p>

CC BY-SA 4.0
2022-05-22T05:36:25.077Z
Just a small tip: to help you see if your results are correct, you can visualize the resulting bam on a genome browser like IGV.
bli
CC BY-SA 4.0
2022-05-22T05:37:19.247Z
If you can program in Python, for this kind of task, I would recommend using pysam.
bli
<rna-seq><sam><samtools>
How to extract reads that map exclusively to a single site with 1 or zero mismatches from BAM files
# How to extract reads that map exclusively to a single site with 1 or zero mismatches from BAM files

(Asked by: username_0 on 2022-05-20)

I generated BAM files (sorted by coordinates) by aligning human RNA reads against the human reference genome using BWA MEM and TopHat2. I would like to know how I can retain only reads that mapped exclusively to a single site in the reference genome with one or zero mismatches for downstream analysis.
For retaining reads mapped to only one site I tried the following two commands, using grep
samtools view -h file1.sorted.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' |
samtools view -b > file1.filtered_unqiue.sorted.bam

or filtering by mapping quality
samtools view -q 10 -b file1.sorted.bam > file1.filtered_unique.sorted.bam

Then I filtered file1.filtered_unique.sorted.bam of either two commands above for mismatches with
samtools view -e '[NM]<=0' -O BAM -o
file1.filtered_unique_mismatches.sorted.bam   
file1.filtered_unique.sorted.bam

I am quite unsure about the correctness of the results.
I also came across a one-liner combining both efforts, where the BAM is first intersected with the coordinates of interest and then filtered for mismatches:
samtools view -h file1.sorted.bam chr:start-end | egrep '(^@|(NM:i:(0|1)))' > file1.filtered_unique_mismatches.sorted.bam

However, here I think this would also include multiple mapped reads since I only look for reads mapped in a specified region?
Does anyone have an advice what the best solution is?

**Tags:** rna-seq, sam, samtools
> **Comments:**
> **username_1**: Just a small tip: to help you see if your results are correct, you can visualize the resulting bam on a genome browser like IGV.  
> **username_1**: If you can program in Python, for this kind of task, I would recommend using pysam.  

---

## 0 answers

