Hi, 

I've done a model design and I hope someone can help out with my understanding of it!

I have an experimental setup that looks something like this:

3 Time points (0hrs, 6hrs, 12hrs)

3 Different Conditions (Treatments A, B and C)

So that makes 9 different combinations of time points and treatments, each is in triplicate. There are therefore27 Samples. 

My design formula is : ~ Treatment + TimePoint + Treatment:TimePoint

My current understanding is this will give me small pValues of Treatment-specific effects over time?

I wanted to further refine this and look at a specific treatment and how it differs between two time points. So I used the following line:

foo <- list(c(“TreatmentA.TimePoint12hrs"), c(“TreatmentA.TimePoint0hrs"))

resMFType <- results(dds, contrast=foo)

EDIT::

as.data.frame(colData(dds))

    TimePoint Treatment
X1        0hr         A
X2        0hr         B
X3        0hr         C
X4        0hr         A
X5        0hr         B
X6        0hr         C
X7        0hr         A
X8        0hr         B
X9        0hr         C
X10       6hr         A
X11       6hr         B
X12       6hr         C
X13       6hr         A
X14       6hr         B
X15       6hr         C
X16       6hr         A
X17       6hr         B
X18       6hr         C
X19      12hr         A
X20      12hr         B
X21      12hr         C
X22      12hr         A
X23      12hr         B
X24      12hr         C
X25      12hr         A
X26      12hr         B
X27      12hr         C

sessionInfo()

R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] DESeq2_1.6.2            RcppArmadillo_0.4.500.0 Rcpp_0.11.3             GenomicRanges_1.18.3   
[5] GenomeInfoDb_1.2.3      IRanges_2.0.0           S4Vectors_0.4.0         BiocGenerics_0.12.1    

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3      annotate_1.44.0      AnnotationDbi_1.28.1 base64enc_0.1-2     
 [5] BatchJobs_1.5        BBmisc_1.8           Biobase_2.26.0       BiocParallel_1.0.0  
 [9] brew_1.0-6           checkmate_1.5.0      cluster_1.15.3       codetools_0.2-9     
[13] colorspace_1.2-4     DBI_0.3.1            digest_0.6.4         fail_1.2            
[17] foreach_1.4.2        foreign_0.8-61       Formula_1.1-2        genefilter_1.48.1   
[21] geneplotter_1.44.0   ggplot2_1.0.0        grid_3.1.2           gtable_0.1.2        
[25] Hmisc_3.14-5         iterators_1.0.7      lattice_0.20-29      latticeExtra_0.6-26 
[29] locfit_1.5-9.1       MASS_7.3-35          munsell_0.4.2        nnet_7.3-8          
[33] plyr_1.8.1           proto_0.3-10         RColorBrewer_1.0-5   reshape2_1.4        
[37] rpart_4.1-8          RSQLite_1.0.0        scales_0.2.4         sendmailR_1.2-1     
[41] splines_3.1.2        stringr_0.6.2        survival_2.37-7      tools_3.1.2         
[45] XML_3.98-1.1         xtable_1.7-4         XVector_0.6.0  

Is this correct?

Thanks, 

"My current understanding is this will give me small pValues of Treatment-specific effects over time?"

The way to get a single p-value for each gene for any treatment-specific effects over time is to use the likelihood ratio test:

dds = DESeq(dds, test="LRT", reduced=~Treatment + TimePoint)
results(dds)

Where the results table will give a LRT statistic and p-value. Note that there are many LFCs involved in this likelihood ratio test (all the interaction effects between treatment and time), but the software has to pick just one to print in the table. It will just print the last one based on the ordering of levels. You can change this column with the 'name' argument. But you will notice that the statistic and p-value do not change, because these are based on all the interaction terms.

"I wanted to further refine this and look at a specific treatment and how it differs between two time points. So I used the following line...list(c(“TreatmentA.TimePoint12hrs"), c(“TreatmentA.TimePoint0hrs"))"

That line is correct if you used just DESeq(dds).

If you used the LRT instead, then the interaction terms for base levels are absorbed by the intercept. This means the Treatment A comparison would be:

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

...where specifying Wald means that we want p-values for this LFC only.

and for Treatment B this would be:

results(dds, contrast=list(c("TimePoint_12hr_vs_0hr","TreatmentB.TimePoint12hr")), test="Wald")

...which is the sum of the 12 vs 0 effect for the base level (treatment A) and the interaction effect of 12 vs 0 in treatment B.

Using DESeq(...,test="LRT") or just DESeq() are equally valid. The first is better to run if you want to test for any effect, and then make individual comparisons using results(). The second is perhaps easier syntax for comparing individual timepoints.