[BioC] DESeqs two-factor two-level, interaction is interested

[Michael Love]

hi Shawn,

Let's keep all replies on the Bioconductor mailing list, so other
users can follow the thread.


On Sat, May 17, 2014 at 11:32 AM, Shucong Li <lishucongnn at gmail.com> wrote:
> Hello Mike,
>
> Thank you for reply on weekend.
>
> Sorry for not make a part of my question properly. What I really meant is to compare :
>
> "Factor A level 1" vs "Factor A  level 2"  within   "Fact B level 1"
> "Factor A level 1" vs "Factor A  level 2"  within   "Fact B level 2"

I don't understand what you mean by these comparisons, specifically
the 'within' part.

There are four terms in this model. I will write in the log2 scale, so
terms are additive (this is multiplication on the scale of counts):

log2(counts) = intercept + treatmenttreat + dayB + dayB.treatmenttreat

the treatmenttreat term is the difference between treatment and
control for day A

the dayB term is the difference between day B over day A for the control group

the final term accounts for the case that the (day B and treatment)
group is not merely the sum of the previous two terms.

If you want to compare treatment vs control for day B, that term is
the treatmenttreat term plus the dayB.treatmenttreat term. This can be
pulled out easily with a numeric contrast (see the help for ?results
for more details).

results(dds, contrast=c(0,1,0,1))

where the numbers correspond to the order in resultsNames(dds)

I'd recommend looking over some introductory material to interaction
models (e.g. http://en.wikipedia.org/wiki/Interaction_(statistics) ).

Mike

>
> Is it possible to do so?  I checked all vignettes and manual of DESeq2 and still could not figure this out.
>
> Shawn
>
> On May 17, 2014, at 9:58 AM, Michael Love <michaelisaiahlove at gmail.com> wrote:
>
>> hi Shawn,
>>
>>
>> On Sat, May 17, 2014 at 2:00 AM, Shawn [guest] <guest at bioconductor.org> wrote:
>>> Hello Mike,
>>>
>>> Please allow me ask a basic question.  What does 'log2FoldChange' in the results of  DESeq2 analysis really mean for the interaction of a two-factor two-level design?
>>
>> The meaning is the same as for other linear models. The interaction
>> term for the generalized linear model tests if the effect of both day
>> and treatment is not simply multiplicative (or additive in the log2
>> space). If this term is significantly non-zero (the default test in
>> DESeq2), then being in the group: (day=B and treatment=treat) is not
>> simply the product of the day=B fold change and the treatment=treat
>> fold change.
>>
>>
>>> Is it possible to compare 'Factor A level 1' to ' Factor A level 2' or other similar comparison?
>>
>> For example, the day B over A effect:
>>
>> results(dds, contrast=c("day","B","A"))
>>
>> or equivalently
>>
>> results(dds, name="day_B_vs_A")
>>
>> Check the help for ?results for more examples.
>>
>> Mike
>>
>>>
>>> Here are the part of codes I used:
>>>
>>>
>>> dds <- phyloseq_to_deseq2(phyloseq.obj, design=~ Treatment*Day)
>>>
>>> colData(dds)$Treatment<- factor(colData(dds)$Treatment,levels=c("Control","Treat"));
>>> colData(dds)$Day<- factor(colData(dds)$Day,levels=c("A","B"))
>>>
>>> dds$Treatment<- relevel(dds$Treatment, "Control")
>>> dds$Day<- relevel(dds$Day, "A")
>>>
>>> dds <- DESeq(dds, fitType="local",betaPrior=FALSE)
>>>
>>> Shawn
>>>
>>>
>>> -- output of sessionInfo():
>>>
>>> sessionInfo()
>>> R version 3.1.0 (2014-04-10)
>>> Platform: x86_64-apple-darwin13.1.0 (64-bit)
>>>
>>> locale:
>>> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
>>>
>>> attached base packages:
>>> [1] parallel  grid      splines   stats     graphics  grDevices utils     datasets  methods   base
>>>
>>> other attached packages:
>>> [1] ecodist_1.2.9           Biostrings_2.32.0       doParallel_1.0.8        foreach_1.4.2
>>> [5] iterators_1.0.7         metagenomeSeq_1.6.0     gplots_2.13.0           limma_3.20.1
>>> [9] Biobase_2.24.0          DESeq2_1.4.5            GenomicRanges_1.16.3    GenomeInfoDb_1.0.2
>>> [13] RcppArmadillo_0.4.300.0 Rcpp_0.11.1             XVector_0.4.0           IRanges_1.22.6
>>> [17] BiocGenerics_0.10.0     locfit_1.5-9.1          phangorn_1.99-7         genefilter_1.46.1
>>> [21] adephylo_1.1-6          scatterplot3d_0.3-35    analogue_0.12-0         rgl_0.93.996
>>> [25] princurve_1.1-12        labdsv_1.6-1            mgcv_1.7-29             indicspecies_1.7.1
>>> [29] biom_0.3.13             ggplot2_0.9.3.1         plyr_1.8.1              phyloseq_1.9.2
>>> [33] pamr_1.54.1             cluster_1.15.2          survival_2.37-7         vegan_2.0-10
>>> [37] lattice_0.20-29         permute_0.8-3           RColorBrewer_1.0-5      matrixStats_0.8.14
>>> [41] MASS_7.3-33             ape_3.1-1               ade4_1.6-2              nlme_3.1-117
>>>
>>> loaded via a namespace (and not attached):
>>> [1] adegenet_1.4-1       annotate_1.42.0      AnnotationDbi_1.26.0 bitops_1.0-6
>>> [5] brglm_0.5-9          caTools_1.17         codetools_0.2-8      colorspace_1.2-4
>>> [9] data.table_1.9.2     DBI_0.2-7            digest_0.6.4         fastmatch_1.0-4
>>> [13] gdata_2.13.3         geneplotter_1.42.0   gtable_0.1.2         gtools_3.4.0
>>> [17] httpuv_1.3.0         igraph_0.7.0         KernSmooth_2.23-12   Matrix_1.1-3
>>> [21] multtest_2.20.0      munsell_0.4.2        phylobase_0.6.8      proto_0.3-10
>>> [25] R.methodsS3_1.6.1    reshape2_1.4         RJSONIO_1.2-0.2      RSQLite_0.11.4
>>> [29] scales_0.2.4         shiny_0.9.1          stats4_3.1.0         stringr_0.6.2
>>> [33] tools_3.1.0          XML_3.98-1.1         xtable_1.7-3         zlibbioc_1.10.0
>>>
>>>
>>> --
>>> Sent via the guest posting facility at bioconductor.org.
>