Not matching Stats number

[satishkumar]

Hi ,

I have a 150 * 150 paired-end MiSeq data on clinical exome.

I performed mapping of the sample with the human genome (from UCSC). using BWA, and here is the stats below of the mapping :

2983558 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
2853021 + 0 mapped (95.62%:nan%)
2983558 + 0 paired in sequencing
1491768 + 0 read1
1491790 + 0 read2
2843530 + 0 properly paired (95.31%:nan%)
2847292 + 0 with itself and mate mapped
5729 + 0 singletons (0.19%:nan%)
3629 + 0 with mate mapped to a different chr
1709 + 0 with mate mapped to a different chr (mapQ>=5)

But the problem is
My total raw reads I fed in BWA was 2,982,430 (Where R1: 1491215, R2:1491215)
whereas after mapping the total reads count is 2983558, which is clearly above than what I fed in.

It will be great if I can get a feed back on why I am seeing different numbers as my total reads?

Thanks,
Satish

[Brian Bushnell]

It's likely that some of these are are secondary (flag 0x100) or supplementary (flag 0x800) alignments. If you remove those, you should end up with the expected number of reads.

[satishkumar]

Hi Brian Bushnell,

Thank you for your suggestion. Can you suggest a way to eliminate those flags? And also it would be great if you can help me to understand why am I seeing those secondary (flag 0x100) or supplementary (flag 0x800) alignments?

Thanks again,

Satish

[Brian Bushnell]

You can do that with samtools. The manual is here.

I think command would be:

samtools view -b -F 0x900 mapped.bam > clean.bam

Those flags indicate that the same read mapped multiple times to different places in the reference.

[satishkumar]

Thanks Brain ...!!

I did figure out the command which runs(yours was giving me the same number probably becuase it was not reading 0*900 properly):

samtools view -F 2048 -b mapped.bam > clean.bam

Thanks again,
Satish