I have a RNA-Seq time series with two conditions (SD and LD) and 3 time points (0, 2, and 4). There are three biological replicates per sample. I have used DESeq2 to do a likelihood ratio test on the following designs: 

full ~ condition + time + condition:time
reduced ~ condition + time

I understand that this tests for the effects of the two conditions on the pattern of expression which is what I want. Here are the names of my coefficients:

> resultsNames(dds_test)
[1] "Intercept"          "condition_LD_vs_SD" "time_2_vs_0"        "time_4_vs_0"        "conditionLD.time2"
[6] "conditionLD.time4"

I would also like to contrast the difference between the conditions at the three time points and between time points within each condition. I was hoping someone could confirm my interpretations and answer a question about comparing time points which are not the baseline.

So am I correct in thinking I can test the difference between the conditions at the three timepoints like so:

results(dds_test, name="condition_LD_vs_SD", test="Wald")
results(dds_test, contrast=list(c("condition_LD_vs_SD","conditionLD.time2")), test="Wald")
results(dds_test, contrast=list(c("condition_LD_vs_SD","conditionLD.time4")), test="Wald")

and between time points within each condition like so:

• Condition SD

results(dds_test, contrast=list(c("time_2_vs_0")), test="Wald") 
results(dds_test, contrast=list(c("time_4_vs_0")), test="Wald")

• Condition LD

results(dds_test, contrast=list(c("time_2_vs_0", "conditionLD.time2")), test="Wald") results(dds_test, contrast=list(c("time_4_vs_0", "conditionLD.time4")), test="Wald")

Is this correct and is there an easy way to calculate and test the fold changes between time points 2 and 4?

This is my sample table:

           sample condition time replicate
LD_ZT4_1 LD_ZT4_1        LD    4         1
LD_ZT4_2 LD_ZT4_2        LD    4         2
LD_ZT4_3 LD_ZT4_3        LD    4         3
SD_ZT4_1 SD_ZT4_1        SD    4         1
SD_ZT4_2 SD_ZT4_2        SD    4         2
SD_ZT4_3 SD_ZT4_3        SD    4         3
SD_ZT0_1 SD_ZT0_1        SD    0         1
SD_ZT0_2 SD_ZT0_2        SD    0         2
SD_ZT0_3 SD_ZT0_3        SD    0         3
SD_ZT2_1 SD_ZT2_1        SD    2         1
SD_ZT2_2 SD_ZT2_2        SD    2         2
SD_ZT2_3 SD_ZT2_3        SD    2         3
LD_ZT0_1 LD_ZT0_1        LD    0         1
LD_ZT0_2 LD_ZT0_2        LD    0         2
LD_ZT0_3 LD_ZT0_3        LD    0         3
LD_ZT2_1 LD_ZT2_1        LD    2         1
LD_ZT2_2 LD_ZT2_2        LD    2         2
LD_ZT2_3 LD_ZT2_3        LD    2         3

I'll jump to the answers. The differences between conditions at each time point are given by the following:

1. condition_LD_vs_SD - for time point 0
2. main effect + conditionLD.time2 - for time point 2
3. main effect + conditionLD.time4 - for time point 4

You can take a look at the diagram of interactions terms in the DESeq2 vignette to see how this is the case.

To build results tables for each of these, you should use:

results(dds, name="condition_LD_vs_SD", test="Wald")

And for the other two you need to do, e.g.:

results(dds, contrast=list(c("condition_LD_vs_SD","conditionLD.time4")), test="Wald")

You can also test between time points for each condition. For the reference level of condition (SD), these are given by the main effect terms. For condition LD, you add the main effect term and the relevant interaction term. So time 4 vs time 0 in LD is:

results(dds, contrast=list(c("time_4_vs_0", "conditionLD.time4")), test="Wald")

While for SD it is simply:

results(dds, name="time_4_vs_0", test="Wald")

For making all these results tables, please note that the correction for multiple testing by the results() function is only *within* a results table across genes. Usually only one or a few results tables are generated for an experiment, where we assume the results of all will be reported. If a user wants to make many results table, but then only selectively report gene sets from a subset of the tables, then we recommend that the user perform multiple test correction over all generated p-values. We don't have a function for this last step in DESeq2 but it can be done with R code.