How to compare peaks (called by MACS2) within a specific genomic area for two samples?

Question

[Amlan]

Hello,

I have two ATAC-Seq peak samples for the same subject. I want to quantify the ATAC-Seq peaks within a defined area around TSS for each of those samples. Then I want to compare the scores to find out which of the two samples has more open chromatin region within that area around TSS?

So, is there any way to find a score for all peaks within a specific genomic area?

Thanks in advance

Answer 1

[Alex Reynolds]

Given a set of intervals and a set of peaks with some score value, you could use BEDOPS bedmap --max-element or bedmap --max to get the maximum scoring peak over the interval, or just the maximum score, respectively:

$ bedmap --echo --max-element --delim '\t' intervals.bed peaks.bed > answer.bed

If you're looking for measured size or other peak features, there are other operators like --echo-map-size, --count, etc. that offer mapping of those results. Feel free to clarify what you're after, or take a look at the documentation or bedmap --help.

Answer 2

[BioinfGuru]

Im working on ATAC-seq data at the moment.

I've used ChipPeakAnno to annotate the peaks in the narrowPeaks file (MACS2 output) for genomic features (TSS etc). You could then do a simple count of the number of peaks annotated for each feature.

Have you considered using DESEQ2 for getting differential peaks first? That is a far more sophisticated package for comparing samples.