[BioC] edgeR, glmLRT error

[Gordon K Smyth]

Hi Cei,

This is now fixed in edgeR 2.99.3 on Bioc-devel.  Your example now runs 
with the results given below.  Note that glmLRT now operates on the fit 
object directly.  The data object is no longer required as an argument.

Best wishes
Gordon

   lrt <- glmLRT(fit, contrast=c(0,0,0,-1,0.5,0.5))
   topTags(lrt)

Output:
Coefficient:  -1*T2_S1 0.5*T2_S2 0.5*T2_S3
                           logFC logCPM     LR   PValue      FDR
CA00-XX-RX1-099-H06-EB.F 16.903   19.2 37.366 9.79e-10 1.08e-08
Contig1405                2.446   18.0 10.440 1.23e-03 6.78e-03
Cc_Contig7187             2.168   17.1  7.449 6.35e-03 2.33e-02
Contig2423               -1.240   16.3  4.252 3.92e-02 1.08e-01
CA00-XX-LP1-022-G05-EB.F -1.064   12.5  2.722 9.90e-02 2.18e-01
Contig10618               0.699   13.1  0.913 3.39e-01 6.22e-01
Cc_Contig2487            -0.506   15.8  0.658 4.17e-01 6.56e-01
Contig6686               -0.387   12.6  0.337 5.62e-01 7.72e-01
Contig15077              -0.279   15.8  0.198 6.56e-01 7.82e-01
Contig7207               -0.247   12.2  0.137 7.11e-01 7.82e-01


On Tue, 17 Jul 2012, Cei Abreu-Goodger wrote:

> Hi guys,
>
> Many thanks for the update, I'll check the devel in a few weeks.
>
> Best,
>
> Cei
>
> On 7/16/12 11:35 PM, Gordon K Smyth wrote:
>> Dear Cei,
>> 
>> Thanks for the nice reproducible example.
>> 
>> We have been working on a programming solution to eliminate this
>> occasionally occuring error in edgeR.  We have a nearly complete
>> solution in our local (non-public) version of edgeR.  Running your data
>> example with our current local code, glmLRT no longer gives an error,
>> but the 9th row of your data generates a NA value for the likelihood
>> ratio test statistic.
>> 
>> That's already an improvement, but I want to eliminate the NA fits as
>> well if possible.
>> 
>> We will push this code out to Bioc-devel soon, but probably not before
>> the Bioconductor conference next week.
>> 
>> Best wishes
>> Gordon
>> 
>>> Date: Sun, 15 Jul 2012 18:36:16 -0500
>>> From: Cei Abreu-Goodger <cei at langebio.cinvestav.mx>
>>> To: bioconductor at r-project.org
>>> Subject: [BioC] edgeR, glmLRT error
>>> 
>>> Hi Gordon,
>>> 
>>> I've been getting the following error upon running glmLRT:
>>> 
>>> Error in beta[k, ] <- betaj[decr, ]
>>>   NAs are not allowed in subscripted assignments
>>> 
>>> which has previously been reported, for instance:
>>> 
>>> https://stat.ethz.ch/pipermail/bioconductor/2011-February/038079.html
>>> 
>>> I subsetted my count data (originally ~14k rows) to narrow down the
>>> culprit row, and I would like to send you this example in case it can be
>>> used to improve the glm fitting code. If it turns out to be impossible
>>> to completely avoid these kinds of errors, perhaps it would be nice to
>>> return the genes that are generating the NAs so that we can more easily
>>> remove them and proceed?
>>> 
>>> I can also avoid this particular case by using only the CommonDisp or
>>> heavily increasing the prior.n, but that kind of defeats the purpose of
>>> the TagwiseDisp.
>>> 
>>> Example code and sessionInfo:
>>> 
>>> ################
>>> library(edgeR)
>>> 
>>> counts <-
>>> read.table("http://datos.langebio.cinvestav.mx/~cei/counts.tab")
>>> grp <- sub("_R.+", "", colnames(counts))
>>> dge <- DGEList(counts=counts, group=grp)
>>> dge <- calcNormFactors(dge)
>>> 
>>> design <- model.matrix(~0 + dge$samples$group)
>>> colnames(design) <- levels(dge$samples$group)
>>> 
>>> dge <- estimateGLMCommonDisp(dge, design)
>>> dge <- estimateGLMTagwiseDisp(dge, design)
>>> 
>>> fit <- glmFit(dge, design)
>>> lrt <- glmLRT(dge, fit, contrast=c(0,0,0,-1,0.5,0.5))
>>> 
>>> Error in beta[k, ] <- betaj[decr, ] :
>>>   NAs are not allowed in subscripted assignments
>>> 
>>> ################
>>> R version 2.15.0 (2012-03-30)
>>> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>>> 
>>> locale:
>>> [1] C/UTF-8/C/C/C/C
>>> 
>>> attached base packages:
>>> [1] grDevices datasets  stats     graphics  utils     methods   base
>>> 
>>> other attached packages:
>>> [1] edgeR_2.6.9         limma_3.12.1        Biobase_2.16.0
>>> BiocGenerics_0.2.0  BiocInstaller_1.4.7
>>> 
>>> ################
>>> 
>>> 
>>> Many thanks,
>>> 
>>> Cei
>>> 
>>> --
>>> Dr. Cei Abreu-Goodger
>>> Langebio CINVESTAV
>>> Tel: (52) 462 166 3006
>>> cei at langebio.cinvestav.mx

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}