MACS2 coverage calculations

Question

[ariel.balter]

This relates to MACS2 Narrow peak (narrowPeak) output format

I'm trying to understand the output macs2 callpeak gives. As you can see, there are some identical coordinate ranges for which macs gives multiple, different fold changes, p-values, and q-values. However, when I calculate their coverage using bedtools multicov I get the exact same coverage (as I would expect). What calculation is macs2 giving me?

$ bedtools multicov -bams ../bam/mock_lsd1_A.bam ../bam/mock_input_A.bam -bed mock_lsd1_A_peaks.narrowPeak > mockAcoverage.bed
$ head -20 mockAcoverage.bed
1       10019   10462   mock_lsd1_A_peak_1a     574     .       5.29438 61.71617        57.47774        69      387     74
1       10019   10462   mock_lsd1_A_peak_1b     542     .       7.67026 58.48692        54.29382        378     387     74
1       11195   11356   mock_lsd1_A_peak_2      27      .       3.65672 5.12036 2.72495 36      24      4
1       46749   46984   mock_lsd1_A_peak_3      232     .       9.81878 26.77591        23.23711        129     72      7
1       180749  180991  mock_lsd1_A_peak_4      132     .       4.15265 16.42585        13.24979        168     101     20
1       184529  184678  mock_lsd1_A_peak_5      23      .       2.92799 4.70182 2.33619 42      42      16
1       359363  359549  mock_lsd1_A_peak_6      33      .       3.21866 5.85335 3.36933 78      37      6
1       450836  450954  mock_lsd1_A_peak_7      17      .       3.29258 3.96944 1.70144 46      24      4
1       588135  588463  mock_lsd1_A_peak_8a     122     .       7.35857 15.38114        12.25107        60      58      7
1       588135  588463  mock_lsd1_A_peak_8b     43      .       4.37685 6.99404 4.35689 257     58      7
1       590976  591412  mock_lsd1_A_peak_9a     37      .       4.24169 6.32635 3.77917 75      113     13
1       590976  591412  mock_lsd1_A_peak_9b     146     .       6.47838 17.83555        14.60371        205     113     13
1       590976  591412  mock_lsd1_A_peak_9c     245     .       10.42606        28.15508        24.58193        360     113     13
1       685661  685809  mock_lsd1_A_peak_10     19      .       2.63428 4.24476 1.93949 82      31      10
1       778587  778870  mock_lsd1_A_peak_11     23      .       3.56365 4.67361 2.31693 105     40      6
1       779189  779324  mock_lsd1_A_peak_12     30      .       3.51999 5.46828 3.00177 61      25      4
1       804698  804993  mock_lsd1_A_peak_13a    33      .       3.88036 5.80592 3.32561 49      45      8
1       804698  804993  mock_lsd1_A_peak_13b    46      .       4.59437 7.28475 4.63472 224     45      8
1       812009  812338  mock_lsd1_A_peak_14     105     .       6.13598 13.58947        10.53462        94      71      13
1       812860  812965  mock_lsd1_A_peak_15     18      .       2.99039 4.11731 1.83291 90      25      6

[ejm32]

The the fold change is calcuated based on the summit, not necessarily the entire width of the peak. Furthermore, the fold change calculated by macs2 is not a simple division of IP/input for a given peak. The denominator (and possibly, the numerator) is a function of the lambda from the peak calling and so looks at a region that is larger than the peak to calculate the "input"

From the macs2 github doc:

fold enrichment for this peak summit against random Poisson distribution with local lambda, -log10(qvalue) at peak summit