DESeq Error: pobs == ps[kAs[i] + 1] is not TRUE
1
0
Entering edit mode
Robert Castelo ★ 3.1k
@rcastelo
Last seen 1 day ago
Barcelona/Universitat Pompeu Fabra
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
DESeq DESeq • 984 views
ADD COMMENT
0
Entering edit mode
Simon Anders ★ 3.7k
@simon-anders-3855
Last seen 3.7 years ago
Zentrum für Molekularbiologie, Universi…
Hi Robert On 2011-11-17 15:20, Robert Castelo wrote: > using the package DESeq i get an error in nbinomTest() after calling > estimateDispersions() with the arguments method="per-condition" and > sharingMode="gene-est-only". Thanks for the bug report. I just corrected the issue; it is now fixed in versions 1.6.1 (release) and 1.7.2 (devel). Please remember to use "gene-est-only" only if you have really many replicates (let's say, ten or more). Otherwise (as the warning issued by the function explains) you will lose type-I error control. Simon
ADD COMMENT

Login before adding your answer.

Traffic: 821 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6