[apredeus]

Oh so I guess you must have R3.1 and later to get the newest DESeq2? Since I've installed the package very recently but I did have R 3.0.1..

[apredeus]

So when I upgraded DESeq2 1.2.10 @ R 3.0.1 to v 1.4.5 (R 3.1.1) the number of genes considered "differentially expressed" went from ~ 2500 to ~ 1300.

every other setting is very much the same. How is this even possible?

[Michael Love]

hi apredeus, I'm note exactly sure why it would have changed so much between 1.2 and 1.4, as there were only slight changes to the gene-wise dispersion estimation routine, but I can say that the code base has stabilized since last fall. Bioc users are always 6 months behind the development cycle, which means that bugs/changes that users request in October have to wait until March to be released.

The other thing is that p-values are tail probabilities, which means they are very sensitive to slight changes in model parameters.

For example, say we have some Z scores:

> set.seed(1)
> z <- rnorm(10000,0,7)

Then we make and adjust p-value according to two models of the data:

> table(p.adjust(2*pnorm(abs(z),0,5,lower=FALSE),method="BH") < .1)

FALSE TRUE
9507 493

> table(p.adjust(2*pnorm(abs(z),0,5.25,lower=FALSE),method="BH") < .1)

FALSE TRUE
9741 259

So a 5% change in a model parameter translates to a nearly halving of adjusted p-values less than 0.1.

[apredeus]

Michael, Simon:

would it be possible to output two mean expression columns (conditions 1 and 2) instead of only one in the results table? It would be a lot more informative I think. You can simply glance at it and tell if there are some dramatic changes.

[kmcarr]

Michael, Simon:

would it be possible to output two mean expression columns (conditions 1 and 2) instead of only one in the results table? It would be a lot more informative I think. You can simply glance at it and tell if there are some dramatic changes.

This post by Michael described how to extract the baseMean per level (condition) from the DESeqDataSet object.

[Michael Love]

I can make it even more concise:

baseMeanPerLvl <- sapply(levels(dds$condition), function(lvl) rowMeans(counts(dds,normalized=TRUE)[,dds$condition == lvl]))

I didn't know until this year that dds$variable is the same as colData(dds)$variable for SummarizedExperiments.