Entering edit mode
Pan Du
▴
440
@pan-du-4636
Last seen 10.2 years ago
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]]