How can I look at find 'open chromatin' from ATAC-seq data using MACS2

Question

[a.rex]

I am new to analysing ATAC-seq data.

As mentioned on ATAC-seq peak calling with MACS, there seem to be two ways to use MACS to analyse ATAC-seq data.

1. Utilising the --shift -100 --extsize 200 command This, I believe, is to find where the cutting sites are.

2. Utilising the --shift 37 --extsize 73 command This is to find the nucleosomes since the DNA is wrapped around nucleosomes is circa 147bp.

I specifically want to document regions of open chromatin, particularly enhancer regions.

I have read in the literature that open chromatin regions will give reads <100bp. If this is correct, do I have to filter my bam file to look at only this length.

Surely if I map all of the my data, not filtering for read size I will end up with genome-wide ATAC-coverage? Looking at my bam file following mapping to the genome, my read sizes range from 50bp-450bp (with a peak at around 70bp).

Answer 1

[spacemorrissey]

You can use MACS, but it is not the best tool for identifying open chromatin regions. I would recommend using Fseq: http://fureylab.web.unc.edu/software/fseq/ It was designed to analyze DNase -seq data, the predecessor of ATAC-seq. When analyzing data from either experiment, you do not really want to look for overlapping reads the way you do for Chip-seq. You want to look for regions with a high density of cut sights or insertion sights, and Fseq is designed to highlight those.

[a.rex]

Hello - how exactly do I use this for ATAC-seq. I have a bed interval file made from my mapped reads bam file. What else is the input? An example command would help.

Answer 2

[ATpoint]

All these macs2 options that people use, combining shift and extsizes in my experience do not make a notable difference. As spacemorrissey mentioned, unlike ChIP-seq (with typical short single-end reads) you want to scan ATAC-seq data for cutting sites within the accessable chromatin. I always use macs2 callpeak -g hs --nomodel -f BAMPE, which skips the macs-typical shifting model and piles up the entire paired-end fragment size. For reliable enhancer calling, ATAC-seq is the wrong assay anyway, because even though active enhancers are accessible, by far not every accessible region is an enhancer (promoters, silencers, insulators etc...). What exactly is the question you are working on?

[a.rex]

Hi ATpoint - I specifically want to identify enhancers (I also have H3K4me1 data for my model system). Problem is, I don't know whether to filter for short reads (i.e. <100bp). I see that some papers do this, others do not.

[ATpoint]

In that case, I would probably use the command I gave and default settings for H3K4me1, and then intersect the two files with bedtools. This then gives you candidate enhancer regions. H3K4me1 stretches are typically larger than ATAC-seq peaks, so I would not worry about subsetting the ATAC data, because the H3K4me1 will inflate the peak size anyway.