Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-05-22T05:36:25.077Z
Comments.Text: 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.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-05-22T05:37:19.247Z
Comments.Text: If you can program in Python, for this kind of task, I would recommend using pysam.
Comments.UserDisplayName: bli
Tags: <rna-seq><sam><samtools>
Title: How to extract reads that map exclusively to a single site with 1 or zero mismatches from BAM files
ThreadText: # 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
  
  