Hi all,
I have 15 filtered BAM files from H3K27ac ChIP-Seq experiments, divided into disease and healthy groups. My goal is to fix the peak regions across all samples, compare these regions, and analyze how the read counts differ between the groups.
Initially, I used MACS2 for peak calling on each sample and then loaded the resulting .narrowPeak files into DiffBind for differential peak analysis. However, I ran into a problem: MACS2 calls peaks for each sample independently, resulting in different peak locations. For example, in sample 1, the first peak might start at chr1:1000 and end at chr1:1100, while in sample 2, the first peak might start at chr1:1050 and end at chr1:1180. This makes it challenging to establish fixed peak regions across all samples (e.g., a consistent region from chr1:1000-1100 for all samples).
I'm now wondering if there's a way to input the BAM files directly into DiffBind and have it both: Identify the peak regions across all samples, and Perform differential analysis between the groups at the same time.
In other words, can I skip the external peak calling step (e.g., using MACS2) and still achieve similar results using only DiffBind? If so, what would be the best approach to accomplish this?
Thank you in advance!
If you call peaks on each sample then you have to derive a consensus peak set. Be it by taking the intersecting ranges, all covered bp positions or anything like that.
You can also merge the bam files and call peaks on the merge. That in my experience will create a lot of nonsense calls, but at least you have s single peakset.
Or you use a group-aware caller such as Genrich. That will still require merging group peaksets.
Or you use csaw from Bioconductor which uses sliding windows for count comparisons (differential analysis).
Many options and no clear consensus, as usual.
This is a comment, and I don't want people to drag this to an answer, thanks.