I am trying to calculate the percentage of reads aligned to a reference genome.
I used samtools flagstat and this is what I have
24887959 + 0 in total (QC-passed reads + QC-failed reads)
18203689 + 0 secondary
210060 + 0 supplementary
0 + 0 duplicates
22937586 + 0 mapped (92.16% : N/A)
6474210 + 0 paired in sequencing
3237105 + 0 read1
3237105 + 0 read2
2537654 + 0 properly paired (39.20% : N/A)
4022194 + 0 with itself and mate mapped
501643 + 0 singletons (7.75% : N/A)
1475914 + 0 with mate mapped to a different chr
736257 + 0 with mate mapped to a different chr (mapQ>=5)
However, I noticed that the total number of reads is different to the input fastq files (3237105 read 1 + 3237105 read 2 = 6474210). After some reading, the numbers reported in samtools flagstat also include multiple mapped reads.
So I try this:
samtools view -c -F 260 999.sorted.PE.bam
options
-c count reads and print the total number
-F bitcode skip reads with bits present in 'bitcode'
-F 260 only output reads that are not unmapped (flag 4 is not set) and only primary alignment (flag 256 is not set)
4733897 is what I got so if I divide this number with the total read number from fastq files (6474210) I have 0.731 (x 100 = 73.1%).
Is this the correct number to report as a percentage of mapped read? Not the one that's reported in samtools flagstat (92.16%)?
•
written
by
rujiporn.sun •
I think you forgot to exclude supplementary alignments. If you do, you should find 69.9% mapping.
Number of input reads = total - secondary - supplementary
= 24887959 - 18203689 - 210060
= 6474210
Number of unmapped reads = total - mapped
= 24887959 - 22937586
= 1950373
Number of mapped reads (primary excluding supplementary and secondary) = input - unmapped
= 6474210 - 1950373
= 4523837
% mapping = mapped (only primary) / input*100
= 4523837 / 6474210 *100
= 69.9%
-F 260 in samtools view should exclude unmapped reads and non-primary reads. At least, I know I am on the right direction. Thank you for your comments and calculations.
Log in to answer this question.
What was the aligner used? It will probably have it's own, more accurate stats, especially about mapping (multimaped, one pair mapped etc etc).
Thank you for your suggestion. I used bwa mem. I will look into it.