Q: ChiP-seq of histone PTMs: Control Peaks not Aligned to Treatment Peaks

[capslockwizard]

Use case I am processing ChiP-seq data of histone PTMs from yeast. The data comes with the input dataset for control. I have performed performed peak calling by executing MACS2 (v2.2.8) with the following command:

macs2 callpeak -B -t treatment.bam -c control.bam -g 12157105 --keep-dup all -n output -q 0.05 --nomodel --shift 67 --extsize 43 --min-length 100 --max-gap 73

The cross-correlation analysis shows the maximum peak value of 134 and the size of the tags is 60. Thus I decided to shift by 134/2=67 (rounded down). Assuming ideal nucleosome size of 147 and the middle of my tags are located at the peak, I extend the tags by 147/2 - (60/2) = 43 (rounded down).

Describe the problem The pileup of the treatment data looks reasonable, however, the peaks of the background lambda derived from the control data are not aligned to the peaks from the treatment data: The first track is my treatment and the second is the background lambda

Describe the solution you tried After going through the MACS2 code, I find that the control data is not shifted at all: https://github.com/macs3-project/MACS/blob/f350e8c3a7a4e5aa890dd31ad73cc7a4472e1d66/MACS2/IO/CallPeakUnit.pyx#L524-L526

I modified the code to become:

            ctrl_pv = self.ctrl.pileup_a_chromosome( chrom, self.ctrl_d_s, self.ctrl_scaling_factor_s,
                                                     baseline_value = self.lambda_bg,
                                                     directional = False,
                                                     end_shift = self.end_shift )

And recompiled MACS2 and reran MASC2 on the data, this time I get: The first track is my treatment, second track is the previous background lambda and the third track is the latest background lambda

I am wondering if what I did here is reasonable?