basemean and log2FoldChange correlation in DESeq2

[gonzalezem]

Hi!

I am using DESeq2 within Phyloseq (as you described here) and I came accross some features that I can't figure out on my own. I used 2 functions to get my DE results:

physeq_ds = phyloseq_to_deseq2(physeq, ~ condition)
diagdds = DESeq(physeq_ds, test="Wald", fitType="parametric")`

From diagdds, I got the results in a table-like output with:

res <- results(diagdds)

The output table has the following columns:
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj

Also from diagdds, I created my own basemean column with:

dds <- estimateSizeFactors(diagdds)
basemean_norm=rowMeans(counts(dds, normalized=TRUE))

And finally, from Michael Love's input here , I calculated the basemean for each factor:

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

The idea behind this was to generate a bar chart to visualize each condition's abundance and compare it with the DE features. I saw some discrepancy between the sign of the log2FoldChange and the factor's basemeans. That's my main problem. I also saw difference in the basemeans calculated separately. I am comparing a factorA (74 samples) and a factorB (10 samples) within a condition. Here is the table that puzzles me ( I rounded the mean values):

taxon	basemeanA	basemeanB	basemean_norm	baseMean	log2FoldChange
taxon1	2480	129	2200	310	5.72577586
taxon2	3204	6753	3627	1218	5.73803002

For many rows, when basemeanB is higher, log2FoldChange is negative. But it's not always the case and the most striking example is on the table. These 2 have close log2FoldChange and really close padj although the normalized mean counts seem to suggest they should be on the opposite side. What am I missing? Thanks!

[gonzalezem]

here and here are the answers to this. The problem seems most likely to come from the data transformation.