Question: Ideas to split peaks into subpeaks - beyond PeakSplitter
gravatar for Sergio Martínez Cuesta
6 days ago by
Cambridge, UK

Here are some thoughts regarding how to break down bed files into smaller (perhaps more concrete) bed files using signal / read counts from e.g. wig or bedgraph files.

One possible scenario would be: you have identified peaks (treatment vs. control) using callers such as macs2, however after having a look at the signal in igv, you would like to further split your resulting peaks bed file into subpeaks e.g. it looks as if some identified peaks actually contain more than one peak.

    *               *
    **              **
   *****     *     ****
  *******   ***   ****** 

We have had some mixed experiences using PeakSplitter.jar to perform this so I am putting forward some possible alternatives:

(1) Split your original peak bed files in sliding windows using bedtools makewindows

    *               *
    **              **
   *****     *     ****
  *******   ***   ****** 

Then count reads in the sliding windows using bedtools coverage. Finally, perform differential analysis on the counts (treatment vs. control) obtained from the sliding windows using limma, edgeR or DESeq2.

(2) Go back to the peak calling step e.g. macs2 callpeak but now using the option --call-summits

macs2 callpeak -h

--call-summits   If set, MACS will use a more sophisticated signal
                 processing approach to find subpeak summits in each
                 enriched peak region. DEFAULT: False

Once you have the summits (1bp), extend them left and right e.g. 5-10bp using bedtools slop

(3) Use tools that work directly on the wig or bedgraph signal (filtered down to your peak bed regions), then identify local maxima and extend left and right as in (2), e.g. the routine bwtool find might be of use here.

Any other ideas?

chip-seq wig peak calling • 87 views
ADD COMMENTlink modified 5 days ago by Rashedul Islam220 • written 6 days ago by Sergio Martínez Cuesta40

Before discussing the technical side, what makes you think that (beyond eyeball inspection) these "subpeaks" are indeed meaningful and not just noise between significant events? What kind of data are this, ATAC/ChIP-seq?

ADD REPLYlink written 6 days ago by ATpoint3.8k

Good point, we do not know if they are meaningful / noise. These are ChIP-seq datasets. I think it comes down to the issue of finding narrower peaks to the ones found so far with existing arguments. Perhaps tweaking arguments in the peak callers would give me the finer detail I am looking for.

ADD REPLYlink written 3 days ago by Sergio Martínez Cuesta40

What's the marker of interest? TFs bind in narrow regions, obviously, whilst other markers 'spread' across a large region, even up to megabases in some cases

ADD REPLYlink modified 3 days ago • written 3 days ago by Kevin Blighe19k
gravatar for Alex Reynolds
6 days ago by
Alex Reynolds23k
Seattle, WA USA
Alex Reynolds23k wrote:

You might take a look at Eric Rynes' hotspot2. It applies a rigorous statistical approach to calling enriched regions.

ADD COMMENTlink written 6 days ago by Alex Reynolds23k
gravatar for Rashedul Islam
5 days ago by
Rashedul Islam220 wrote:

You can use FindER peak caller. This tool has the options for selecting minimum and maximum peak size. Along with it you can select the distance between two peaks to be merged together.

ADD COMMENTlink written 5 days ago by Rashedul Islam220
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 847 users visited in the last hour