Hello,

I have a small dataset with 3 replicates per condition. One condition is a gene over-expression, the other is a control. I'm used to run DESeq2 with betaPrior = True to have a comparability with the 'old' DESeq2 behavior. I can't apply this old workflow because the expression differences between my conditions it quite small and one gene, namely the over-expressed  gene, would stand out and with betaPrior = T it's shrunk to a very small log2 FoldChange. According to the manual this is to be expected which is why I used betaPrior = F and the several lfcShrink methods (normal, ashr and apeglm).

However I was wondering if I should apply a shrinkage or just use betaPrior = F? Secondly, I find the shrinkage quite drastically and not comparable to the manual (https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#alternative-shrinkage-estimators). I know you can't compare datasets but still it's good to use it as a reference. It's also explained here that imprecise fold changes are shrunk and large fold changes are kept (https://support.bioconductor.org/p/77461/). And that is valid for the data set. The DEGs have a similar log2Foldchange over the used methods.

Here is the data set with 3 replicates and the several methods. From the left to the right it is the following and I used coef in the lfcShrink function:

betaPrior = F and method = normal, betaPrior = F, betaPrior = F and method = ashr, betaPrior = F and method = apeglm

https://ibb.co/cwERg8

If I use two replicates, it's even more visible, especially for apeglm and I wonder what's happening here. Especially most of the not DE are reduced to log2Foldchange of < 0.001 and even some DEG have a log2Foldchange < 0.001

https://ibb.co/fjz6ET

Thanks

Mathias

dear Mathias,

Thanks for posting these images. What you observe is consistent with what we see in testing on the benchmarking data and on simulation data. If you just compare method="normal" to method="apeglm" or "ashr", the differences you are likely to see is that normal will shrink large effects even if they have high precision (so shrinking too much) and allow small effects to float around 0, while apeglm/ashr will not shrink the precise, large effects much at all and the small effects which are indistinguishable from 0 will be shrunk to 0. By indistinguishable I mean that the SE(beta) >> |beta|. This kind of estimator is very good at ranking genes and it also has low error compared to the true LFC (which we know in the benchmarking datasets and the simulation datasets). Also all of these estimators converge to MLE as the sample size gets larger. It's up to you which you want to use for plots. For ranking genes by effect size, I'd recommend apeglm or ashr, as these perform the best in our hands.

One more note: the DEG according to the Wald test can be shrunk a lot by apeglm or ashr, they are different approaches statistically, and you can have a very small p-value for small effect size which nevertheless is not equal to 0. If you want to have more consistent error-bounded sets (FDR/FSR sets) with the shrunken effect sizes, you can use svalue=TRUE when you call lfcShrink(). Then the error-bounded set is consistent with the shrunken effect size, because the same posterior is being used for both.