DESeq2 multivariate analysis - retrieving certain stats

[Boel]

Dear All,

I am running a DESeq2 within-subject treatment response analysis. I have RNA seq data on 75 individuals, for most of these I have data before and after treatment (but not all). I also have a variable indicating successful treatment response across all individuals. I want to model gene expression as dependent on treatment and treatment response (res01) within each subject. counts is an R object with counts across all samples and genes.

The design matrix:

> head(A_design)
     sampleID subject treatment res01
A10a     A10a     A10         0     0
A10b     A10b     A10         1     1
A11a     A11a     A11         0     0
A11b     A11b     A11         1     0
A12a     A12a     A12         0     0
A12b     A12b     A12         1     1

The analysis is ran as follows:

dds <- DESeqDataSetFromMatrix(countData = counts[,A_design[,1]], colData = as.data.frame(A_design), design = ~ subject + treatment + res01)
Atreat.dds <- DESeq(dds)

Originally I though that I could extract the effect of res01 on gene expression as:

> res <- results(Atreat.dds,name='res01',pAdjustMethod='BH')
Error in results(Atreat.dds, name = "res01", pAdjustMethod = "BH") : 
  cannot find appropriate results in the DESeqDataSet.
possibly nbinomWaldTest or nbinomLRT has not yet been run.

At this point I realized that the results looked differently than I expected:

> resultsNames(Atreat.dds)
 [1] "Intercept"  "subjectA01" "subjectA02" "subjectA03" "subjectA04"
 [6] "subjectA05" "subjectA06" "subjectA07" "subjectA08" "subjectA09"
.....
[76] "subjectA92" "treatment0" "treatment1" "res010"     "res011"

I have two questions regarding this:
1. Why do I get two res01 results and two treatment results? This seems to be due to the within-subject design. I don't understand which model that has been tested here explicitly.
2. Should I run the following code to extract the effect that res01 has on gene expression given all other covariates?

Ares_res01=results(Atreat.dds, contrast=list("res010",'res011'))

Thanks in advance for any help,
Boel

My sessionInfo:

> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Scientific Linux release 6.7 (Carbon)

locale:
 [1] LC_CTYPE=sv_SE.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=sv_SE.UTF-8        LC_COLLATE=sv_SE.UTF-8    
 [5] LC_MONETARY=sv_SE.UTF-8    LC_MESSAGES=sv_SE.UTF-8   
 [7] LC_PAPER=sv_SE.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=sv_SE.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] DESeq2_1.8.1              RcppArmadillo_0.5.400.2.0
[3] Rcpp_0.12.0               GenomicRanges_1.20.5     
[5] GenomeInfoDb_1.4.2        IRanges_2.2.7            
[7] S4Vectors_0.6.3           BiocGenerics_0.14.0      

loaded via a namespace (and not attached):
 [1] RColorBrewer_1.1-2   futile.logger_1.4.1  plyr_1.8.3          
 [4] XVector_0.8.0        futile.options_1.0.0 tools_3.2.1         
 [7] rpart_4.1-10         digest_0.6.8         RSQLite_1.0.0       
[10] annotate_1.46.1      gtable_0.1.2         lattice_0.20-33     
[13] DBI_0.3.1            proto_0.3-10         gridExtra_2.0.0     
[16] genefilter_1.50.0    stringr_1.0.0        cluster_2.0.3       
[19] locfit_1.5-9.1       nnet_7.3-10          grid_3.2.1          
[22] Biobase_2.28.0       AnnotationDbi_1.30.1 XML_3.98-1.3        
[25] survival_2.38-3      BiocParallel_1.2.20  foreign_0.8-66      
[28] latticeExtra_0.6-26  Formula_1.2-1        geneplotter_1.46.0  
[31] ggplot2_1.0.1        reshape2_1.4.1       lambda.r_1.1.7      
[34] magrittr_1.5         scales_0.2.5         Hmisc_3.16-0        
[37] MASS_7.3-43          splines_3.2.1        xtable_1.7-4        
[40] colorspace_1.2-6     stringi_0.5-5        acepack_1.3-3.3     
[43] munsell_0.4.2