DESeq versions: mimic 1.4.1 with 1.6.1 settings?

[daler]

I am unable to replicate DESeq v1.4.1 results using v1.6.1, even when using the settings that -- as far as I can tell from the docs -- should replicate the old behavior. Here's a self-contained working example . . . but it needs parallel installations of R 2.13.1 and R 2.14.0 in order to work.

First, I created data using only v1.6.1 and saved it to file:

library(DESeq)
cds <- makeExampleCountDataSet()
write.table(counts(cds), file='example.counts')

In R 2.13.1 I ran DESeq v1.4.1:

library(DESeq)
x <- read.table('example.counts')
conds <- c('A', 'A','B','B','B')
cds <- newCountDataSet(x, conds)
cds <- estimateSizeFactors(cds)

cds <- estimateVarianceFunctions(cds, method='normal')

res <- nbinomTest(cds, 'A', 'B')
write.table(res, file='old.results', sep='\t', row.names=F)

Then, over in R 2.14.0, I ran DESeq v1.6.1. Note that everything except the "estimateDispersions" line is the same:

library(DESeq)
x <- read.table('example.counts')
conds <- c('A', 'A','B','B','B')
cds <- newCountDataSet(x, conds)
cds <- estimateSizeFactors(cds)
cds <- estimateSizeFactors(cds)

cds <- estimateDispersions(cds, sharingMode='fit-only', 
                           fitType='local', method='per-condition')

res <- nbinomTest(cds, 'A', 'B')
write.table(res, file='new.results', sep='\t', row.names=F)

When I compare new.results with old.results, basemeanA, basemeanB, and the fold change columns are identical.

However, the pval and padj columns are different; plotting them results in two straight-ish lines on either side of the 1:1 (see attached PNG):

new = read.table('new.results', header=T)
old = read.table('old.results', header=T)
plot(new$pval, old$pval)
abline(0, 1, col='red')

What could be causing this discrepancy? Are there other parameters to estimateDispersions that I'm missing? Has something changed in nbinomTest between versions?

-ryan

[Simon Anders]

Yes, we did change 'nbinomTest'. It used to employ an approximation that usually would only be a few percent off (which, for p values, does not matter; one is only interested ion the magnitude, after all), but gave in a few rare cases drastically wrong results. Realizing that this approximation was not really necessary anyway, we removed it. See the end of the new vignette, by the way, for a summary of this and related changes.

[daler]

Ah, your explanation here plus the note in the vignette (which I somehow missed) clears it up -- thanks.