Dear all,

I am working with count data from a sequencing experiment and need some help with the design. 

I have 4 treatment groups which can be accounted by a factorial design. Defining 2 variables, A and B, the control group would be 1,1; the A group would be 2,1; the B groups would be 1,2; and the AB group would be 2,2 in the design matrix.

For all the groups I have a measurement before (T0) and after the treatment (T1). Every group had 10 individuals, so I have in total 80 measurement, 10 for each Treatment*Time.

I have used de design ~A*B*Time, but can not retrieve the contrast I look for.

I want to contrast for different times within a treatment, and for a given time (T0 or T1) the contrast between treatments. 

For example, contrast="A1.B1.TimeT0" vs "A1.B1.TimeT1"

or contrast="A2.B1.TimeT1" vs "A1.B1.TimeT1"

but when I try to ask R for that it tells me that I can only contrast the parameters shown in resultsName where the combination of interaction are not to be found.

 resultsNames(diagdds0)
[1] "Intercept"     "A_2_vs_1"      "B_2_vs_1"      "Time_T1_vs_T0"
[5] "A2.B2"         "A2.TimeT1"     "B2.TimeT1"     "O2.P2.TimeT1"

 I think I could get the contrast that I want if I the variables them selves were parameters, i.e "A" and "B" "Time" and I could use contrast=list(c(), c()),  and concatenate different levels of the variables. But there not there! were they absorbed by the intercept? what does that mean? or is it that I can only access to the information of the the second level because the first is considered as the reference?  

When I try to use the design ~A+B+Time+A:B+A:B:Time it gives the error, some of the variables are lineal combination of others. (This would also happen if I do the model ~A:B. I don´t get why they would be a lineal combination if there are different combination of 1 and 2 for each treatment and have different individuals having those different combination).

I have created a new variable where every treatment is a factor and used the design ~Treatment*Time, and it gave me the combinations I want to contrast, but I think that the other design should account better for the variance and wish to use it instead. Hope some one can give me a hand.

Thank you very much in advance, cheers

Isabel.

     
Design matrix
SampleID   A B Time
13.1 2 2   T0
25.2 2 2   T1
34.1 1 2   T0
38.1 1 2   T0
2.2  2 1   T1
1.1  1 2   T0
10.2 1 1   T1
8.1  2 1   T0
26.2 1 1   T1
2.1  2 1   T0
15.1 1 2   T0
9.1  2 2   T0
39.1 1 2   T0
29.2 2 2   T1
1.2  1 2   T1
19.1 2 1   T0
34.2 1 2   T1
14.2 2 1   T1
18.2 1 1   T1
35.2 2 2   T1
14.1 2 1   T0
11.1 2 1   T0
12.2 2 2   T1
18.1 1 1   T0
28.1 1 2   T0
13.2 2 2   T1
20.1 1 1   T0
9.2  2 2   T1
6.2  1 2   T1
15.2 1 2   T1
27.1 1 2   T0
36.1 1 1   T0
8.2  2 1   T1
10.1 1 1   T0
36.2 1 1   T1
37.2 2 1   T1
39.2 1 2   T1
31.1 2 1   T0
30.1 2 1   T0
17.2 1 2   T1
35.1 2 2   T0
31.2 2 1   T1
32.2 2 1   T1
6.1  1 2   T0
27.2 1 2   T1
32.1 2 1   T0
33.1 1 1   T0
37.1 2 1   T0
19.2 2 1   T1
17.1 1 2   T0
3.2  1 1   T1
11.2 2 1   T1
26.1 1 1   T0
29.1 2 2   T0
30.2 2 1   T1
12.1 2 2   T0
20.2 1 1   T1
3.1  1 1   T0
23.1 1 1   T0
7.1  1 2   T0
16.2 2 2   T1
41.1 2 2   T0
16.1 2 2   T0
7.2  1 2   T1
4.2  1 1   T1
41.2 2 2   T1
38.2 1 2   T1
21.1 2 2   T0
22.1 1 2   T0
5.2  2 2   T1
25.1 2 2   T0
22.2 1 2   T1
4.1  1 1   T0
40.2 1 1   T1
40.1 1 1   T0
21.2 2 2   T1
23.2 1 1   T1
5.1  2 2   T0
24.2 2 1   T1
24.1 2 1   T0
28.2 1 2   T1


> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq2_1.6.3              RcppArmadillo_0.4.600.4.0 Rcpp_0.11.4               GenomicRanges_1.18.4     
 [5] GenomeInfoDb_1.2.4        IRanges_2.0.1             S4Vectors_0.4.0           BiocGenerics_0.12.1      
 [9] plyr_1.8.1                phyloseq_1.10.0          

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3      ade4_1.6-2           annotate_1.44.0      AnnotationDbi_1.28.1 ape_3.2             
 [6] base64enc_0.1-2      BatchJobs_1.5        BBmisc_1.9           Biobase_2.26.0       BiocParallel_1.0.3  
[11] biom_0.3.12          Biostrings_2.34.1    bitops_1.0-6         brew_1.0-6           checkmate_1.5.1     
[16] cluster_2.0.1        codetools_0.2-10     colorspace_1.2-4     data.table_1.9.2     DBI_0.3.1           
[21] digest_0.6.8         evaluate_0.5.5       fail_1.2             foreach_1.4.2        foreign_0.8-62      
[26] formatR_1.0          Formula_1.2-0        genefilter_1.48.1    geneplotter_1.44.0   ggplot2_1.0.0       
[31] grid_3.1.2           gtable_0.1.2         Hmisc_3.15-0         igraph_0.7.1         iterators_1.0.7     
[36] knitr_1.9            lattice_0.20-29      latticeExtra_0.6-26  locfit_1.5-9.1       MASS_7.3-39         
[41] Matrix_1.1-5         mgcv_1.8-4           multtest_2.22.0      munsell_0.4.2        nlme_3.1-119        
[46] nnet_7.3-9           permute_0.8-3        plotly_0.5.24        proto_0.3-10         RColorBrewer_1.1-2  
[51] RCurl_1.95-4.5       reshape2_1.4.1       RJSONIO_1.3-0        rpart_4.1-9          RSQLite_1.0.0       
[56] scales_0.2.4         sendmailR_1.2-1      splines_3.1.2        stringr_0.6.2        survival_2.37-7     
[61] tools_3.1.2          vegan_2.2-1          XML_3.98-1.1         xtable_1.7-4         XVector_0.6.0       
[66] zlibbioc_1.12.0  

hi Isabel,

Recently I answered another question on the forum about designs with 2nd order interactions

DESeq2: with multiple factors and interaction terms won't show all effects

Basically, the DESeq2 code should have prevented you from running a model with a beta prior and a second order interaction (because we haven't implemented how the prior should work in this case), and suggested you use either 1) a single grouping variable with all unique combinations of the variables or 2) betaPrior=FALSE. In v1.6 this was not being caught, but I have fixed this in v1.7 so it gives a helpful message.

My recommendation is to combine all 3 variables into one, and then contrasts will be easy:

dds$group <- factor(paste0(dds$A, dds$B, dds$Time))
design(dds) <- ~ group
dds <- DESeq(dds)
results(dds, contrast=c("group","11T1","11T0")) # for T1 vs T0