Friday, March 1, 2013

Using ENCODE methylation data (RRBS) in R

ENCODE project has generated reduced-representation bilsulfite sequencing data for multiple cell lines. The data is organized in an extended bed format with additional columns denoting % methylation and coverage per base. Luckily, this sort of generic % methylation information can be read in by R package methylKit, which is a package for analyzing basepair resolution 5mC and 5hmC data.

The code snippets below show how to read RRBS bed file produced by ENCODE. But, let's first download the files.
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/wgEncodeHaibMethylRrbsA549Dm002p7dHaibSitesRep1.bed.gz

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/wgEncodeHaibMethylRrbsAg04449UwstamgrowprotSitesRep1.bed9.gz

Unfortunately, methylKit currently can not read them directly because the track definition line causes a problem. It should be deleted from each bed file. Ideally, methylKit should be able to skip over lines (this issue should be fixed in later versions)
For now, we have to use some unix tools to remove the first lines from the bed files. You run the code below in your terminal. This set of commands will delete the first line from every *.gz file in the directory so be careful.
for files in *.gz
do
  gzip -dc "$files" | tail +2 | gzip -c > "$files".tmp
  if [ "$?" -eq 0 ]; then
    mv "$files".tmp "$files"
  fi
done

Now we can read the files using methylKit. The pipeline argument defines which columns in the file are corresponding to chr,start,end,strand, percent methylation and coverage:
You can also read multiple files at a time:

Since we have read the files and now they are methylKit objects, we can use all of the methylKit functionality on these objects. For example, the code below plots the distribution of methylation % for covered bases.

getMethylationStats(obj[[1]], plot = TRUE)


You can check the methylKit vignette and the website for the rest of the functionality and details.




7 comments:

  1. I tried the above. But I get the following error message

    > library(methylKit)
    > obj = read("wgEncodeHaibMethylRrbsA549Dm002p7dHaibSitesRep1.bed.gz",
    + sample.id = "test",assembly = "hg19",header = FALSE,
    + context = "CpG",resolution ="base",
    + pipeline = list(fraction = FALSE,chr.col = 1, start.col = 3, end.col = 3,
    + coverage.col = 5,strand.col = 6, freqC.col = 11))
    Error in .local(location, sample.id, assembly, pipeline, header, context, :
    unknown 'pipeline' argument, supported alignment pipelines: amp
    In addition: Warning message:
    In if (pipeline %in% c("amp", "bismark")) { :
    the condition has length > 1 and only the first element will be used
    >

    ReplyDelete
    Replies
    1. Try the latest version of the package with the appropriate version of R.

      Delete
  2. Thanks for the post. Maybe it's better to use coverage.col = 10, which is the reads count or coverage for given site. The 5th column is the score of the track, which is equal to coverage when coverage is not greater than 1000.

    ReplyDelete
  3. Great post! Thank you so much for sharing..

    For those who want to learn R Programming, here is a great new course on youtube for beginners and Data Science aspirants. The content is great and the videos are short and crisp. New ones are getting added, so I suggest to subscribe.
    https://www.youtube.com/watch?v=BGWVASxyow8&list=PLFAYD0dt5xCzTQHDhMPZwBoaAXWeVhZzg&index=19

    ReplyDelete
  4. Hi my friend, Can you show me how to insert R or Perl Script into the blogs. Thanks

    ReplyDelete
  5. Hi my friend, Can you show me how to insert R or Perl Script into the blogs. Thanks

    ReplyDelete