[BioC] DESeqs two-factor two-level, interaction is interested
Shucong Li
lishucongnn at gmail.com
Sun May 18 06:32:32 CEST 2014
Hi Mike,
Thank you. With your help, I finally figured it out, sort of :-)
Enjoy the rest of your weekend.
Shawn
On May 17, 2014, at 11:07 AM, Michael Love <michaelisaiahlove at gmail.com> wrote:
> 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.
>>
More information about the Bioconductor
mailing list