Hello Mike,
Follow up your post here
C: DESeq2 paired and multi factor comparison
I'm working on a project with a similar design.
I'm interested to know the effect of "Day" across both groups ('Treatment' and 'Control') of animals.
Thanks in advance.
Here are my experimental design. 'Calf_nested' was created following the trick from the edgeR user guide, as you suggested in the post DESeq2 paired multifactor test.
Treatment | Calf | Calf_nested | Day |
Control | 30 | 1 | d36 |
Control | 30 | 1 | d54 |
Control | 31 | 2 | d36 |
Control | 31 | 2 | d54 |
Control | 36 | 3 | d36 |
Control | 36 | 3 | d54 |
Control | 37 | 4 | d36 |
Control | 37 | 4 | d54 |
Control | 38 | 5 | d36 |
Control | 38 | 5 | d54 |
Treatment | 32 | 1 | d36 |
Treatment | 32 | 1 | d54 |
Treatment | 33 | 2 | d36 |
Treatment | 33 | 2 | d54 |
Treatment | 34 | 3 | d36 |
Treatment | 34 | 3 | d54 |
Treatment | 35 | 4 | d36 |
Treatment | 35 | 4 | d54 |
Treatment | 39 | 5 | d36 |
Treatment | 39 | 5 | d54 |
sessionInfo()
R version 3.1.3 (2015-03-09)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)
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] stats4 parallel grid splines stats graphics grDevices utils datasets methods base
other attached packages:
[1] ecodist_1.2.9 Biostrings_2.34.1 doParallel_1.0.8 foreach_1.4.2 iterators_1.0.7 metagenomeSeq_1.8.3
[7] interactiveDisplay_1.4.1 limma_3.22.7 Biobase_2.26.0 DESeq2_1.6.3 GenomicRanges_1.18.4 GenomeInfoDb_1.2.5
[13] RcppArmadillo_0.5.000.0 Rcpp_0.11.5 XVector_0.6.0 IRanges_2.0.1 S4Vectors_0.4.0 BiocGenerics_0.12.1
[19] locfit_1.5-9.1 genefilter_1.48.1 adephylo_1.1-6 scatterplot3d_0.3-35 analogue_0.16-0 rgl_0.93.1098
[25] princurve_1.1-12 labdsv_1.6-1 mgcv_1.8-6 biom_0.3.12 ggplot2_1.0.1 reshape2_1.4.1
[31] plyr_1.8.2 phyloseq_1.10.0 pamr_1.55 cluster_2.0.1 survival_2.38-1 vegan_2.2-1
[37] lattice_0.20-31 permute_0.8-3 RColorBrewer_1.1-2 matrixStats_0.14.0 MASS_7.3-40 ape_3.2
[43] ade4_1.7-2 nlme_3.1-120
loaded via a namespace (and not attached):
[1] acepack_1.3-3.3 adegenet_1.4-2 annotate_1.44.0 AnnotationDbi_1.28.2 base64enc_0.1-2
[6] BatchJobs_1.6 BBmisc_1.9 BiocParallel_1.0.3 bitops_1.0-6 brew_1.0-6
[11] brglm_0.5-9 Category_2.32.0 caTools_1.17.1 checkmate_1.5.2 chron_2.3-45
[16] codetools_0.2-11 colorspace_1.2-6 data.table_1.9.4 DBI_0.3.1 digest_0.6.8
[21] fail_1.2 foreign_0.8-63 Formula_1.2-1 gdata_2.13.3 geneplotter_1.44.0
[26] gplots_2.16.0 graph_1.44.1 gridSVG_1.4-3 GSEABase_1.28.0 gtable_0.1.2
[31] gtools_3.4.2 Hmisc_3.15-0 htmltools_0.2.6 httpuv_1.3.2 igraph_0.7.1
[36] interactiveDisplayBase_1.4.0 KernSmooth_2.23-14 latticeExtra_0.6-26 Matrix_1.2-0 mime_0.3
[41] multtest_2.22.0 munsell_0.4.2 nnet_7.3-9 phylobase_0.6.8 proto_0.3-10
[46] R6_2.0.1 RBGL_1.42.0 RJSONIO_1.3-0 rpart_4.1-9 RSQLite_1.0.0
[51] scales_0.2.4 sendmailR_1.2-1 shiny_0.11.1 stringr_0.6.2 tools_3.1.3
[56] XML_3.98-1.1 xtable_1.7-4 zlibbioc_1.12.0
Here are the codes I tried:
dds <- phyloseq_to_deseq2(phyloseq.obj,design= ~ Treatment*Day)
design(dds) <- ~Treatment+Treatment:Calf_nested+Treatment:Day
design(dds)
colData(dds)$Calf_nested<- factor(colData(dds)$Calf_nested)
colData(dds)$Treatment<- factor(colData(dds)$Treatment,levels=c("Control","Treatment"));
colData(dds)$Day<- factor(colData(dds)$Day,levels=c("d36","d54"))
dds$Treatment<- relevel(dds$Treatment, "Control")
dds$Day<- relevel(dds$Day,"d36")
dds <- removeResults(dds)
metaDF <- data.frame(sample_data(phyloseq.obj))
mm1 <- model.matrix(~ Treatment + Treatment:Calf_nested + Treatment:Day, metaDF)
mm0 <- model.matrix(~ Treatment + Treatment:Calf_nested, metaDF)
design(dds) <- ~ Treatment + Treatment:Calf_nested + Treatment:Day
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds, modelMatrix=mm1)
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds, modelMatrix=mm1)
dds <- nbinomLRT(dds, full=mm1, reduced=mm0)
### To test for taxa which respond to day effect in control group:
resCtrst_day_in_control <- results(dds,name="TreatmentControl.Dayd54",test="Wald")
### To test for taxa which respond to day effect in treatment group:
resCtrst_day_in_treatment <- results(dds,name="TreatmentTreatment.Dayd54", test="Wald")
### To test for taxa where the day effect differs between two groups:
resCtrst_interaction = results(dds, contrast=list("TreatmentTreatment.Dayd54","TreatmentControl.Dayd54"), test="Wald")
Thanks in advance.
Shucong
hi Shucong,
I'm sticking to my answer above. You don't have individual calves across treatment. You can't estimate the calf effect and the treatment effect because they are collinear. To account for sample correlations you can use duplicateCorrelation from the limma package.