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.

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,

You do not have the same experimental design as in the previous link, so the design recommendation doesn't work for you.

The experimental design in the link had treatment and control for each culture. You do not have treatment and control for each calf (the effect you want to control for). Calves are either in treatment group, or in control group.

I don't think you can control for the calf effect using fixed effects in DESeq2. The calf terms will be collinear with the treatment term. 

You can use limma instead, which has an option for informing about correlated samples. The function to look up would be limma's duplicateCorrelation.

Or you can either continue with DESeq2, and not control for calf effect, in which case, I would recommend the following design for all unique levels of day and treatment:

dds$group <- factor(paste0(dds$Day, dds$Treatment))
design(dds) <- ~ group
dds <- DESeq(dds)

Then you can contrast any of the groups you like using simple contrasts.

The average treatment effect is:

results(dds, contrast=list(c("TreatmentTreatDay36","TreatmentTreatDay54"),
c("TreatmentControlDay36","TreatmentControlDay54")),
listValues=c(.5, -.5))

Hi Mike,

I forgot to mention that the factor "day"  means the calves were sampled on day 36 and day 54. And the data I'm analyzing are NGS 16s RNA illumina data.

Basically, we spited calves  into two pens (groups) first. After a couples of weeks of adaption, all calves were samples once (d36) . Then one group of calves went into an abrupt weaning program (treatment) and another group of calves had a gradual weaning program (control). A couple of weeks after all calves were completely weaned, we samples them again (day54).   What we are interested in are 1) effect of weaning (treatment effect), 2) effect of day ( pre weaning vs post weaning), 3) interaction between treatment and day.

I really like to use DESeq2 to  analyze the data and don't want to explore Limma for now.  However, the variation of gut microbiota among individual calves is rather significant.  Please suggest a model ,if possible, in which "calf" can be taken into account.

Shucong

Hello Mike,

The post I reference is actually this one.

C: DESeq2 paired and multi factor comparison

Sorry for the mistake and the confusion it caused.  I have revised my original post. Please check it again and help me out.

Shucong