Hi all,

I'm reanalyzing a gene expression data set with DESeq2 that I had already analyzed in DESeq and have some questions about how which models would be most appropriate.  My sample data is below.  I have a virus treatment (virus "V" or control "C) and a heat treatment (heat "H" or ambient "A").  There is also a genetic effect (gen, a combination of dam and sire) that I'd like to control for.  

virus    heat     gen

C         A         1C

C         H         1C

V         A         1C

C         A         1D

V         H         1D

C         H         2A

V         H         2A

C         A         2C

C         H         2C

V         A         2C

V         H         2C

C         A         2D

C         H         2D

V         A         2D

V         H         2D

I'm wanting to test for the effects of (a) virus alone, (b) heat alone, and (c) the interaction of virus and heat treatments.  I've tried using both likelihood ratio tests and Wald tests with some sort of result, but I'm not sure which would be most appropriate.  

Using the LRTs I ran the models three times:

dds.h <- DESeq(dds.h, test="LRT", full=~gen+heat, reduced=~gen)

dds.v <- DESeq(dds.v, test="LRT", full=~gen+virus, reduced=~gen)

dds.i <- DESeq(dds.i, test="LRT", full=~gen+virus*heat, reduced=~sire+virus+heat)

res.h<-results(dds.h)

res.v<-results(dds.v)

res.i<-results(dds.i)

Using the Wald test I only ran the model one time:

design <- ~ gen + virus +heat + virus:heat

dds<-DESeqDataSetFromMatrix(countData=countData, colData=colData, design=design)

and specified contrasts to get pvalues:

res.h<-results(dds,contrast=c("heat", "A", "H"))

res.v<-results(dds,contrast=c("virus", "C", "V"))

res.i<-results(dds)

Both of these methods seem to work just fine in the sense that I don't receive error messages and they return largely the same DEGs as well.  Which method (LRT/Wald) do you think would be most appropriate?

Thank you!

Rachel

#############

> sessionInfo()

R version 3.1.1 (2014-07-10)

Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] grid      parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

 [1] VennDiagram_1.6.9         gplots_2.14.2             RColorBrewer_1.0-5        DESeq2_1.4.5              RcppArmadillo_0.4.450.1.0

 [6] Rcpp_0.11.3               GenomicRanges_1.16.4      GenomeInfoDb_1.0.2        IRanges_1.22.10           BiocGenerics_0.10.0      

loaded via a namespace (and not attached):

 [1] annotate_1.42.1      AnnotationDbi_1.26.1 Biobase_2.24.0       bitops_1.0-6         caTools_1.17.1       DBI_0.3.1           

 [7] gdata_2.13.3         genefilter_1.46.1    geneplotter_1.42.0   gtools_3.4.1         KernSmooth_2.23-13   lattice_0.20-29     

[13] locfit_1.5-9.1       RSQLite_0.11.4       splines_3.1.1        stats4_3.1.1         survival_2.37-7      tools_3.1.1         

[19] XML_3.98-1.1         xtable_1.7-4         XVector_0.4.0   

hi Rachel,

I would recommend using the Wald test, as it simplifies the results extraction and reduces the number of models to fit. Actually the LRT calls above are not good, as you also need to change the design(dds) before running DESeq(), rather than just changing 'full'. This is mentioned in the help for ?DESeq, but I could see how one would miss it: "full : the full model formula, this should be the formula in design(object)".  This is because the dispersion steps use design(dds), and this should be the same as the design used by the final testing step. I should add at the least a warning when the design in full is not the same as the design(dds).