2
$\begingroup$

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?

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

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.