The problem: calling peaks with MACS2 from ChIP-seq data
What did it happen? See below.
- I got bam files after mapping to the reference genome, then removed duplicates using samtools with 'rmdup'parameter.
Then I got bigwiggle files using the command below.
bamCoverage -b test.uq.st.rmdup.bam -o test.uq.st.rmdup.bw --smoothLength=200 --ignoreForNormalization chrX chrM --normalizeTo1x 2451960000
- I uploaded the .bw file into the Genome Browser and I saw four peaks following HBE1 gene as showed in the attached picture.
But when I called peaks with the command below,
macs2 callpeak -t test.uq.st.rmdup.bam -f BAMPE -g hs -n test --outdir ./and I didn't call the four peaks out.
But when I set p-value as 0.0001 instead of default q-value0.05 with the command below,
macs2 callpeak -t test.uq.st.rmdup.bam -f BAMPE -g hs -p 0.0001 -n test --outdir ./
only the first peak was called out and the rest of the three peaks were missed. It seems weird.
What am I expecting? MACS2 can call all the four peaks out which is consistent with results the Big wiggle track shows.
Actually, the fold_enrichment from big wiggle file for the peak2 is very high compared to the peak1. I don't know why it does not call the second one out. Is it because the peak1 meets Poisson Distribution and the others do not? Could you please help me with this issue and guide me on how to figure out this problem? I would greatly appreciate. Thanks.