Answers.Body: <p>The key function is <a href="https://github.com/taoliu/MACS/blob/master/MACS2/IO/CallPeakUnit.pyx#L1443" rel="nofollow noreferrer"><code>call_broadpeaks</code></a>:</p>
  
  <p>The description attached to the function says:</p>
  
  <blockquote>
    <p>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.</p>
    
    <p>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']</p>
    
    <p>lvl1_cutoff_s:  list of cutoffs
    at highly enriched regions, corresponding to scoring functions.</p>
    
    <p>lvl2_cutoff_s:  list of cutoffs at less enriched regions,
    corresponding to scoring functions.</p>
    
    <p>min_length :minimum peak length,
    default 200.</p>
    
    <p>lvl1_max_gap   :  maximum gap to merge nearby enriched
    peaks, default 50.</p>
    
    <p>lvl2_max_gap   :  maximum length of linkage
    regions, default 400.</p>
    
    <p>Return both general PeakIO object for highly
    enriched regions and gapped broad regions in BroadPeakIO.</p>
  </blockquote>
  
  <p>To give some basic explanation, the algorithm (briefly) appears to be as follows:</p>
  
  <ol>
  <li><p>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 <code>-p</code> and level 2 is controlled by <code>--broad-cutoff</code>. When each peakset is called, they are immediately linked by the max gap parameter for each set.</p></li>
  <li><p>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.</p></li>
  </ol>
  
  <p>...</p>
  
  <p>This has a few implications:</p>
  
  <ol>
  <li><p>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).</p></li>
  <li><p>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 <code>--broad-cutoff 0.1</code> in broad peak mode, and the <code>-p 0.1</code> for narrow peak mode)</p></li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-21T12:09:59.927Z
Answers.Text: Could you clarify in your answer that "higher pval" is actually "lower" or "more significant" p value and that "lower pval" is actually "higher" or "less significant" p value? I foresee someone getting confused by the wording.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-22T03:08:50.783Z
Answers.Text: Added a clarification - thanks for the catch!
Answers.UserDisplayName: Daniel Kim
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-20T21:38:54.743Z
Answers.LastActivityDate: 2017-05-24T10:01:42.090Z
Answers.LastEditDate: 2017-05-24T10:01:42.090Z
Answers.UserDisplayName: Daniel Kim
Body: <p>The peak calling tool MACS2 can call peaks in either narrow peak mode (for focused signals like transcription factor ChIPseq) or broad peak mode (for more defuse signals, like certain histone modifications). </p>
  
  <p>The algorithm for narrow peak calling is well described in the MACS publication. But I don't find much documentation for how peak calling is different in broad peak mode. The manual only contains the following:</p>
  
  <blockquote>
    <p>--broad</p>
    
    <p>When this flag is on, MACS will try to composite broad regions in
    BED12 ( a gene-model-like format ) by putting nearby highly enriched
    regions into a broad region with loose cutoff. The broad region is
    controlled by another cutoff through --broad-cutoff. The maximum
    length of broad region length is 4 times of d from MACS</p>
  </blockquote>
  
  <p>But this doesn't really describe exactly how this is performed. </p>
  
  <p>So what is the algorithm MACS uses for calling broad peaks?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-20T19:36:54.420Z
Comments.Text: I found the [algorithm](https://github.com/taoliu/MACS/blob/24a1eab9fe7e885c27a37dbac2efb99d6da8dc74/MACS2/IO/BedGraph.pyx#L594), but honestly I could not really get the try/except construction on line 629...
Comments.UserDisplayName: Kamil S Jaron
Tags: <algorithms><chip-seq>
Title: How are MACS2's narrow peak and broad peak algorithms different?
ThreadText: # How are MACS2's narrow peak and broad peak algorithms different?
  
  (Asked by: username_0 on 2017-05-20)
  
  The peak calling tool MACS2 can call peaks in either narrow peak mode (for focused signals like transcription factor ChIPseq) or broad peak mode (for more defuse signals, like certain histone modifications).
  
  The algorithm for narrow peak calling is well described in the MACS publication. But I don't find much documentation for how peak calling is different in broad peak mode. The manual only contains the following:
  
  > --broad
  >
  > When this flag is on, MACS will try to composite broad regions in
  > BED12 ( a gene-model-like format ) by putting nearby highly enriched
  > regions into a broad region with loose cutoff. The broad region is
  > controlled by another cutoff through --broad-cutoff. The maximum
  > length of broad region length is 4 times of d from MACS
  
  But this doesn't really describe exactly how this is performed.
  
  So what is the algorithm MACS uses for calling broad peaks?
  
  **Tags:** algorithms, chip-seq
  
  > **Comments:**
  > **username_1**: I found the [algorithm](https://github.com/taoliu/MACS/blob/24a1eab9fe7e885c27a37dbac2efb99d6da8dc74/MACS2/IO/BedGraph.pyx#L594), but honestly I could not really get the try/except construction on line 629...  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2017-05-20. Score: 10 | Accepted answer
  
  The key function is [`call_broadpeaks`](https://github.com/taoliu/MACS/blob/master/MACS2/IO/CallPeakUnit.pyx#L1443):
  
  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:
  
  1. 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.
  2. 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:
  
  1. 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).
  2. 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)
  
  > **Comments:**
  > **username_3**: Could you clarify in your answer that "higher pval" is actually "lower" or "more significant" p value and that "lower pval" is actually "higher" or "less significant" p value? I foresee someone getting confused by the wording.  
  > **username_2**: Added a clarification - thanks for the catch!  
  
  ---
  
  