How to break broad peaks from H3K9me2 ChIP-Seq?

Question

[YLS1010]

Hi all,

I am having trouble for breaking broad peaks from H3K9me2 ChIP-Seq. I am considering to use PeakSplitter to split the broad peaks into subpeaks, but I don't know how to perform the differential expression analysis for individual peaks after that. Can anyone give me some advices? Thanks a lot.

Yan

Answer 1

[Sinji]

Not entirely sure what you mean by breaking broad peaks, but you could try calling peaks using a narrow peak caller like MACS2. Differential expression analysis can be done using DiffReps or any of R's many ChIP-seq differential expression packages.

[YLS1010]

Thanks a lot for the response. I may not clearly explain my question. Actually I got my result from our bioinformatic analyst, she use MACS2 for peak calling. From there I got narrow peaks and broad peaks results. She then use edgeR for differential expression analysis. The final result I got including the start and end position of each peaks in the chromosome, the fold change, p value and et al.. But some of those peaks from this file contain more than one gene. The p value and fold change data come from several genes. I can't make any conclusion base on this. This is why I will need to further break those peaks into individual genes and get differential expression result for each genes. Do you think this is possible? Thanks again.

[SP]

One peak belong to multiple gene, you have to be clear (may be paste small part of list here). while annotating a peak, most of the people consider the nearest TSS for gene annotation. you have to ask her, is there a reason for annotation with two genes? She might have given two genes for intragenic peaks because it is hard to annotate them with one gene.

[YLS1010]

Chr1_13371811-13417113 logFC: -0.16795205 logCPM:10.02651203 LR: 39.86880816 pValue: 2.72E-10 FDR: 1.64E-09

Chr1_13417920-13449994 logFC: -0.157807128 logCPM:9.787226265 LR: 33.24872197 pValue: 8.11E-09 FDR: 4.32E-08

Thanks again for the reply. Here are two examples for the peaks I mentioned. Both coordinates contain 8 genes respectively. The peaks cover entire genes...

[SP]

45kb and 32kb wide peaks, those are very big regions. I don't know what is your factor of chipping (transcription factor or histone). Of course, if you have such big region then you will have multiple genes associated with it. Easy fix could be to find TSS of a gene which is closest to the peak summit and annotate peak with this single gene.

[YLS1010]

Thanks a lot! It is H3K9me2 ChIP. Will try accordingly.

Answer 2

[SP]

I agree with sinji, you call peaks with MACS2 and without --broad parameter it will give you narrow peaks. Regarding calling for differential peaks, there are multitude of options. This paper sums it up quite well. Just to give you an idea, it is really important to have replicates of chip-seq for differential peak calling. If you have them one method could be DESeq2. For this you extract tag count for associated peak in all the chip-seq and run DESeq2 on it. Aware of: for normalization DESeq2 assumes that a lot of peaks will not change (at least 50%), if this condition is satisfied for your data then go ahead and try it. Hope it helps.

[YLS1010]

Thanks, will check with the reference.