Friday, October 5, 2012

How to read BSMAP methylation ratio files into R via methylKit

BSMAP is an aligner for bisulfite sequencing reads. It outputs aligned reads as well as methylation ratios per base (via script). The methylation ratios can be read into R via methylKit package and regular methylKit analysis can be performed using the BSMAP data.

First,  you need to separate BSMAP methylation ratio file based on methylation context. Cs in CpG, CHH and CHG context should be separated into different files. This can be achieved with a perl or awk one-liner. If you want to analyze Cs from all contexts then of course you don't need to separate the file. For example, to get all Cs in the CpG context you will need to do something like following:

awk '($3=="-" && $4~/^.{1}CG/ ) || ($3=="+" &&  $4~/^.{2}CG/)' BSMAPexample.txt > CpG.txt

Here is how methylation ratio file from BSMAP looks like (file includes only Cs in CpG context):

chr  pos     strand context ratio total_C methy_C CI_lower CI_upper
chr1 3121589 +      CGCGT   0.000 56      0       0.000    0.064
chr1 3121597 +      ATCGG   0.000 56      0       0.000    0.064
chr1 3121599 +      GTCGT   0.000 56      0       0.000    0.064
chr1 3121605 +      CTCGG   0.000 56      0       0.000    0.064
chr1 3121606 +      TGCGC   0.000 56      0       0.000    0.064
chr1 3121607 +      GGCGC   0.000 56      0       0.000    0.064
chr1 3121611 +      CTCGA   0.000 56      0       0.000    0.064
chr1 3121614 +      TACGC   0.000 56      0       0.000    0.064
chr1 3121631 +      CTCGT   0.000 56      0       0.000    0.064

You can read the methylation ratio file by using "pipeline" argument in read() function. You need to provide a list of column numbers corresponding to chr,start,end,strand,coverage and ratio of methylation. Actually, you can read any generic methylation ratio or percentage file using this method. The file needs to have the location information (chr,start,end and strand), coverage information and methylation percentage or ratio information.

Thanks to Maxime Caron for sharing the BSMAP methylation ratio file.

Reference for methylKit:
Altuna Akalin, Matthias Kormaksson, Sheng Li, Francine E. Garrett-Bakelman, Maria E. Figueroa, Ari Melnick, Christopher E. Mason.(2012)"methylKit: A comprehensive R package for the analysis of genome-wide DNA methylation profiles." Genome Biology , 13:R87.

PS: For the newer releases of BSMAP (~v2.73), you may need to coerce the CT counts to an integer. BSMAP ratio extractor may output effective CT counts as float but they should be coerced to integer before they are read in by methylKit. methylKit (>= v0.5.6) will automatically coerce those values to nearest integer. However, you can also use awk to coerce float to integer values. The example below coerces the 6th column of the the file to integer

