I have performed peak calling on a number of separate ChIP-seq experiments and would like to harmonise the peaks from each of the experiments in order to convert my data into a matrix for further analysis. I would like to do this by taking the union of the called peaks (available in BED / narrowPeak / MACS2 format) and retrieving the enrichment score / q-value for all of the experiments on each peak. Is this possible without writing my own extension to MACS2?
You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like:
macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7
You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise.
But you can then use these "non-significant" values in combination with a
bedtools map command, to map them to your predefined union of "actual" peaks.