Hello,

I have a complex time-series RNA-seq expression dataset which I'm currently trying to analyse using DESeq2.

I have two plant species, grown under two nitrate levels along six days of a heat wave experiment (a day before the heat exposure, three days of heat exposure, and two days of recovery). I also have the data from the last day, of plants (both species, with either nutrient levels) that were not exposed to heat at all (as a second control), so I tried analysing the data with and without these second control samples.   

Here is my coldata table (I'm presenting only the samples of species "At" to reduce the text, but it's the same for species "Ah"):

As one species is more tolerant to the treatment than the other, my main biological questions are

1. Which genes expressed differently between the two species under different nutrient treatments (species:nitrate interaction).

2. Which genes are differently regulated during the different stages (control, heat and recovery), hence the species:nitrate:time interaction with different contrasts.

However, when I'm running the test, it seems that some of the resultsNames levels are missing. 

Here is my code, when I'm ignoring the last day control samples (excluded from both counts and coldata tables):

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~species + nitrate + time_point + species:time_point + nitrate:time_point + species:nitrate + species:nitrate:time_point)

dds <- DESeq(dds)

## As a result I get:

resultsNames(dds)

[1] "Intercept"                         "species_At_vs_Ah"                 
 [3] "nitrate_low_vs_high"               "time_point_2_vs_1"                
 [5] "time_point_3_vs_2"                 "time_point_4_vs_1"                
 [7] "time_point_5_vs_2"                 "time_point_6_vs_1"               
 [9] "speciesAt.time_point2"             "speciesAt.time_point3"            
[11] "speciesAt.time_point4"             "speciesAt.time_point5"            
[13] "speciesAt.time_point6"            "nitratelow.time_point2"           
[15] "nitratelow.time_point3"            "nitratelow.time_point4"           
[17] "nitratelow.time_point5"            "nitratelow.time_point6"          
[19] "speciesAt.nitratelow.time_point1"  "speciesAt.nitratelow.time_point2" 
[21] "speciesAt.nitratelow.time_point3"  "speciesAt.nitratelow.time_point4" 
[23] "speciesAt.nitratelow.time_point5"  "speciesAt.nitratelow.time_point6"

I was expecting to get also interaction results for the "Ah" species and high nitrate levels, e.g. speciesAh.nitratelow.time_point2  and speciesAt.nitratehigh.time_point2, to contrast the two species. Am I missing something here?

I also tried to combine the species and nitrate factors to one species_nitrate factor, and running the test with the following design:

design(dds)
~species_nitrate + time_point + species_nitrate:time_point

## and then I got:

resultsNames(dds)

[1] "Intercept"                           "species_nitrate_Ah.low_vs_Ah.high"  
 [3] "species_nitrate_At.high_vs_Ah.high"  "species_nitrate_At.low_vs_Ah.high"  
 [5] "time_point_5_vs_2"                   "time_point_7_vs_2"                  
 [7] "time_point_8_vs_2"                   "time_point_9_vs_2"                  
 [9] "time_point_10_vs_2"                  "species_nitrateAh.low.time_point5"  
[11] "species_nitrateAt.high.time_point5"  "species_nitrateAt.low.time_point5"  
[13] "species_nitrateAh.low.time_point7"   "species_nitrateAt.high.time_point7" 
[15] "species_nitrateAt.low.time_point7"   "species_nitrateAh.low.time_point8"  
[17] "species_nitrateAt.high.time_point8"  "species_nitrateAt.low.time_point8"  
[19] "species_nitrateAh.low.time_point9"   "species_nitrateAt.high.time_point9" 
[21] "species_nitrateAt.low.time_point9"   "species_nitrateAh.low.time_point10" 
[23] "species_nitrateAt.high.time_point10" "species_nitrateAt.low.time_point10"

Where here, I'm missing the species_nitrate_At.low_vs_At.high and species_nitrate_At.low_vs_Ah.low results, as well as the results for the Ah.high interaction with time.

Can you please help me sort it out?

best,

Gil

R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)
DESeq2_1.12.4

hi,

The levels are not missing. When you have interaction terms in a formula and main effects, you will not have interaction terms for the reference level of the factors (e.g. species=Ah, nitrate=high, time point=2). The interaction terms are *additional* effects for the non-reference levels beyond the effect seen for the reference level. This is diagramed in the DESeq2 vignette section on interactions.

Given that you have a fairly complex experimental design with potential interactions between three factors, I'd recommend that you collaborate with a statistician to help with a design formula and how to extract the comparisons you'd like to make. The choice of design, and so the terms that show up in the model, is not one-to-one with the layout of the sample table, but involves some discussion and assumptions.