DESeq2 design formula: how to use the opposite interaction term?

Question

[salamandra]

Hi, a question regarding deseq2 design formula.

I have two variables time and treatment and I know how to get genes which fold-change of treated vs untreated varies with time:

Table <- data.frame(sampleName = sampleNames, fileName = sampleFiles, time = time, Treatment = Treatment, condition = condition)
Table
  sampleName         fileName time Treatment condition
1       A1_1 1A_ATCACG_counts   t0     Treat   Treatt0
2       A1_2 2A_CAGATC_counts   t0     Treat   Treatt0
3       B1_1 2B_ACTTGA_counts   t0   Untreat Untreatt0
4       B1_2 3B_CGATGT_counts   t0   Untreat Untreatt0
5       A2_1 2C_GATCAG_counts   t1     Treat   Treatt1
6       A2_2 3C_TTAGGC_counts   t1     Treat   Treatt1
7       B2_1 2E_GGCTAC_counts   t1   Untreat Untreatt1
8       B2_2 3F_GCCAAT_counts   t1   Untreat Untreatt1

dds3 <- DESeqDataSetFromHTSeqCount(sampleTable = Table, design= ~ Treatment + time + Treatment:time)
ddsHTSeq3 <- DESeq(dds3)
resultsNames(ddsHTSeq3)
"Intercept"                  "Treatment_Treat_vs_Untreat" "time_t1_vs_t0"              "TreatmentTreat.timet1"
resHTSeqb <- results(ddsHTSeq3, name = "TreatmentTreat.timet1")

But how do to get genes which fold-change of time1 vs time0 varies with treatment? cause

results(ddsHTSeq3, name = "timet1.TreatmentTreat")

doesn't work

Answer 1

[Kevin Blighe]

If you had just specified this:

DESeqDataSetFromHTSeqCount(sampleTable = Table, design= ~ Treatment + time)

...then, the main time effect would look at time while controlling for Treatment. This is what you want, I believe.

Kevin

[salamandra]

In the example we test whether treatment depends on the time: results(ddsHTSeq3, name = "TreatmentTreat.timet1") what I wanted to know is how to test whether time depends on the treatment, that is if time1 vs time0 fold-change is the same for treated and untrated samples or if it changes.

[Kevin Blighe]

In the example, 'TreatmentTreat.timet1' relates to the time effect for TreatmentTreat Versus TreatmentUnTreat. That is, it tests whether the effect of time is different between treated and untreated. I feel that this indirectly answers your question.

Another way to look at it:

In your current example, 'time_t1_vs_t0' relates to timet1 versus timet0 for TreatmentUntreat. If you re-level Treatment to have Treat as the reference level, then 'time_t1_vs_t0' would relate to timet1 versus timet0 for TreatmentTreat. This would tell you if treated and untreated change in the same way independently at the 2 time-points.

-----------------------------------------------------

Also be sure to read the DESeq2 vignette, in particular:

The key point to remember about designs with interaction terms is that, unlike for a design ~genotype + condition, where the condition effect represents the overall effect controlling for differences due to genotype, by adding genotype:condition, the main condition effect only represents the effect of condition for the reference level of genotype (I, or whichever level was defined by the user as the reference level). The interaction terms genotypeII.conditionB and genotypeIII.conditionB give the difference between the condition effect for a given genotype and the condition effect for the reference genotype.

[salamandra]

Thank you for your time. From what you pointed to in theDESeq2 vignette I intrepert as time_t1_vs_t0 is time t1 vs t0 for the refference level (Untreat)

Looking at the example in the vignette (easier to understand). To test if the condition effect (A vs B) varies with genotype (genotype II vs I) we would do:

dds3 <- DESeqDataSetFromHTSeqCount(sampleTable = Table, design= ~ genotype + condition + genotype:condition)
ddsHTSeq3 <- DESeq(dds3)
resultsNames(ddsHTSeq3)    
[1] "Intercept"              "genotype_II_vs_I"       "genotype_III_vs_I"      "condition_B_vs_A"       "genotypeII.conditionB"  "genotypeIII.conditionB"
results(ddsHTSeq3, name="genotypeII.conditionB")

