How to count mapping rates for ChIP-seq?

Question

[baijiangshan9726]

Hi, so how should we count for mapping rates for ChIP-seq? I used "samtools flagstat", but it always showed 100% mapping. For example, I have a fastq file, there are 136,192,060 lines in the file, so there should be 136,192,060/4=34,048,015 reads in this file. After mapping the reads to reference genome, Iused 'samtools flagstat' to see the mapping rate, it showed "11,892,792 + 0 mapped (100.00% : N/A)". But the mapped reads(11,892,792) is far lower than the original reads(34,048,015), how could the mapping rate be 100%? Thanks! And also, after trim-galore, the QC result showed the total sequences of the file is 9,477,948. So actually where does this number come from? Thanks so much!

[Ian]

Hello. What did you use to map the reads? I use Bowtie2, the log file of which gives useful information.

[baijiangshan9726]

yes, you are right. Thanks!

Answer 1

[Istvan Albert]

It all depends on how the file was filtered.

Various alignments may be filtered out - for example, unmapped reads are often removed or not even reported if you pass certain flags to the aligner.

When the unmapped reads are filtered out then the alignment rate becomes 100% :-)

[baijiangshan9726]

I see, thanks so much! So how do you usually count the mapping rate?

[Istvan Albert]

It is the number of primary alignments divided by the number of reads.

[baijiangshan9726]

Thanks! I mean what kind of tools you use to calculate this?

[Istvan Albert]

for example with samtools you would need to filter out unmapped (4), supplementary (2048) and secondary alignments (256) like so

samtools view -c -F 4 -F 2048 -F 256   alignments.bam

will produce a number of primary alignments