Hi Tim

The algorithm for methylation status estimation has not been well tested
yet. That's why I did not put the details on the vignette. Now I am working
on other parts of 450K methylation analysis, like DMR detection and
annotation. I will work on methylation status estimation later when I have
time. Sorry about it.


Pan


Date: Wed, 8 Jun 2011 12:38:58 +0100
From: Tim Rayner <tfrayner@gmail.com>
To: Bioconductor <bioconductor@stat.math.ethz.ch>
Subject: [BioC] lumi, Illumina Methylation 450k, and robust
       methylation calls
Message-ID: <BANLkTi=7TU-5WTrOxhk0qRaUcneVhnQqRw@mail.gmail.com>
Content-Type: text/plain; charset="ISO-8859-1"

Hi,

I have recently started to use the lumi package to analyse some
Illumina Human Methylation 450k data, and I have run into some
problems which seem to revolve around division by zero in the
gammaFitEM() function. I have adjusted the colour balance and quantile
normalised as suggested in the vignette, but when I call gammaFitEM()
the function complains (see the end of this email for a session dump).

I've traced the likely cause of the error to zero values returned by
these calls in gammaFitEM():

       f1 <- dgamma(x - s[1], shape = k[1], scale = theta[1])
       f2 <- dgamma(s[2] - x, shape = k[2], scale = theta[2])

The problem is that the z1 variable subsequently contains divisions by
these zero values:

       z1 <- p[1] * f1/(p[1] * f1 + p[2] * f2)

When z1 is later used in calls to sum() in many places, this obviously
returns NaN which causes the function to raise an exception. I think
I've got around this by editing the function and putting na.rm=TRUE in
each of the relevant calls to sum(), and the generated plots look
quite believable, but I can't be sure if that's actually a valid
approach. Is there a better way to address this problem?

Many thanks,

Tim


--
Tim Rayner
Bioinformatician, Smith Lab,
CIMR, University of Cambridge, U.K.



## session dump follows:
> library(lumi)
Loading required package: Biobase

Welcome to Bioconductor

 Vignettes contain introductory material. To view, type
 'browseVignettes()'. To cite Bioconductor, see
 'citation("Biobase")' and for packages 'citation("pkgname")'.

Loading required package: nleqslv
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009

Attaching package: 'lumi'

> # data.c.quantile is the normalised MethyLumiM object.
> fit<-gammaFitEM(exprs(data.c.
quantile)[,1], plotMode=TRUE, verbose=TRUE)

Initial estimation:
k: 28 8
s: -10.17401 5.376681
theta: 0.2424133 0.4134306
p: 0.4264381 0.5735619
logLikelihood: -1135672
Error in if (abs(p.new[1] - p[1]) < tol) break :
 missing value where TRUE/FALSE needed

> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lumi_2.4.0     nleqslv_1.8.5  Biobase_2.12.1

loaded via a namespace (and not attached):
 [1] affy_1.30.0           affyio_1.20.0         annotate_1.30.0
 [4] AnnotationDbi_1.14.1  DBI_0.2-5             grid_2.13.0
 [7] hdrcde_2.15           KernSmooth_2.23-5     lattice_0.19-26
[10] MASS_7.3-13           Matrix_0.999375-50    methylumi_1.8.0
[13] mgcv_1.7-6            nlme_3.1-101          preprocessCore_1.14.0
[16] RSQLite_0.9-4         tools_2.13.0          xtable_1.5-6

	[[alternative HTML version deleted]]

