[BioC] DESeq Error: pobs == ps[kAs[i] + 1] is not TRUE

[Robert Castelo]

dear list, and particularly DESeq developers,

using the package DESeq i get an error in nbinomTest() after calling
estimateDispersions() with the arguments method="per-condition" and
sharingMode="gene-est-only".

this does not happen if the call to estimateDispersions() is performed
with method="per-condition" and sharingMode="fit-only".

i'm pasting below the code, adapted from the vignette of DESeq that
reproduces this error, jointly with the output of the error and of
sessionInfo().


thanks!!!!
robert.

=======================================================================

library(DESeq)
library(pasilla)

data(pasillaGenes)
pairedSamples <- pData(pasillaGenes)$type == "paired-end"
countsTable <- counts(pasillaGenes)[ , pairedSamples ]
conds <- pData(pasillaGenes)$condition[ pairedSamples ]
cds <- newCountDataSet( countsTable, conds )

cds <- estimateSizeFactors( cds )
cds <- estimateDispersions( cds, method="per-condition",
sharingMode="fit-only" )
res <- nbinomTest( cds, "untreated", "treated" )

## so far, everything OK, problem arises now

cds <- estimateDispersions( cds, method="per-condition",
sharingMode="gene-est-only" )
Warning message:
In .local(object, ...) :
  in estimateDispersions: sharingMode=='gene-est-only' will cause inflated numbers of false positives unless you have many replicates.
res <- nbinomTest( cds, "untreated", "treated" )
Error: pobs == ps[kAs[i] + 1] is not TRUE
> traceback()
7: stop(paste(ch, " is not ", if (length(r) > 1L) "all ", "TRUE", 
       sep = ""), call. = FALSE)
6: stopifnot(pobs == ps[kAs[i] + 1])
5: FUN(1:14470[[1L]], ...)
4: lapply(X = X, FUN = FUN, ...)
3: sapply(1:length(kAs), function(i) {
       if (kAs[i] == 0 & kBs[i] == 0) 
           return(NA)
       ks <- 0:(kAs[i] + kBs[i])
       ps <- dnbinom(ks, mu = mus[i] * sum(sizeFactorsA), size = 1/sumDispsA[i]) * 
           dnbinom(kAs[i] + kBs[i] - ks, mu = mus[i] * sum(sizeFactorsB), 
               size = 1/sumDispsB[i])
       pobs <- dnbinom(kAs[i], mu = mus[i] * sum(sizeFactorsA), 
           size = 1/sumDispsA[i]) * dnbinom(kBs[i], mu = mus[i] * 
           sum(sizeFactorsB), size = 1/sumDispsB[i])
       stopifnot(pobs == ps[kAs[i] + 1])
       if (kAs[i] * sum(sizeFactorsB) < kBs[i] * sum(sizeFactorsA)) 
           numer <- ps[1:(kAs[i] + 1)]
       else numer <- ps[(kAs[i] + 1):length(ps)]
       min(1, 2 * sum(numer)/sum(ps))
   })
2: nbinomTestForMatrices(counts(cds)[, colA], counts(cds)[, colB], 
       sizeFactors(cds)[colA], sizeFactors(cds)[colB], rawScvA, 
       rawScvB)
1: nbinomTest(cds, "untreated", "treated")

sessionInfo()
R version 2.14.0 (2011-10-31)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=C                 LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] pasilla_0.2.10      BiocInstaller_1.2.1 DESeq_1.6.0        
[4] locfit_1.5-6        lattice_0.20-0      akima_0.5-4        
[7] Biobase_2.14.0     

loaded via a namespace (and not attached):
 [1] annotate_1.32.0      AnnotationDbi_1.16.2 DBI_0.2-5           
 [4] genefilter_1.36.0    grid_2.14.0          IRanges_1.12.1      
 [7] RSQLite_0.10.0       splines_2.14.0       survival_2.36-10    
[10] tools_2.14.0         xtable_1.6-0