The q-value in broadPeak output does not reflect the "--broad-cutoff " setting?

[aliceZhu]

Hello @taoliu ,

While running macs2 on broad mark mode, I set the q-value cutoff as 0.1 for the broad region, and 0.05 for the narrow region, by:

# Command line: callpeak -t treatment.bam -c control.bam --broad -g hs --broad-cutoff 0.1 -n H3K27me3_vs_input --outdir output/
# ARGUMENTS LIST:
# name = H3K27me3_vs_input
# format = AUTO
# ChIP-seq file = ['treatment.bam]
# control file = ['control.bam']
# effective genome size = 2.70e+09
# band width = 300
# model fold = [5, 50]
# qvalue cutoff for narrow/strong regions = 5.00e-02
# qvalue cutoff for broad/weak regions = 1.00e-01
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 1000 bps and 10000 bps
# Broad region calling is on

I would thus expect the column 9 (9th: -log10qvalue) in the output broadPeak file to be all >= 1. However, I still see values as low as 0.5 in the 9th column, suggesting the filtering does not occur.

Is that a bug or is the q-value somehow postoperatively re-computed once after calling peaks using 0.1 for broad, and 0.05 for narrow? How should I interpret the 9th column, if I want to get all peaks with q-value<=0.01 ?

When I change the value of n in "--broad-cutoff n" , the number of peaks outputted does change.

Thanks a lot for your time and help!

[Zepeng-Mu]

I noticed the same thing. Basically change qvalue cutoff does not alter the number of peaks in output. I'm wondering how to solve this?

[taoliu]

@Zepeng-Mu Please use --broad-cutoff to control the behavior in broad mode.

$ macs2 callpeak -h
...
  --broad               If set, MACS will try to call broad peaks using the
                        --broad-cutoff setting. Please tweak '--broad-cutoff'
                        setting to control the peak calling behavior. At the
                        meantime, either -q or -p cutoff will be used to
                        define regions with 'stronger enrichment' inside of
                        broad peaks. The maximum gap is expanded to 4 * MAXGAP
                        (--max-gap parameter). As a result, MACS will output a
                        'gappedPeak' and a 'broadPeak' file instead of
                        'narrowPeak' file. Note, a broad peak will be reported
                        even if there is no 'stronger enrichment' inside.
                        DEFAULT: False

  --broad-cutoff BROADCUTOFF
                        Cutoff for broad region. This option is not available
                        unless --broad is set. If -p is set, this is a pvalue
                        cutoff, otherwise, it's a qvalue cutoff. Please note
                        that in broad peakcalling mode, MACS2 uses this
                        setting to control the overall peak calling behavior,
                        then uses -q or -p setting to define regions inside
                        broad region as 'stronger' enrichment. DEFAULT: 0.1

[Zepeng-Mu]

I just figured this out, thanks!

[chen5027]

Hello, all, I need some help to optimize my broad peak calling please - I am trying to call CUT&Tag peaks for H3K27me3 but keep seeing narrow peaks are called instead of broad peaks/domains, which are typical for H3K27me3. On tracks, it is clear to me that there are broad peaks but they were just called as broken narrow peaks. Could you suggest what I can do please?

Thanks!

Dahong

macs2 callpeak --broad -t file.bam -f BAMPE -g dm -n H3K27me3_broad -B -q 0.1 --broad-cutoff 0.1 --max-gap 80 --outdir ../MACS2/ &> ../MACS2/logs/log.txt

[taoliu]

@chen5027 My suggestion is to use p-value cutoff rather than q-value cutoff since the p-value is less strigent, and considering H3k27me3 is a fairly weak mark. Also, you may visualize the signal track in the genome browser to see if you can observe the 'peak'.