Hi, 

according to my script here below, why DESEq2 calculations and my own are not identical? What am I missing?

Any halp will be most welcome, thanks.

David

dds<- makeExampleDESeqDataSet(n=100,m=18) 
dds <- DESeq(dds, betaPrior = FALSE)
resultsNames(dds)
res <- results(dds)

res[1,2]
# [1] -0.3811565
my.log2FoldChange <- log(rowMeans(counts(dds[,10:18],normalized=TRUE))[1]/rowMeans(counts(dds[,1:9],normalized=TRUE))[1],2)
my.log2FoldChange
# -0.3792306

res[1,3]
# [1] 0.2197007
my.lfcSE <- sd(log(counts(dds[1,10:18],normalized=TRUE)/counts(dds[1,1:9],normalized=TRUE),2))/sqrt(8)
# [1] 0.2036169

For the log2FoldChange, the reason for the small difference is that we have in our model some bounding at the very low end of the estimated mean value which is conceptually equivalent to using a fractional pseudocount. This helps with stability of inference on fold changes and their standard error. You would only notice a difference when the counts are very small. The random data generator makeExampleDESeqDataSet() by default is making fairly small counts, mostly because I use this function to generate data with a number of rows with all 0's for software testing purposes.

For the standard error, that formula is not correct for the calculation of standard error of a coefficient for any model I know of. We give the formula we use for the standard error of a GLM coefficient in the DESeq2 paper.