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
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
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?