DESeq factor order
Entering edit mode
4 months ago
Chris ▴ 10


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))


keep <- rowSums(counts(dds)) >= 1

dds <- dds[keep,]

res <- results(dds, name = "age_Young_vs_Old")


res <- results(dds, name = "hearing_Deaf_vs_Control")


int <- results(dds, name = "ageYoung.hearingDeaf")


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


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,


DESeq2 • 369 views
Entering edit mode

This code is not reproducible. It does not contain any test data, nor the definition of what mod_mat is (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 with makeExampleDESeqDataSet() (using set.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:

Entering edit mode


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


Login before adding your answer.

Traffic: 3185 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6