I am looking more for clarification about how DESeq2 handles multiple conditions/genotypes. We are interested in looking at a WT/Mut with treatment based on the analysis for two comparisons, two genotypes with interaction outline here. Our data is WT:Treat, WT:NoTreat, Mut:Treat, Mut:NoTreat. I want to know what genes are upregulated/important for each of the different versions of our experiment; Mut, WT, Treat, NoTreat, but also what is different based off the interactions (do the genes upregulated in the treatment change in the mutant compared to WT). Can DESeq2 handle this type of question and is the following code how we should approach this? Does DESeq2 account for the condition or genotype when making the reciprocal comparison?

I am not doing a great job explaining this, but I really want to be able to confidently say that when I compare our mutant to WT the treatment alters X number of genes differently than when looking at just WT and all permutations within these (WT vs Mut NO treatment and also Treatment vs no Treatment Mut etc.).

dds <- DESeqDataSetFromMatrix(as.data.frame(rawdata),
                              as.data.frame(design_deseq),
                              ~ genotype + condition + genotype:condition)

res = results(dds, contrast=c("condition","Geno1","Geno2"))

res = results(dds, contrast=c("genotype","Cond1","Cond2"))

Our data is WT:Treat, WT:NoTreat, Mut:Treat, Mut:NoTreat. I want to know what genes are upregulated/important for each of the different versions of our experiment; Mut, WT, Treat, NoTreat

Hi, you may need to explain this, specifically what is meant by "each of the different versions" - you mean you want to find genes uniquely expressed in, e.g., Mut? This cannot really be addressed via differential expression comparisons. You may want to instead go the route of transforming the dataset to Z-scores.

Regarding the interaction, have you checked the examples accessible via ?DESeq2::results?

 ## Example 2: two conditions, two genotypes, with an interaction term

 dds <- makeExampleDESeqDataSet(n=100,m=12)
 dds$genotype <- factor(rep(rep(c("I","II"),each=3),2))

 design(dds) <- ~ genotype + condition + genotype:condition
 dds <- DESeq(dds) 
 resultsNames(dds)

 # the condition effect for genotype I (the main effect)
 results(dds, contrast=c("condition","B","A"))

 # the condition effect for genotype II
 # this is, by definition, the main effect *plus* the interaction term
 # (the extra condition effect in genotype II compared to genotype I).
 results(dds, list( c("condition_B_vs_A","genotypeII.conditionB") ))

 # the interaction term, answering: is the condition effect *different* across genotypes?
 results(dds, name="genotypeII.conditionB")