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:
In R 2.13.1 I ran DESeq v1.4.1:
Then, over in R 2.14.0, I ran DESeq v1.6.1. Note that everything except the "estimateDispersions" line is the same:
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):
What could be causing this discrepancy? Are there other parameters to estimateDispersions that I'm missing? Has something changed in nbinomTest between versions?
-ryan
First, I created data using only v1.6.1 and saved it to file:
Code:
library(DESeq) cds <- makeExampleCountDataSet() write.table(counts(cds), file='example.counts')
In R 2.13.1 I ran DESeq v1.4.1:
Code:
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)
Code:
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)
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):
Code:
new = read.table('new.results', header=T) old = read.table('old.results', header=T) plot(new$pval, old$pval) abline(0, 1, col='red')
-ryan
Comment