[BioC] lumi, Illumina Methylation 450k, and robust methylation calls

[Tim Rayner]

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