ATAC-seq peak calling

[igordot]

There has been a lot of discussion in the last few years regarding the correct approach to ATAC-seq peak calling with MACS (see #145). It still seems like the issue is not perfectly resolved. As suggested by @taoliu I am starting a discussion (new GitHub feature) which may be a better way to organize any updates or best practices.

[hoondy]

adding a reference to tweet: https://twitter.com/XiChenUoM/status/1336658454866325506

[hoondy]

FYI, here is a setting I use for ATAC-seq:

macs2 callpeak -t {BAM or BED files} -n {NAME} --outdir {OUTDIR} -f {BAM|BED} -g hs -q 0.01 --nomodel --shift -75 --extsize 150 --keep-dup all -B --SPMR

how are you guys processing your ATAC-seq?

[nservant]

I agree. So here, you are starting from a BED file to avoid the issue discussed before with Macs in BAMPE mode ?

[hoondy]

Yes, you can either convert to BED file or modify BAM and subtract 1 from samflag to treat alignment like single-ended. However, I read this post yesterday and I am rethinking about removing paired info: https://twitter.com/epigeneticsnerd/status/1337081679589101576

[nservant]

Thanks for the link. "(...) So if you are going to use MACS1/2 for ATAC-seq, and I will instead recommend using HMMRATAC or MACS3 (which will use the HMMRATAC algorithm for ATAC analysis), then the best settings are -BAMPE --call-summits. (...)"

So I guess I miss something, is there anything new in Macs3 now ?

[hoondy]

I don't think HMMRATAC is included in the current version (3.0.0a5) yet, #320

[dbrg77]

Thanks, all. This is a great place to put different thoughts together.

This is my current workflow:

# mapping
hisat2 -X 2000 -p 12 --no-temp-splicesite --no-spliced-alignment \
            -x hisat2_index/hg38 -1 {input.r1} -2 {input.r2} -3 1 \
            --summary-file {output.stats} | \
            samtools view -ShuF 4 -f 2 -q 30 - | \
            samtools sort - -T {wildcards.cell}_tmp -o {output.bam}

# remove dup using picard
java -jar -Xmx16g picard-2.20.3/picard.jar MarkDuplicates \
            INPUT={input.bam} OUTPUT={output.bam} \
            REMOVE_DUPLICATES=true \
            ASSUME_SORTED=true \
            METRICS_FILE={output.metric} \
            1> {log.out} 2> {log.err}

# Convert bam to bed
bamToBed -i {input.bam} > {output.bed}

# peak calling
macs2 callpeak -t {input.bed} -g hs -f BED -q 0.01 \
            --nomodel --shift -100 --extsize 200 \
            --keep-dup all -B --SPMR \
            --outdir macs2_pk -n ATAC \
            1> {log.out} 2> {log.err}

I think it is worth mentioning this ref. from the Liu lab: https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-9-537

It is not really ATAC-seq, but check Figure 1, they found 75 bp has better resolution. Maybe --shift -37 --extsize 74 would be better??

[dbrg77]

I don't know. Intuitively, the expected fragment length around nucleosomes should be 100 - 200 bp, right? I guess smaller extsize will give better res, but you need fragments of certain length to get smooth signals. 75 bp just turns out to be good.

[nservant]

The size of DNA around a nucleosome is of 146bp. So I guess the idea of 75 (or even 73 is some papers) is to centered in the middle of the nucleosome. So to me --shift -75 --extsize 150 would be perfect.

[nservant]

Coming back to this question of using --shit -37 --extsize 73 vs --shift -75 --extsize 150, as expected, using --shift 37 allows to identify smaller peaks, especially when we want to detect the precise Tn5 insertion sites. But overall, the effect remains very limited.
Here is an example ;

Of note, I also tried here the --call-summits option which is something used in some studies. It allows to find multiple summits for a given peaks, but honestly, I do not think this is really informative as the summit with the highest score is still the one we have without the option.

[alexyfyf]

This is pretty much what I used in my pipeline as well.
https://github.com/alexyfyf/atac_nf/blob/7f996b7de0e349c5a10dbbd75b2c266339517a3b/atac.nf#L341

[DineshRavindraRaju]

Thanks for the thread

[nservant]

Hi all. Thanks for the nice discussion.
One more question, does it make sense to shift the reads in the BAM file to correct for transposase insertion (+5/-4) before doing the peak calling with --shift and --extend ?
Or is it something which is actually also included in the macs peak calling ?

[hoondy]

I am not sure how others approach this problem, but I think it is sensible to account for Tn insertion. I do read trimming (5' only) by 9bp to account for Tn occupancy.

[nservant]

Thanks. So you are doing this trimming in any case before any analysis ... mapping, peak calling, etc. ?

[hoondy]

Trimming is done after the mapping.

Here's my workflow: 1) mapping using bwa-mem with mostly default parameter (this can be modified depending on your project) 2) alignment QC (remove mitochondrial, chimeric (mapped to different chromosomes or too far apart), duplicates, MAPQ<30, and trim 5' 9bp, etc) 3) throw processed alignments to macs2 for peak calling.

[harshini-c]

Hi all,

Thanks for this discussion. It really helps for those of us new to analysing ATACseq data.

I had a question related to calling Nucleosomes from ATAC-seq. The approaches described in this Discussion Thread -- they do not really filter the Nucleosme free regions (NFR) from the nucleosomes. But we say that the called peaks represent open chromatin. Some papers have filtered based on the fragment lenghts and calling peaks only on the NFR reads and converting the non-NFR bams to bigwig and visualising them. what do you think is a good approach to get nucleosome positions?

would greatly benefit from your thoughts.