Hello,

I'd like to know how the log2 fold change is calculated between target and comparison population in DEXSeq. Going over the estimateExonFoldChanges function in an older version (0.12.1) of the package, I realize the interaction coefficient is taken from the model: count ~ condition * exon and fold change is calculated by applying a variance stabilizing transformation and then transformed to a log2 scale:

alleffects <- do.call(rbind, alleffects)
alleffects <- vst(exp(alleffects), object)
alleffects <- log2(alleffects/alleffects[, denoCol]).  ###foldChange <- effects[,"target"] - effects[,"comparison"]

alleffects data frame looks like the following:

Gene target comparison
Ptma 5.36425504487572   5.10532234811512
Ptma 4.43604234783272   4.7521435893567
Ptma 4.30355887270297   4.72294913039353
Tmpo 2.12202872975088   1.08346386248873
Msn 1.86941999824138    2.34083780006062

However, in a newer version that I am currently using (1.28.3), looks like the target and comparison values being used for fold change calculation are not vst transformed. Also, being divided by log(2).

alleffects <- rbind(alleffectsBM, alleffectsSM)
alleffectsVst <- vst(exp(alleffects), object)
alleffects <- alleffects/log(2)
alleffects <- alleffects - alleffects[, denoCol] ###foldChange

The calculation in the older version makes more sense for a log2 fold change than in version 1.28.3.

An explanation for this discrepancy in fold change calculation would be appreciated.

Thanks in advance!

Hi @tmvarsha,

Thanks for your detailed report. The first version actually had a bug: vst data is already log-like, so the code was wrongly calculating a log2 fold change from a log-like data. In the second version, you are right that it is not variance-stabilized transformed. The appropiate approach would be to use the shrinkage approaches of DESeq2: I've been thinking on how to implement this but I need to find the time to do this.

Alejandro