Factorial Design with DESeq2; contrast problem.
1
0
Entering edit mode
sifujio • 0
@sifujio-7453
Last seen 6.5 years ago
Chile

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  
deseq2 multiple factor design complex-design • 2.1k views
ADD COMMENT
1
Entering edit mode
@mikelove
Last seen 18 hours ago
United States

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
ADD COMMENT

Login before adding your answer.

Traffic: 432 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