I'm trying to perform a differential binding ChIP-seq experiment and am struggling with the best way to incorporate my replicates using MACS2. I have two sets of samples, control rep 1-3 and control input, as well as treated rep 1-3 and treated input. If I run macs2 callpeak as recommended I run:
macs2 callpeak -B -t ctrl.rep1.bam ctrl.rep2.bam ctrl.rep3.bam -c ctrl.input.bam -f BAM -g hs -n ctrl.repAll.peaks --nomodel --extsize 146 --buffer-size 1000000 macs2 callpeak -B -t treat.rep1.bam treat.rep2.bam treat.rep3.bam -c treat.input.bam -f BAM -g hs -n treat.repAll.peaks --nomodel --extsize 146 --buffer-size 1000000
This will generate output files in the bedgraph format that I can then use to run macs2 bdgdiff:
macs2 bdgdiff --t1 treat.repAll_treat_pileup.bdg --t2 ctrl.repAll_treat_pileup.bdg --c1 treat.repAll_control_lambda.bdg --c2 ctrl.repAll_control_lambda.bdg -l 146 -d1 spikeIn_treat -d2 spikeIn_ctrl --o-prefix treat_vs_ctrl
However, these bedgraph files are the POOLED replicates. I want to be able to utilize my individual replicates in order to get some statistical power.
1) Does anyone have experience running macs2 bdgdiff across multiple replicates for treated/untreated samples?
2) How would I combine the log odds ratio in the macs2 bdgdiff output files to get a meaningful measurement of the difference in binding pre- and post-treatment?
I've been searching online and can't find anything. Any help would be appreciated.