DESeq2: interaction controlling for genetic effects
1
0
Entering edit mode
@rachelwright8-7146
Last seen 9.6 years ago
United States

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   

 

deseq interactions • 2.1k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 17 hours ago
United States

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).

ADD COMMENT
0
Entering edit mode

Wald it is.  Thanks, Michael!

ADD REPLY
0
Entering edit mode

also, I'd recommend updating to v1.6 which was released in October if possible. http://bioconductor.org/install/#update-bioconductor-packages

ADD REPLY

Login before adding your answer.

Traffic: 615 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6