Question: DESeq2: interaction controlling for genetic effects
0
gravatar for rachelwright8
4.6 years ago by
United States
rachelwright810 wrote:

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 • 1.2k views
ADD COMMENTlink modified 4.6 years ago by Michael Love24k • written 4.6 years ago by rachelwright810
Answer: DESeq2: interaction controlling for genetic effects
0
gravatar for Michael Love
4.6 years ago by
Michael Love24k
United States
Michael Love24k wrote:

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 COMMENTlink written 4.6 years ago by Michael Love24k

Wald it is.  Thanks, Michael!

ADD REPLYlink written 4.6 years ago by rachelwright810

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

ADD REPLYlink written 4.6 years ago by Michael Love24k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 309 users visited in the last hour