Hello,
I'm a bit of a novice to DE analyses, or perhaps am misunderstanding the negative binomial model and what information it can tell us, but I'm confused on how to extract some bits of important data from my analyses. I have RNASeq count data and a design with three factors: age at 5 levels, sex at 2 levels, and condition at 2 levels. My design is set as ~0+age*sex*condition. code sample:
age <- (rep(c("P1","P2","P3","P4","P5"),each=8))
sex <- (rep(c("M","M","M","M","F","F","F","F"), length=40))
condition <- (rep(c("Exp","Exp","Base","Base"),length=40))
x <- as.data.frame(cbind(age,sex,condition))
mm<-model.matrix(~0+age*sex*condition, data=x)
I'm interested in looking at the intersection of age and condition, both with and without sex included. That is, for example: first, what is the difference between conditionExp and conditionBase at age P5 across both sexes; then, asking if there is a sex difference in that interaction.
I understand how to extract the additive effect of M on the condition:age interaction, but I am unsure how to ask the first question, to use contrasts to extract the effect of condition:age across both sexes, if that's clear. We don't conceptualize Female as a baseline to Male, and I understand the utility of needing a reference level, but we would like to first interrogate the other factors before asking if there is an additional sex difference.
For the effect of condition:age:male (i.e. what is the Exp/Base difference in P5 Males), my coding would be originally the main effect terms for ageP5, conditionExp, sexM, and their two-way and three-way interaction terms as the numerator, divided by ageP5, sexM, ageP5:sexM as the denominator, which reduces to: contrast=list(c("ageP5:conditionExp","sexM:conditionExp","ageP5:sexM:conditionExp")
My initial thoughts on coding the age:condition interaction without respect to sex was to use P5, Exp, P5:Exp as the numerator and P5 as the denominator, reducing to:
contrast=list(c("conditionExp","ageP5:conditionExp")
However, I believe this would in fact tell me the condition:age effect for the reference level of sex (F), not for all levels together. Is this correct? And if so, how might I extract the condition:age effect without regard to sex? Is it possible with my current design?
I've read through the vignette and understand that the interaction terms should be interpreted as the interaction effect of those 2/3 terms when any terms left out are at their reference level. Therefore, my second contrast example above would tell me Exp/Base within P90 when sex is set to F (reference), not Exp/Base for all P90, correct?
One suggestion I had received from a statistician colleague was, to recode the model matrix to represent -0.5, 0.5 instead of 0,1 so that with a single call of an interaction term we might be able to extract the contrast of interest. I'm not necessarily sure if this is advisable for ngs data, or even possible, with DESeq2? She was unfamiliar with the package and its capabilities, but suggested this was what she might do for a singular dataset (i.e. one set of counts instead of many genes' worth of counts).
Thank you in advance for any help! Please let me know if I can be more descriptive or provide additional helpful detail.
