Read on Twitter
ATAC peak calling with MACS2: I know this is a recurring problem but I got asked a few times recently, so I just put all info here for my own ref. We now routinely convert paired BAM to simple BED file and use "-f BED --shift -100 --extsize 200" for the peak calling. Why? 1/7
Like many other peak callers, MACS2 is for ChIP. In ChIP, the reads are flanking the actual binding of the protein. Therefore, many peak callers shift or extend reads towards the mid. of the fragments to reflect the actual binding. MCS2 uses extension, the "--ext" flag. 2/7
In ATAC/DNase, the mid of the fragments is not really what we are interested in. Instead, we are interested in the blue and red dots in the fig. which are the cutting sites of the enzyme. Those dots are the start (5' end) of your reads, so default of MACS2 doesn't fit here 3/7
Instead, we want to call peaks with fragments centred on the 5' end of your reads. We had this discuss about 6 years ago in the MACS2 google usergroup. Both @fooliu and @anshulkundaje provided excellent info. Check this link if haven't seen it before: https://groups.google.com/g/macs-announcement/c/4OCE59gkpKY/m/v9Tnh9jWriUJ 4/7
An update of MACS2 (ver 2.1.0 20140616) was made after the discussion, and you could freely manipulate the read positions with the combination of the "--shift" and the "--extsize" flags. Then, why covert paired BAM to simple BED to use "-f BED" ? 5/7
According to Issue #145 from the MACS2 GitHub https://github.com/macs3-project/MACS/issues/145 . When you set "-f BAM" or "-f BAMPE", MACS2 only takes the left read, ignoring the other. However in ATAC the 5' end of both R1 and R2 are of interest. Convert to BED will solve this problem, I guess?. 6/7
Not sure if all those modifications will make huge difference, but the results are different. I haven't tried HMMRATAC, Genrich, MACS3 etc. yet. Will do in future. Hope this helps to those who are new to this type of analysis. 7/7
