This is a test version of Biostars. For the public version, visit https://www.biostars.org.
ChIP-seq differential peak calling (Any way to set the window size for calling differential enriched peaks?)

Hi everyone. I'm analysing the H3K9me3 ChIP from the ChIP experiment I've conducted. However, the bioinformatic parts were did by a labmate who understand programming languages but lack of experience in analysing ChIP-seq.

We are now very struggling in calling the differential enriched peaks of my H3K9me3 ChIP between 2 conditions. The peaks were called by MACS2 either using default or --broad. Other parameters are as below. However, after calling the differential enriched peaks using bdgdiff, I found those "differential peaks" called are not as what I expected.

Question 1: Is there any way we can call the differential peaks based on the length of the broad peaks? In this way, it can tell if the histone mark is differentially enriched over a wide region (such as gene body).

Question 2: Sometimes, ChIP signal is actually enriched over a broad region and having a similar intensity within that region as shown in IGV. However, bdgdiff still called out some narrow regions within the broad one in one of the conditions. I'm afraid that this kind of enriched peaks are not biologically meaningful as the histone enrichment is actually similar across the region. Is there any way to eliminate this kind of "false positive" differential enriched peaks? Thanks!!

enter image description here enter image description here enter image description here enter image description here

chip-seq differential peak calling macs2 bdgdiff
• 6,431 views • link
written by chiefcat chiefcat

I am not an expert on this, but I attended a class taught by the author of MMDiff and it looked great for differential peak calling. Maybe worth a try? The paper is here, the software is a bioconductor package.

• 0 views • link
written by Fabio Marroni Fabio Marroni

I only played with the exercises they gave us in the class, which - if I well remember - were H3K4me3 data. It is possible that things stand as you say, this is beyond my knowledge, sorry.

• 0 views • link
written by Fabio Marroni Fabio Marroni

Do you have replicates per condition ?

• 0 views • link
written by GouthamAtla GouthamAtla

Hi geek, yes I've got biological replicates but they were pooled (probably during/ just before peak calling).

• 0 views • link
written by chiefcat chiefcat

If you have replicates, DiffBind is a really useful R package that make differential peak calling pretty straightforward. I was really impressed with it. If you don't have replicates, MAnorm is also useful for directly comparing two samples. Both will derive a consensus peak set to normalize and compare signal at these peaks across the samples/conditions/treatments. Sounds like they may be worth a look in your case.

• 0 views • link
written by jared.andrews07 jared.andrews07

1 answer

What I would do is:

  1. Pool all Chip bam files from one condition.

  2. Call broad/narrow peaks with relaxed cutoffs using MACS2. Like p-value of 0.01

  3. Repeat 1,2 for other condition.

  4. Make a union set of peaks from both conditions using bedtools mergeBed.

  5. Convert the bed file to SAF file and quantify the number of reads falling in the peaks using featureCounts for all samples.

  6. Feed the resulting matrix to DESeq2 and call the differential peaks.

• 0 views • link
written by GouthamAtla GouthamAtla

Thanks for the reply geek! Could you explain the purpose of point 4? And why you choose DESeq2 over bdgdiff from MACS2? Thanks!

• 0 views • link
written by chiefcat chiefcat

4th point makes sure you don't miss any regions from any condition. Ideally you want to consider all the regions enriched in any condition.

I don't know if bdgdiff considers the replicates.

• 0 views • link
written by GouthamAtla GouthamAtla

Log in to answer this question.

More posts like this