For DESeq2, I am trying to understand the interpretation of with/without an interaction term in the design formula. Assume we have two conditions (Treatment vs Control) and two genotypes (Mutant vs WT). Here is a simulation study and my interpretation.
> library("DESeq2") # make example dataset and fix the reference level > dds <- makeExampleDESeqDataSet(n=10000,m=12) > dds$condition <- as.factor(c(rep("Control",6), rep("Treatment",6))) > dds$genotype <- as.factor(rep(c(rep("WT", 3), rep("Mutant",3)),2)) > dds$genotype <- relevel(dds$genotype, "WT") > as.data.frame(colData(dds)) > head(assay(dds))
Case1: With interaction
I could include the interaction term in the design, as I showed in the following code. Then I would retrieve the result of conditionTreatmentvs_Control. From what I understand, this is comparing Treatment and Control under WT, or comparing sample 7-9 with sample 1-3.
> design(dds) <- ~ condition + genotype + condition:genotype > dds <- DESeq(dds) > results(dds, name = "condition_Treatment_vs_Control")[1:5,]
Case2: Without Interaction
I could write the design formula without the interaction term and still, I want to retrieve the result of conditionTreatmentvs_Control. This is technically the same as model the batch effect. DESeq2 will regress out the effect of genotype. In other words, DESeq2 will take all samples into consideration.
> design(dds) <- ~ condition + genotype > dds <- DESeq(dds) > results(dds, name = "condition_Treatment_vs_Control")[1:5,]
But if you look at the results, they are very different. I want to validate my point by manually calculate the log2FoldChange column.
Here I try to calculate the log2foldchange for gene 1, with the interaction term included in the design formula
> log2(mean(counts(dds, normalized = T)[1,7:9]) / mean(counts(dds, normalized = T)[1,1:3]))
This is almost the same as the result of DESeq2 (with interaction).
So here is my question: 1. Is my interpretation of the design formula correct? 2. How to manually calculate the log2foldchange for case2 (without interaction)?