Create BSseq object for BSmooth

[ge2sasag]

Hello,
I got whole genome bisulfite sequencing methylation data as bed files (output from Bowtie2) with 3 columns for genomic position and M (methylation) and Cov (Coverage) columns. I´m trying to create the BSseq object needed to run BSmooth pipeline but I´m not able to do it. I´m very newbie in R so is hard for me to understand the steps from the user guide, probably I´m missing something.

Any help with R code for this would be EXTREMELY helpful.

Thank you very much in advance

Gema

[ge2sasag]

please... anyone?

[kalyankpy]

Welcome to the R world! I do not have a straight forward answer for you. But I share my personal experience of working with a specific package with out knowledge of R

My experience with R started in a similar way. I didn't know R but I wanted to analyse. After fiddling for few days, I realized that I took hours and hours understanding a specific function from the package. I started questioning if this is the best way moving forward. Because, I don't know R and don't know how to get help, it was bleak ahead.

Then I shifted to learning basics of R (Datacamp, Coursera, QuickR pages). It is time taking but it helped me understand the very basics of R. It helped me understand how to get help, how to read help pages and interpret it. Surprsingly, when I started my anlaysis again, it was a smooth ride. I could get help easiy since I know R terms (eg: library in R is different from any other library)

R is wonderful. Learning the basics would speed up things rather than starting with a particular package. Also, start using R with Rstudio. It is a different world of R: much more easier to work with R than from terminal!

[ge2sasag]

Thank you very much for your answer and advice!
You described exactly how I feel... Actually I enroled in R course from Coursera 2 weeks ago but I haven´t started with it yet.... But after your post, I will do asap for sure, I need to learn but I always delay it until I really need and I got stuck

[raphael123]

First you should get the strand information for each position !
Could you show the first lines of your bed file ?

[ge2sasag]

Hi raphael123, thanks for answering,

The file I was given is the output from Bowtie2 and is looking like this:

chr1 564495 564496 1 1
chr1 564501 564502 1 1
chr1 565013 565014 1 1
chr1 565040 565041 1 1
chr1 565262 565263 0 8
chr1 565397 565398 1 5
chr1 565469 565470 2 2

4th column is methylation and 5th column is coverage, I have one file like this for each sample, i.e, untreated and treated.

[raphael123]

Okay, I am surprise bowtie2 can call methylation ... I would advise you to use BisSNP to call your methylation from aligned file, then you will still have the strand information and can format it to bssmoth format :

assembly position strand class mc h
1 10497 + CG 31 31
1 10525 + CG 6.9998 31
1 10542 + CG 30 30
1 714168 - CG 0 67
1 714170 - CG 0 67
1 714178 - CG 0 68
1 714182 - CG 0 68
1 714199 - CG 0 68
1 714203 - CG 0 5

[ge2sasag]

With BSmooth you can align with Bowtie2 (not included in BSmooth, you have to run separately) or with an align method that is included and is called Merman (ref. Hansen et al. Genome Biology 2012). My sequencing data is WGBS. Is it possible that the strand info was already taken into account?

[ge2sasag]

Bismark software is also using Bowtie2

[raphael123]

Ok I don t know about BSmooth aligner ...

What we do is :
trimming + quality filter
bismark using bowtie2
BisSNP to improve alignment and call methylation
meth diff using bssmooth( we are testing it now)

[ge2sasag]

OK I see, thanks for sharing your pipeline

I´m trying to get the fastq files from our collaborators, because I don´t really know how to get the bsseq object that I need for calling methylation with BSmooth with the files that they sent me...

[dpryan]

Why don't you just use the read.bsmooth() function from bsseq? If nothing else, BSseq objects aren't particularly difficult to create yourself (just post here and I can tell you how).

[ge2sasag]

Hi, thanks for answering I would be very grateful if you could help me to create the BSseq object to run the analysis.

The files that I have are the output from Bowtie2 and look like this (one file per sample):

chr1 564495 564496 1 1
chr1 564501 564502 1 1
chr1 565013 565014 1 1
chr1 565040 565041 1 1
chr1 565262 565263 0 8
chr1 565397 565398 1 5
chr1 565469 565470 2 2

$4 is methylation and $5 is coverage (according to what I´ve been told)

Is it better to create the BSseq object or to run the read.bsmooth function? I don´t know which is the difference sorry

[dpryan]

Assuming that that's the output from BSmooth (I don't use it, it's too slow and inaccurate), the read.bsmooth() function will return a BSseq object.

[ge2sasag]

This is the output from Bowtie2 alignment, it is not included in BSmooth, the one included in BSmooth is Merman