The key function is call_broadpeaks
:
The description attached to the function says:
This function try to find enriched regions within which, scores are
continuously higher than a given cutoff for level 1, and link them
using the gap above level 2 cutoff with a maximum length of
lvl2_max_gap.
scoring_function_s: symbols of functions to calculate
score. 'p' for pscore, 'q' for qscore, 'f' for fold change, 's' for
subtraction. for example: ['p', 'q']
lvl1_cutoff_s: list of cutoffs
at highly enriched regions, corresponding to scoring functions.
lvl2_cutoff_s: list of cutoffs at less enriched regions,
corresponding to scoring functions.
min_length :minimum peak length,
default 200.
lvl1_max_gap : maximum gap to merge nearby enriched
peaks, default 50.
lvl2_max_gap : maximum length of linkage
regions, default 400.
Return both general PeakIO object for highly
enriched regions and gapped broad regions in BroadPeakIO.
To give some basic explanation, the algorithm (briefly) appears to be as follows:
Two separate levels of peaks are called, level 1 (a higher pval, ie more significant) and level 2 (a lower pval). Level 1 is controlled by -p
and level 2 is controlled by --broad-cutoff
. When each peakset is called, they are immediately linked by the max gap parameter for each set.
Then, assuming that all level 1 peaks should be inside level 2 peaks (this is an explicit assumption by MACS2), the algorithm groups level 1 peaks inside level 2 peaks to output a broad peak.
...
This has a few implications:
The broad peak calls really come from the level 2 peaks alone (+ linking). The level 1 peak calls allow you to distinguish sub peaks (so that you can have gapped peaks).
Aside from the linking, the broad peak calls would be the same as narrow peak calls, if you called both with the same pval threshold (for example, if you set --broad-cutoff 0.1
in broad peak mode, and the -p 0.1
for narrow peak mode)