Entering edit mode
Hi,
I'm using DESeq2 to analyse some RNA-seq data. I have different treatments (XXX, YYY and ZZZ) at different time points (2h, 4h, 8h), with DMSO used as a control for the treatment. 'se' is my SummarizedExperiment object containing the raw counts.
## create DESeqDataSet
dds <- DESeqDataSet(se, design = ~ treatment + time + treatment:time)
# make sure DMSO is the base level of the treatment factor
dds$treatment <- relevel(dds$treatment, "DMSO")
## differential expression analysis
dds<- DESeq(dds)
## get results
# let's check the effect of treatment XXX
res <- results(dds, contrast=c("treatment", "XXX", "DMSO"))
## MA plot
plotMA(res)
Now, this is what my MA plot looks like:
Does this look right? It almost looks like the fold changes have not been shrunken! Is that the case?
Am I doing something wrong?
Thank you.
Thanks Michael. Why is the prior only applied to the interaction terms?
I very often do designs with interaction terms because I'm interested in, say, condition A in group 1, but I might also be interested in condition A over all groups. Since the prior is not applied to the second example, how reliable are its fold changes?
Thank you.
It takes a bit of geometric intuition, but for certain settings it was possible to observe that both shrinking the main effect terms and the interaction term resulted in the interaction term growing (despite the prior) to account for the difference between the observed counts and the predicted value using only the main effect terms. This could be demonstrated when the interaction term was truly 0, that the distribution of Wald or likelihood ratio statistic for the estimated interaction effects was too wide or too long tailed. This does not occur in the current implementation.
The unshrunken fold changes are perfectly reliable, if you filter on adjusted p-value.