Question: DESeq2 nested design, multifactor test
3
0
Entering edit mode
Shucong Li ▴ 60
@shucong-li-6266
Last seen 8.7 years ago
Canada

 

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

deseq2 experimental design • 2.1k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 1 day ago
United States

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))

 

ADD COMMENT
0
Entering edit mode
Shucong Li ▴ 60
@shucong-li-6266
Last seen 8.7 years ago
Canada

 

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

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode
Shucong Li ▴ 60
@shucong-li-6266
Last seen 8.7 years ago
Canada

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

ADD COMMENT

Login before adding your answer.

Traffic: 645 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6