awk '{OFS="\t";print $1,$2,$3,$4,$5,int($6);}' sample.BSMAP.txt


  1. Dear Altuna,

    Referred to this sentence:
    "you need to separate BSMAP methylation ratio file based on methylation context. Cs in CpG, CHH and CHG context should be separated into different files. This can be achieved with a perl or awk one-liner."

    Could you please write the Perl or awk one liner example, so that we can use it for the task you refer to?

    With best wishes

  2. You can use something like this to get Cs in the CpG context
    awk '($3=="-" && $4~/^.{1}CG/ ) || ($3=="+" && $4~/^.{2}CG/)' BSMAPexample.txt > CpG.txt

    for CHH and CHG context you will probably need to write a slightly more complex regular expression.

  3. There's another problem I've come across in this workflow. Once I have the BSMPA2.7 results in bsp format, I get the following error while running

    Traceback (most recent call last):
    File "", line 81, in
    meth[cr] = array.array('I', [0]) * len(ref[cr])

    To my understanding (I don't know any Python yet) it's due to a Python module but I have no clue which module is involved.
    Has anybody else experienced this error using

    Thanks in advance you for your attention

    1. Hi Jose Luis,
      I think you should take this to BSMAP developers or their e-mail group, they should be able to help you on this.


  4. I think it should be awk '($3=="-" && $4~/^.{2}GC/ ) || ($3=="+" && $4~/^.{2}CG/)' BSMAPexample.txt > CpG.txt ( for minus strand it should be ($3=="-" && $4~/^.{2}GC/ ))

  5. I think not, BSMAP readme.txt documents says that the strand of the 5bp sequence is the plus strand irrespective of the strand of the covered CpG . So a CpG on the minus strand will appear as ".{1}CG". For example, If the strand="-" the following string denotes a CpG on the minus strand: ACGAA. The minus strand of the same string will appear as: TTCGT

  6. Hi, Thanks for the info on reading in the BSMAP files. It was really helpful. I

    was wondering if you have any advice for generating a methylRawList from individual files that were read in this way.

    The example given in the read{methylKit} help file only shows how to generate a methylRaw object from 'generic' files, I can't figure out how to rework the example to generate a methylRawList from these individual objects.


    1. you should be able to give a list of file locations, keep the "pipeline" argument as shown above, give "" as a list, also use "treatment" argument, And that will read a list of files and produce a methylRawList

      some thing like the following (file.list is a list of file locations)
      myobj=read( file.list,pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,
      coverage.col=6,strand.col=3,freqC.col=5 ),"test1","test2","ctrl1","ctrl2"),assembly="hg18",treatment=c(1,1,0,0))

    2. Wonderful! Thank you

  7. Hello,

    Firstly, thanks for posting details on how to use methylKit on BSMAP.

    I have recently installed BSMAPv2.74 and the program outputted a file with the following headers:

    context ratio

    I can't find documentation that will help me relate these 11 headers to the 9 headers you have written about in your initial post. Some of them are of course obvious, but I'm quite confused trying to link 'total_C' and 'methy_C' with 'eff_CT_count', 'C_count' and 'CT_count'.

    Do you know the answer, or could you post a link to where this might be described?

    Kind regards

    1. Hi,
      I don't have the most recent version of the BSMAP but I think they usually describe the output in a readme file or in the script help. I think eff_CT_count or CT_count is now replaced the total_C. The best I can suggest is to find the readme or help page where the output is described and/or take it to the authors of BSMAP. and I would appreciate it if you can post the answer as a comment here.


    2. I've contacted the BSMAP team and will respond as soon as i get an answer.

    3. Hi,

      So i contacted the authors. They said that the eff_CT_count is new and is calculated by adjusting the methylation ratio for the C/T SNP, using the reverse strand mapping information. It's also possible to turn off these columns when you use by using this parameter '-i no-action'.

      Therefore, the 'C_count' and 'CT_count' columns correspond to the 'total_C' and 'methy_C' column headers respectfully.

    4. Thanks for the update Claire!!

  8. This comment has been removed by the author.

  9. 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.

  10. This comment has been removed by the author.

  11. This comment has been removed by the author.

  12. Why is the position of the CG important? And why is the position of the CG different on the - vs + strand? (^.{2}CG vs ^.{1}CG).

    I'm trying to calculate Cs in CpG, CHH and CHG context and it would be nice to know what is going on.

    1. Because the position will change depending upon the strand which the cytosine is present.

      strand context type
      + CTCGT CGH
      - TCGTT CGH
      - GTGCT CHH
      - CGGTT CHG

      The context is printed in relation to only the + strand.
      This regular expression will only find cytosines in the CG context. When the cytosine is on the + strand it is in the 3rd position. When the cytosine is on the - strand there should also be a another cytosine in the 2nd position (complementary to the G on the - strand).

      It matters that you not just specify 'CG' as a CG can be in a different location for some lines (see the 4th line in my example, a CG is present but it is a CHG context).

  13. For those interested in what and awk statement would look like for all C contexts:

    awk '(NR>1){if(($3=="-" && $4~/^.CG../ ) || ($3=="+" && $4~/^..CG./)) print $1"\t"$2-1"\t"$2"\t"$3"\t""CG""\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12; else if(($3=="-" && $4~/^C[AGT]G../ ) || ($3=="+" && $4~/^..C[ACT]G/)) print $1"\t"$2-1"\t"$2"\t"$3"\t""CHG""\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12; else if(($3=="-" && $4~/^[AGT][AGT]G../ ) || ($3=="+" && $4~/^..C[ACT][ACT]/)) print $1"\t"$2-1"\t"$2"\t"$3"\t""CHH""\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12; else print $1"\t"$2-1"\t"$2"\t"$3"\t""CNN""\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12}' BSMAPratio/infile.txt > BSMAPratio/BSMAP_outfile.txt

    1. I should note that a column was added for start and stop information ($2 and $2-1) this is in case you want to visualize the data downstream as BED format file require start and stop information.