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
•
written
by
salamandra
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
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.
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.
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...
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?
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)
Yes, that is all that relevel does, but that is key to know.
For what you want, if genotype 'I" is the reference level:
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
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
Log in to answer this question.