Tool to compare whole genebody enrichment between different samples (MeDip-seq)

Question

[florian.noack]

Hi everyone, I have some MeDip datasets and did already comparisions of my datasets based on peaksets (MACS2). However I would like to find out if there is a average increase/decrease in methylation of, for example, genebodys (and if yes which genes) between celltyp A and celltyp B (normalized by liberay size and genebody length). Somehow similare to what is done in this paper (http://hmg.oxfordjournals.org/content/23/3/657.long Figure 4C):

... then estimated the 5 hmC and 5 mC densities in 200 bp-sized bins within 5 kb of the TSS and within 5 kb of the transcription end site (TES) of each gene using the refGene annotation of UCSC Genome Bioinformatics (http://genome.ucsc.edu/). A normalized methylation enrichment score for 5 hmC or 5 mC was calculated from the number of mapped 5 hmC and 5 mC reads in each bin divided by the bin size (200 bp) and the total number of mapped sequence reads on the genome for each hES or NP cell multiplied by the human genome size (3 × 10^8 bp).

I have a rough idea how I can do this calculations using bedtools:

1. get gene coordinates from UCSC
2. create bins using bedtools makewindows
3. count reads within each window of celltyp a and b using bedtools multicov
4. divied readcounts from celltype a and b by the total number of reads to normalize for libary size
5. associate the windows with the corresponding gene (i think should be possible with the UCSC bed file)
6. calculate the average methylation of a gene by ave all windows falling into this gene (should be there a additional normalization for different length of genes ?)

However Iam not a informatician so I would prefer to use a premade tool or r.script to do this kind of calculation to avoid doing mistakes.

Thanks for any suggestion, Flo

Answer 1

[Devon Ryan]

You can do all of this except the per-gene averaging with deepTools, we even have a public Galaxy server so you don't have to install anything. For the averaging, you can load stuff into R if you really want to and average that way, though we find it more informative to just plot the results (e.g., with plotProfile).

[florian.noack]

Just tried deepTools and its a great tool, thanks a lot ! However one question. To generate bigwig files i use bamcompare with -b1 ab treated sample -b2 input control --ration substraction and --normalizeUsingRPKM . I do this for celltype A and celltype B. Iam now confused if the bigwig files of celltype A and celltype B are normalized to each other as well (RPKM) and i can visually compare them in ucsc ? I ask since I have not only celltyp A and B but also C, D, E ... ;) so i can not simply use bamcompare with A and B.

Or in this case its better to use bamcoverage to generate bigwig files since I know (because of other tools) that the enrichment worked quite well.

[Devon Ryan]

I'd normally use bamCoverage in cases like this.

[florian.noack]

Great I will check this. However a "simple" plot will not help that much since I would like to correlate the changes of average DNA methylation of the gene body with the gene expression of the corresponding gene (like they did in the paper mention above).