Recently, I reanalyzed some gene expression data using DESeq2. To my surprise, the resulting volcano plots looked quite different from the previous analysis. During the troubleshooting phase of the analysis, I realized that the data was previously generated using a Galaxy-based version of DESeq2 v 1.34.0 and I used version 1.40.1 in RStudio (see below).

Close inspection of the results files showed that the calculation of the log2FoldChange was the major reason for the discrepancy. I am including (see below) as an example, information about gene A. I am also including the raw numbers for the gene, the normalized counts coming from the dds object (RStudio), and the VST normalized counts for both Galaxy and RStudio files.

1) Would you please comment which fold change will be closer to reality?

2) Why am I getting such a discrepancy in the log2FoldChange between the analyses?

Thank you for your help!

Platform    baseMean    Log2FoldChange  lfcSE       stat        pvalue      padj
Galaxy      17.102747   -0.44276    0.15465304  -2.86295    0.004197    0.0227
Rstudio     17.1026     -5.92047    1.49447     -3.961587   7.45E-05    0.000820

Samples             CntRNA-1    CntRNA-2    CntRNA-3    KORNA-1 KORNA-2 KORNA-3
Raw counts          1       78      6       1   0   1
Normalizeddds-RStudio       1.120692    92.55239    7.274768    0.68380 0   0.983963
VST-Counts-RStudio      9.446782    9.920987    9.538052    9.43386 9.3877  9.443066
VST-Counts-Galaxy       7.477439    8.417805    7.661466    7.45136 7.3582  7.469936

I guess that Galaxy may be using some kind of LFC shrinkage, but without seeing their code I wouldn't know for sure.

Shrinkage generally gives much more accurate LFC. I would recommend using apeglm or ashr methods in lfcShrink.