Hi,
I have an experiment with both age and experimental condition. The only age samples are young and old. The only experimental conditions are control and conditioned.
I apply the following code:
dds <-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable,directory = Counts, design = ~ hearing + age + age:hearing)
dds <- DESeq(dds)
mod_mat <- model.matrix(design(dds), colData(dds))
resultsNames(dds)
keep <- rowSums(counts(dds)) >= 1
dds <- dds[keep,]
res <- results(dds, name = "age_Young_vs_Old")
summary(res)
res <- results(dds, name = "hearing_Deaf_vs_Control")
summary(res)
int <- results(dds, name = "ageYoung.hearingDeaf")
summary(int)
young_control <- colMeans(mod_mat[dds$age == "Young" & dds$hearing == "Control", ])
young_deaf <- colMeans(mod_mat[dds$age == "Young" & dds$hearing == "Deaf", ])
old_control <- colMeans(mod_mat[dds$age == "Old" & dds$hearing == "Control", ])
old_deaf <- colMeans(mod_mat[dds$age == "Old" & dds$hearing == "Deaf", ])
resyc_vs_oc <- results(dds, contrast = young_control - old_control)
summary(resyc_vs_oc, alpha = 0.05)
resyc_vs_yd <- results(dds, contrast = young_control - young_deaf)
summary(resyc_vs_yd, alpha = 0.05)
resoc_vs_od <- results(dds, contrast = old_control - old_deaf)
summary(resoc_vs_od, alpha = 0.05)
resyd_vs_od <- results(dds, contrast = young_deaf - old_deaf)
summary(resyd_vs_od, alpha = 0.05)
My problem is with the model. If I change it from:
design = ~ hearing + age + age:hearing
to
design = ~ age + hearing + age:hearing
My standard age and condition results are the same, but the number of differentially expressed genes in the contrasts switches around (although the numbers are the same). So for example old control vs old deaf gives 1000 DEGs in the first model. When the model is switched around this becomes 9. Meanwhile Young control vs old control has 9 in the first model and this changes to 1000 in the second model.
Am I performing the contrasts wrong? If which factor (age or hearing) should I start the model with?
Many thanks,
Christian
This code is not reproducible. It does not contain any test data, nor the definition of what
mod_matis (I guess it's simply a model matrix but who knows...) and also not the part where you switch design and how you let DESeq2 know about it. Please try to make a MRE withmakeExampleDESeqDataSet()(usingset.seed(1)for reproducibility) and then evaluate this with a reprex and share that.But actualy...why not combining hearing and age into a single factor and then use contrasts as in the vignette rather than this custom way of extracting it from the model matrix? See:
Apologies.
I made a typo which excluded resetting the mod_mat part each time. The model works excluding this. I've edited this in the script above in case it's of any use to anyone going forward. Thank you for your help ATpoint