macs2 bdgdiff calculation without genomic background control?

Question

[biocyberman]

We do not have genomic background control. We ran treated (LT1) samples and non-treated (control, LC1) samples. My task is to identify regions that are differentially methylated between control and treated. The commands adapted from here. However, this case we do not have the genomic background control so how does macs2 generate the b3tTime0_control_lambda.bdg? And I am doing it right in order to get the differently methylated areas?

macs2 callpeak -B -t LC3.bam -n b3cTime0 -f BAM -g 2105347242 --outdir b3
macs2 callpeak -B -t LT1.bam -n b3tTime0 -f BAM -g 2105347242 --outdir b3

macs2 bdgdiff --t1 b3/b3tTime0_treat_pileup.bdg\
 --c1 b3/b3tTime0_control_lambda.bdg\
 --t2 b3/b3cTime0_treat_pileup.bdg\
 --c2 b3/b3cTime0_control_lambda.bdg\
 --d1 13141281 --d2 17278204 -g 60 -l 120 --o-prefix diff_

[ablanchetcohen]

It may just be me, but I don't understand your post.

How can you use MACS2 to determine methylation? MACS2 is used to find enriched regions (peaks) in immunoprecipitated DNA, and does not check the methylation status of the bases, as far as I know.

To compare the immunoprecipitated regions in two ChIP-Seq experiments, without any background controls, I just run MACS2 callpeaks on both samples. I don't run MACS2 bdgdiff. So, just run MACS2 callpeak on both samples at the same time. You'll have to run it twice, inverting the treatment and the control, to get peaks enriched in both conditions.

But, this has nothing to do with methylation.

[biocyberman]

@ablanachetcohen My fault. I should mention that my libraries were generated with MeDIP-seq protocol.

[ablanchetcohen]

Interesting. No, I guess it's just my ignorance. We always use bisulfite conversion, occasionally combined with ChIP-Seq (BisChIP-Seq). I hadn't heard of MeDIP.

My answer though would still be the same. Just run MACS2 callpeak twice, putting each sample as either treatment or control. If the background is the same in both samples, this approach is valid. I sometimes use this approach even when I have background controls for the samples, because MACS2 bdgdiff returns only a log10 probability ratio that is difficult to interpret. MACS2 callpeak returns much more detailed results that make it easier to subsequently filter the peaks.

[biocyberman]

I actually did think about using this approach, then I saw the tutorial of bdgdiff. And I thought it may be the better tool to do the task. Anyway, thanks for the suggestion, I will try the "flipping peakcall" approach as well.

Answer 1

[mmarins]

Dear All,

I got interested in this questioning because I am also starting using MACS2 for peak calling and differential binding. I have three files file 1 is the control of cells treated with DMSO. File 2 is the cells tretaed with my drug. Both these samples were imunoprecipitated with the same transcription factor. ANd I have a third file that is INPUT DNA.

So please, in my case, using bdgdiff, I should replace in the command line cond1_CHIP.bam with my file containing the DMSO sample and CONtrol.bam for my INPUT DNA sample? And the same for the treatment with the drug I also use INPUT DNA sample as the control?

You suggest to use MACS2 only. But then how you get the common peaks and the unique peaks? You filter?

I apologize if this is too basic. Att, Mozart

[meer.k.mustafa]

If you want to find enriched peaks of drug treatment over the DMSO control then you would run a command like so:

macs2 bdgdiff --t1 TFdrug.bam --c1 TFdrug_input.bam --t2 TFDMSO.bam --c2 TFDMSO_input.bam --d1 --d2 -l <read lengths> --o-prefix TF

What you'll also want to do is specify --d1 & --d2 which is the reads filtered value from running macs2 predictd -i TFdrug.bam and macs2 predictd -i TFDMSO.bam and -l is your read lengths e.g. 36, 101, etc...

This macs2 bdgdiff should output 3 BEDs: one being t1 peaks enriched compared to t2, one being t2 peaks being enriched compared to t1, and one being peaks that are commonly enriched. You'll likely want to use the first BED file here.

Best of luck, Meer

Answer 2

[Vanilla]

For the question "how does macs2 generate the b3tTime0_control_lambda.bdg", I found the answer here: https://groups.google.com/forum/#!msg/macs-announcement/9N-sMFc8784/SxlSU8joSkwJ. In brief, MACS2 estimates control lambda from the ChIP sample itself when the control sample is absent.

Answer 3

[Stefan.Kurtenbach]

Try this differential peak caller I have written and used successfully. This first version is not very quick yet, so I recommend palatalization. Let me know what you think!

https://github.com/StefanKurtenbach/Ucallpeak