I am trying to have a bam file containing just the unique mapped reads from BWA.
I've used to check the mapping results:
samtools flagstat input.bam
and I got
151399675 + 0 in total (QC-passed reads + QC-failed reads) 1299223 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 150524060 + 0 mapped (99.42% : N/A) 150100452 + 0 paired in sequencing 75050226 + 0 read1 75050226 + 0 read2 141573428 + 0 properly paired (94.32% : N/A) 148964374 + 0 with itself and mate mapped 260463 + 0 singletons (0.17% : N/A) 6491734 + 0 with mate mapped to a different chr 5052530 + 0 with mate mapped to a different chr (mapQ>=5)
I know 150524060 means the unique and multiple mapped reads, but I want to make sure how many unique mapped I have.
So I used:
samtools view -F 0x40 -c input.bam
and I got 75664312
I am not sure if this number above is right.
Could someone please help me with this? Which command should I use to count unique mapped reads? and which to extract in a new bam file just the unique mapped (excluding the multi and unmapped ones)?
Thanks in advance.