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 To: Bioconductor Subject: [BioC] lumi, Illumina Methylation 450k, and robust methylation calls Message-ID: 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]]