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

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?

  • $\begingroup$ 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. $\endgroup$
    – bli
    May 22, 2022 at 5:36
  • $\begingroup$ If you can program in Python, for this kind of task, I would recommend using pysam. $\endgroup$
    – bli
    May 22, 2022 at 5:37


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.