bdgpeakcall vs callpeak

[Aarthi Mohan]

Hi All,
I have called peaks using MACS2 callpeak and got the peak file and SPMR normalized bedgraphs. I then used this bedgraph files to do peak calling with bdgpeakcall (after doing the bdgcmp step with ppois) command withe same cutoff I used in callpeak (0.05).But the number of peaks reported by bdgpeakcall is lesser than callpeak by 1000. Has anyone encountered this?

I have a paired-end ChIP-seq of RNAPol2. I am desperately trying to use the bdgpeakcall, so that I can disable the global background from lambda. So, I wanted to test if callpeak and bdgpeakcall gives same results, which is not the case from my files.

Any help will be greatly appreciated.

Thanks,
Aarthi

[Timothy Parnell]

Since you’re using SPMR bedgraph files, make sure you’re using the -S option in the bdgcmp module to scale the reads back to normal depth (not per million). If you don’t, then the p and q values are diminished. Alternatively, run the callpeak module without the --SPMR option to get just the scaled (relative to treatment and control) bedgraphs.

I don’t believe global background lambda is calculated unless you give the --nolambda option in callpeak. Callpeak always prefers local lambda, and that itself can be controlled with the --slocal and --llocal options. The bdgpeakcall module has no options for controlling background because the control lambda track used in bdgcmp is itself the calculated background. Use a control ChIP pileup if you truly want to disregard all calculated background (not entirely recommended for most scenarios).

Tim

> --
> You received this message because you are subscribed to the Google Groups "MACS announcement" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to macs-announcem...@googlegroups.com.
> To post to this group, send email to macs-ann...@googlegroups.com.
> Visit this group at https://groups.google.com/group/macs-announcement.
> For more options, visit https://groups.google.com/d/optout.

[Aarthi Mohan]

Hi Timothy,

Thanks for the reply. The default behaviour of callpeak is to use the maximum value from the regions of length d, 1kb, 10kb, and whole genome sequencing depth. That's why the control pileup never has a value '0' (from this tutorial https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands#Step_2_Decide_the_fragment_length_d). 

And yes, If i use the bedgraph file generated without --SPMR, number of peaks are consistent between callpeak and bdgpeakcall. 

It will be nice if Tao can clarify on the following,

What will be the steps to generate the local lambda in case of a paired-end sequencing? Especially, what will be the 'fragment size' (in case of paired-end input) used to generate the 'd background' as done in this tutorial (https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands#The_d_background).

Thanks,

Aarthi