Entering edit mode
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
