This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Inconsistency between number of reads in Samtools flagstat and Macs2 total tags.

I am trying to call peaks on my ATAC data using MACS2 and realized that there is a huge inconsistency within the number of reads in my bam file and the total tags reported by MACS2.

I have seen previous posts that mention using "--keep-dup 'all'" option but I have deduplicated my bam using picard before macs. Below is the command am using and the outputs for macs and samtools flagstat.

macs2 callpeak -t PrimaryKeratinocytes.suf.bam -g hs -n test --nomodel --shift -51 --extsize 102 --q 0.05 --keep-dup 'all'

macs2

INFO  @ Fri, 17 Aug 2018 14:53:58: #1 tag size is determined as 51 bps 
INFO  @ Fri, 17 Aug 2018 14:53:58: #1 tag size = 51 
INFO  @ Fri, 17 Aug 2018 14:53:58: #1  total tags in treatment: 37171037

flagstat

52164544 + 0 in total (QC-passed reads + QC-failed reads)
14993507 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
52164544 + 0 mapped (100.00% : N/A)

Why is macs using only 37 million reads? How can I include all the reads for peak calling?

Thank you

chip-seq macs2 atac peak calling
• 2,709 views • link
updated by GenoMax GenoMax • written by ahmadammarmalik • ahmadammarmalik

The dedup bam file has 52164544 reads which are all mapped. I preprocess the bam files to keep only the unique reads that are mapped to the chromosome.

The flagstat shows this "52164544 + 0 mapped (100.00% : N/A)"

The way I preprocess the bam I know for sure it has 52 million reads unique mapped reads but the macs is only seeing 37171037.

• 0 views • link
written by ahmadammarmalik • ahmadammarmalik

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized.

This comment belongs under @swbarnes2's answer.

• 0 views • link
written by GenoMax GenoMax

1 answer

How many reads did you actually have after deduping?

Did you observe that 52164544-1499307 = 37171037

I'm guessing you really only have 37171037 reads, and 14993507 aligned to multiple places, so are in your file multiple times

• 1 views • link
written by swbarnes2 swbarnes2

The dedup bam file has 52164544 reads which are all mapped. I preprocess the bam files to keep only the unique reads that are mapped to the chromosome.

The flagstat shows this "52164544 + 0 mapped (100.00% : N/A)"

The way I preprocess the bam I know for sure it has 52 million reads unique mapped reads but the macs is only seeing 37171037.

• 1 views • link
written by ahmadammarmalik • ahmadammarmalik

samtools view myfile.bam | cut -f 1 | sort | uniq -c | wc -l Do you get 37M or 52M?

• 1 views • link
written by swbarnes2 swbarnes2

I get 41121155 reads. Still a few million more than what macs2 is reporting

• 1 views • link
written by ahmadammarmalik • ahmadammarmalik

I checked another bam and still, the number of reads are very inconsistent:

flagstat: 9613009

samtools view myfile.bam | cut -f 1 | sort | uniq -c | wc -l: 5987026

MAC2 total tags: 3252338

Is there a way to keep all the undup tags in my analysis?

• 1 views • link
written by ahmadammarmalik • ahmadammarmalik

Log in to answer this question.

More posts like this