I am trying to take the "difference of differences" in contrasts from two factors (
group). We have male and female animals (
sex factor) that were untrained or trained for 1, 2, 4, or 8 weeks (
group factor, i.e., "control", "1w", "2w", "4w", "8w"). I want to know the difference in the training effect at each time point between the two sexes.
Conceptually, this is what I want at the 1-week time point, for example:
(male:1w - male:control) - (female:1w - female:control)
However, I can't figure out how to do this in DESeq2. Specifying this contrast in limma would be straightforward, but I want to use DESeq2 for consistency with other arms of the differential analysis. I have looked at other answers here and here as well as the Interactions section of the DESeq2 tutorial, but I am still perplexed.
I have tried specifying the contrast with an interaction term (i.e.
~ sex*group) as well as with a new factor that pastes together the
group factors. I feel like I got the closest with the interaction term:
dds <- makeExampleDESeqDataSet(n=50,m=100) dds$sex <- factor(rep(c("male","female"), each=25)) dds$group <- factor(rep(c("control","1w","2w","4w","8w"), 10)) # define reference levels dds$sex <- relevel(dds$sex, ref="male") dds$group <- relevel(dds$group, ref="control") design(dds) <- ~ group*sex dds <- DESeq(dds) resultsNames(dds)
 "Intercept" "group_1w_vs_control" "group_2w_vs_control" "group_4w_vs_control" "group_8w_vs_control" "sex_female_vs_male" "group1w.sexfemale" "group2w.sexfemale"  "group4w.sexfemale" "group8w.sexfemale"
From this post, I thought
group1w.sexfemale would test if the 1-week effect (i.e.,
1w - control) is different between females and males. However, the results aren't consistent with that. See below the plot of a gene that had a log fold-change of -23.2 and an adjusted p-value of 1e-4 for this specific contrast name. From the TMM-normalized data, I would expect a marginally significant result of roughly 3-fold.
I would appreciate any help!