To test if the genotype effect (genotype II vs I) varies with condition effect (A vs B) I believe we would do:

dds4 <- DESeqDataSetFromHTSeqCount(sampleTable = Table, design= ~ condition + genotype + condition:genotype)
ddsHTSeq4 <- DESeq(dds4)
resultsNames(ddsHTSeq4)
[1] "Intercept"              "condition_B_vs_A"       "genotype_II_vs_I"       "genotype_III_vs_I"      "conditionB.genotypeII"  "conditionB.genotypeIII"
results(ddsHTSeq4, name="conditionB.genotypeII")

But then noticed that results(ddsHTSeq3, name="genotypeII.conditionB") is same as results(ddsHTSeq4, name="conditionB.genotypeII") Maybe, both tests are actually testing for the same...

[Kevin Blighe]

Yes, they would be looking at the same thing. Changing the order of the interaction term will not make any difference (e.g. Condition:Time or Time:Condition). What you need to do is re-level the factors in your metadata in order to change the category that is the reference level.

Take a look at Devon's answer, here: A: DESeq2 doesn't compare the right ctrl vs treatment?

[salamandra]

sorry, i'm not understanding we have a table like this:

Table
   sampleName         fileName genotype condition group
1       AI_s1 A1_ATCACG_counts        I         A    AI
2       AI_s2 A2_TGACCA_counts        I         A    AI
3      AII_s1 A3_CAGATC_counts       II         A   AII
4      AII_s2 A4_TAGCTT_counts       II         A   AII
5     AIII_s1 B1_CGATGT_counts      III         A  AIII
6     AIII_s2 B2_ACAGTG_counts      III         A  AIII
7       BI_s1 B3_ACTTGA_counts        I         B    BI
8       BI_s2 B4_GGCTAC_counts        I         B    BI
9      BII_s1 C1_TTAGGC_counts       II         B   BII
10     BII_s2 C2_GCCAAT_counts       II         B   BII
11    BIII_s1 C3_GATCAG_counts      III         B  BIII
12    BIII_s2 C4_CTTGTA_counts      III         B  BIII

and we do:

dds3 <- DESeqDataSetFromHTSeqCount(sampleTable = Table, design= ~ genotype + condition + genotype:condition)
dds3 <- dds3[rowSums(counts(dds3)) > 1, ]
ddsHTSeq3 <- DESeq(dds3)
results(ddsHTSeq3, name="genotypeII.conditionB")
levels(ddsHTSeq3$genotype) 
[1] "I"   "II"  "III"    # I comes first so it's the reference level
levels(ddsHTSeq3$condition) 
[1] "A" "B"    # A comes first so it's the reference level

if we relevel like you suggested, for e.g. like this:

dds3$genotype <- relevel(dds3$genotype, "II")
dds3$condition <- relevel(dds3$condition, "B")
resultsNames(ddsHTSeq3)

and when do:

ddsHTSeq3 <- DESeq(dds3)
resultsNames(ddsHTSeq3)
[1] "Intercept"              "genotype_I_vs_II"       "genotype_III_vs_II"     "condition_A_vs_B"       "genotypeI.conditionA"   "genotypeIII.conditionA"

There is no term "conditionB.genotypeII" available. The only thing relevel seems to do is change reference level in each factor.

what I wanted was to know how condition effect (B vs A with A as ref level) varies with genotype (II vs I with I as ref level) I do not want to know how genotype effect (II vs I with I as ref level) varies with condition (B vs A with A as ref level)

[Kevin Blighe]

Yes, that is all that relevel does, but that is key to know.

For what you want, if genotype 'I" is the reference level:

• genotypeII.conditionB = gives the condition effect for genotype II vs I
• genotypeIII.conditionB = gives the condition effect for genotype III vs I

If you want something more complex, then consider merging multiple factors into the same factor (e.g. your group factor)). If you are still not content with my explanation, then please try the Bioconductor support forum

[Kevin Blighe]

For example, if you abandon the interaction terms in your model and instead just use ~ group, you could then generate different results with pairwise comparisons:

results(dds, contrast=c("group", "IIB", "IA"))
results(dds, contrast=c("group", "IIIB", "IA"))

Please take time to read the vignette for interactions: https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions