Samtools flagstat number of reads do not match actual total number of reads.
2
1
Entering edit mode
5.4 years ago

samtools flagstat results:

some number + 0 in total (QC-passed reads + QC-failed reads)

my this "some number" does not match actual number of reads in my paired fastq file. Is this expected?

samtools flagstat • 5.0k views
ADD COMMENT
1
Entering edit mode

Coming back to my post after quite sometime

This command should give the exact number of reads in your sample:

expr $(samtools view -f 4 your.bam -c) + $(samtools view -F 2308 your.bam -c)

Explanation

Count the number of unmapped reads

samtools view -f 4 your.bam -c ( -f : accept the reads with the flag 4 , this flag is for unmapped reads)

unmapped


Count of the number of reads which do not have flag 2308( -F : reject the reads with the flag 2308 , flag explanation in below image)

samtools view -F 2308 your.bam -c

sam-flags

IMAGE COURTESY : Explain SAM FLAGS

ADD REPLY
7
Entering edit mode
5.4 years ago

The number of reads reported is the number actually in the file. If your aligner produced secondary alignments then this will often be higher than the original number in the fastq files.

ADD COMMENT
0
Entering edit mode

Also, the aligner may throw out unaligned reads.

ADD REPLY
5
Entering edit mode
5.4 years ago

My 2p: It's good to keep in mind that SAM/BAM stores alignments not reads (in fact, sam stands for sequence alignment/map) even if unmapped reads can be present.

For this reason the simple question "how many reads have been aligned?" can be tricky to answer. A simple strategy is to count the reads that have not been aligned and get the difference with the raw read count from fastq. But if you want to know how many reads have been aligned with certain criteria (e.g. mapq > x, alignment score > y, properly paired etc) than you should consider also split reads.

I suspect the use of the word "read" in samtools flagstat causes a lot of misunderstandings in this respect.

ADD COMMENT

Login before adding your answer.

Traffic: 2237 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6