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:
- get gene coordinates from UCSC
- create bins using bedtools makewindows
- count reads within each window of celltyp a and b using bedtools multicov
- divied readcounts from celltype a and b by the total number of reads to normalize for libary size
- associate the windows with the corresponding gene (i think should be possible with the UCSC bed file)
- 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