Extract LINE-1 methylation status from WGBS data

Question

I've been trying to figure out the correct way of extracting methylation status of particular elements from WGBS data, and found several possible approaches, but would like to discuss which one could be the best:

1. using software dedicated to the extraction of these elements on WGS (like TraFiC or MELT and then quantify methylation (like Bismark)
2. using bedtools to scan bam file with a filtered version of UCSC RepeatMasker table in bed format) and then quantify methylation
3. I've seen the use of custom perl scripts on fastq files like in this paper, in which the Supplemental Methods Section states: "A custom Perl script was written to assay methylation of CpGs within a previously published motif of LINE-1 elements (Yang et al., 2004). Methylation levels for each of the four CpGs within the consensus LINE-1 element were determined globally for each WGBS sample prior to alignment using Line1_FASTQ.pl. Dup15q and control were compared and tested for significance using an unpaired t-test."

Any thoughts on this?

Answer 1

Have you checked out WGBS Tools on github? All of the custom scripts used in the paper that you referenced are in that repository. Each tool has a brief description of how to use it.

Answer 2

After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:

element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.

Answer 3

Via BEDOPS:

$ bedmap —-echo —-echo-map-id-uniq <(bam2bed < reads.bam) <(rmsk2bed < repeats.rmsk) > answer.bed