Tim Rayner <tfrayner at="" ...=""> writes: > > 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 > Hi Tim, I have your same problem even if with the Illumina Human Methylation 27k. You wrote about editing any sum() function with na.rm=TRUE. It really worked for you ? I also made the same changes but then R says that did not find the ".initialGammaEstimation function" Such function is internal to the gammaFitEM package ... How you solved the problem ? In the end the results are valid or really we can't use this function at all ? Thanks, Pietro

Hi Pietro, You can replace the call to .initialGammaEstimation() with lumi:::.initialGammaEstimation(). It's a fairly shameless hack which may well fail to work with subsequent versions of lumi, but with any luck that point is moot. Having looked at our results I'd say I'm pretty confident that the function is doing at least roughly what one would expect. I do have some concerns about applying this function to the 450k data, since the presence of both Infinium I and II assays on the chip generates a distribution with more than two density peaks. It might be worth separating the probes into I and II and applying the function to each set separately. In our case we're only using these calls for a quick look at the gross methylation changes between samples; to calculate DMRs we're using limma as suggested in the lumi vignette. Cheers, Tim On 4 August 2011 01:24, zoppoli pietro <zoppoli.pietro at="" gmail.com=""> wrote: > Tim Rayner <tfrayner at="" ...=""> writes: > >> >> 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 >> > Hi Tim, > > I have your same problem even if with the Illumina Human Methylation 27k. > You wrote about editing ?any sum() function with na.rm=TRUE. > It really worked for you ? > I also made the same changes but then R says that did not find the > ".initialGammaEstimation function" > Such function is internal to the gammaFitEM package ... > How you solved the problem ? > > In the end the results are valid or really we can't use this function at all ? > > Thanks, > Pietro > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >