How peaks are calculated ?

[biocyberman]

I noticed that after callpeak command with treat vs control MeDIP samples, I don't have peaks with fold-changes < 1, which, I suppose, are demethylated spots. Is it correct to expect fold-change < 1 for demethylated peaks? In general, a plain English explanation of how peaks and fold-changes are calculated will be very helpful.

I tried looking for that information in the publication http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3120977 but could not find anything. The last resort could be to look into the source code, but I prefer to have the authors' words first.

[taoliu]

It’s impossible to call ‘depleted’ regions from ‘callpeak’ command since it is designed to find ‘enriched’ regions only. If you have MeDIP data from two conditions, you can consider ‘bdgdiff’ command after you get four bdg files from two separate runs of ‘callpeak’ on both conditions.

To analyze MeDIP data, you may consider more specific tools such as MEDIPS or MnM.

Tao

On Jun 23, 2014, at 7:28 AM, biocyberman notifications@github.com wrote:

I noticed that after callpeak command with treat vs control MeDIP samples, I don't have peaks with fold-changes < 1, which, I suppose, are demethylated spots. Is it correct to expect fold-change < 1 for demethylated peaks? In general, a plain English explanation of how peaks and fold-changes are calculated will be very helpful.

I tried looking for that information in the publication http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3120977 but could not find anything. The last resort could be to look into the source code, but I prefer to have the authors' words first.

—
Reply to this email directly or view it on GitHub.

[biocyberman]

Hello Tao,
Thank you for very helpful hints. I will take a closer look at MEDIPS and MnM. However, just to think about a way (around) to use MACS for 'depleted' peaks also, would it still be OK to switch 'treated' and 'control' samples to call these peaks? That means, I first run

macs2 callpeak -t treat_sample -c control_sample

and then run macs2 callpeak -t control_sample -c treat_sample

Afterward I would do some preprocessing of peak_names, fold_changes before JOINING the two peak_list together so I will have both 'enriched' and 'depleted' peaks.

My concern is that it may not be correct to join the two output together because may be normalized differently. What is your thoughts about this?

Thanks
Vang

[taoliu]

Sorry for late reply! I am cleaning issues nowadays and just saw your question...

Yes. You can do so. By default, MACS2 uses size of smaller dataset to normalize, so in your case, this step is fine. But MACS2 will use control data to calculate so-called local background in 1,000bps, 10,000 bps surrounding regions. This may not be a good solution for your methylation data since you may want to compare two libraries at the same position. You can try to disable this by --slocal 0 --llocal 0.

I close this issue for now, since it's more about using MACS2 not about coding. If you still have question, please send email to MACS google user group or email me directly.

-T