DESeq2 shrinking log fold-change

Question

[JACKY]

I'm conducting an RNA-seq gene differential analysis using DESeq2. I'm trying to shrink LFC using apeglm, ashr or even normal, but something strange is happening when plotting MA plot, the log fold change is going bigger rather than shrinking. Even when using ashr or normal, it does a similar thing, which is very weird. I am using of course raw data without cleaning it (just removed duplicates), so the low reads and zeros are still there. Then I tried to remove them before hand, cleaned the all zero rows, but I got similar result.

And the exact same thing happens when using ashr. Why is this happening ? Should I just use the orinigal DESeq results and not do shrinkage at all ? Maybe just remove the lowly expressed genes by hand and not do the shrinking ?

Answer 1

[Michael Love]

Did you demonstrate a LFC that is larger after posterior inference compared to MLE?

[JACKY]

hmm I don't realy understand your question but comparing the results before and after shrinkage.. they are really similar. Looking at the first 20 genes in the shrinking results, yes there is some shrinkage in the LFC and some genes just stayed the same. I tried addMLE = TRUE but it gives me an error. I made a plot of LFC and shrunk LFC, does this help ?

[Michael Love]

The "without shrinkage" seems to be with shrinkage. Not sure how you are getting that out of DESeq().

Try the following (assuming this is two group comparison):

# very minimal filtering
keep <- rowSums(counts(dds) >= 5) >= 3 
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds)
lfc <- lfcShrink(dds, coef=2, type="apeglm")
plot(res$log2FoldChange, lfc$log2FoldChange); abline(0,1)

[JACKY]

Yes it's a two group comparison. This is the output of your code, I also did the MA plot.. I'm a bit lost of what these results mean though. I guess it looks better than before, but I don't think the fold change is shrunk as it should be, right ?

[Michael Love]

Something seems strange about your DESeq(). What versions are you using sessionInfo().

Above it appeared to be shrinking though the default value is betaPrior=FALSE. It is as if your defaults are not set the same.