[BioC] bsseq and BSmooth with RRBS

[Alex Gutteridge]

On 21.02.2013 17:52, Kasper Daniel Hansen wrote:
> On Thu, Feb 21, 2013 at 11:31 AM, Alex Gutteridge
> <alexg at ruggedtextile.com> wrote:
>> I'm looking over the documentation for bsseq:
>>
>> http://www.bioconductor.org/packages/2.11/bioc/html/bsseq.html
>>
>> And the associated Genome Biology paper:
>>
>> http://genomebiology.com/content/pdf/gb-2012-13-10-r83.pdf
>>
>> From what I can see the method should work on reduced representation
>> bisulfite sequencing (RRBSeq) datasets as well as WGBSeq, though the
>> documentation only mentions WGBSeq. The only issue I can see being 
>> the
>> smoothing procedure which I thought might struggle in regions with 
>> scattered
>> CpGs in the RRBSeq.
>>
>> I will have data in hand shortly to test this myself, but was just 
>> curious
>> if anyone on the list had run bsseq (successfully or otherwise) on 
>> RRBS data
>> before?
>
> I don't have any hands on experience with RRBS, but I agree with your
> assessment that in principle it should work, depending on how large
> regions are actually captured in your experiment.  Mathematically, it
> should run out of the box.  Pay close attention to the maxGap
> parameter which sets a limit for how far away two neighboring CpGs 
> can
> be before not smoothing.
>
> I have had inquires like yours about using it for RRBS data, but I
> never heard back about their experience (and as I recall it, it was
> much earlier in the development phase, so the software ought to work
> much better now).
>
> I would be happy to hear your experience.
>
> Note that currently, (some of) the statistics are really meant for 
> the
> case where you have biological replicates, but that question is 
> really
> orthogonal to the choice of assay.
>
> Best,
> Kasper

Thanks Kasper,

I will give it a go and report back. I was also looking at modelling 
this on an individual locus basis with a binomial glm. I.e. With made up 
data:

> M
[1] 10 10 10 30 40 10
> U
[1] 10 20 15  3 10  0
> y = cbind(M,U)
> group = factor(c(0,0,0,1,1,1))
> model = glm(y~group,binomial)
> summary(model)

Call:
glm(formula = y ~ group, family = binomial)

Deviance Residuals:
       1        2        3        4        5        6
  0.9036  -0.7537   0.0000   0.8569  -1.1656   1.7354

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4055     0.2357  -1.720   0.0854 .
group1        2.2225     0.3808   5.837 5.31e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

     Null deviance: 46.8213  on 5  degrees of freedom
Residual deviance:  6.4888  on 4  degrees of freedom
AIC: 28.198

Number of Fisher Scoring iterations: 4

Where M is the number of methylated reads at a single locus across 6 
samples and U is the number of unmethylated reads at the same locus in 
the same
samples. Group is a factor giving the design of the experiment (here 
two groups both with n=3).

Is the conclusion from your work that pre-smoothing the data gives 
better performance than this naive way for detecting DMRs? If we were 
interested in testing specific CpG loci (as opposed to broader regions) 
would this approach make sense?

-- 
Alex Gutteridge