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?
