Entering edit mode
Tim Rayner
▴
270
@tim-rayner-2913
Last seen 10.2 years ago
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