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?
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?
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.
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