Hi All,

I'm trying to put together a proper contrast matrix for DESeq2, and failing miserably. 

My experimental design involves using RNA-seq to look at two regions of the brain (Region1 and Region2) at two developmental stages (E and P). To get enough RNA I have pooled RNA from genetically "identical" mice to generate my replicates, so each sample contains 4-6 mice. However, for N=1 and N=2 I used the same animals to isolate region1 and region 2, while for N=3 I didn't have enough RNA, and had to use 1-2 "extra" mice for Region1 (but not Region2).

I am interested in the pairwise comparisons:

• Region1 vs Region2 at E
• Region1 vs Region2 at P
• E vs P in Region1
• E vs P in Region2

as well as the more general comparisons of:

• Region1 vs Region2  (irrespective of stage)
• E vs P (irrespective of region)

I've tried using design = ~ Animals + Stage + Region and design = ~ Animals + Stage + Region + Stage:Region

which gives me the error:

"Error in DESeqDataSet(se, design = design, ignoreRank) :   the model matrix is not full rank, so the model cannot be fit as specified.  one or more variables or interaction terms in the design formula are linear combinations of the others and must be removed

What am I doing wrong? And how do I extract the results for the comparisons I am interested in (and doing the combined comparisons make sense, right? Or should I just be looking at the overlap between the two paired comparisons?)? Or is it just impossible to construct a proper contrast matrix when I have only one replicate of the "Animals" factor that is == to 7 and 8?

Thanks in advance,
Darya

R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)

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

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

other attached packages:
 [1] ggplot2_1.0.0             DESeq2_1.6.1              RcppArmadillo_0.4.450.1.0
 [4] Rcpp_0.11.3               GenomicRanges_1.18.1      ReportingTools_2.6.0     
 [7] AnnotationDbi_1.28.0      GenomeInfoDb_1.2.0        IRanges_2.0.0            
[10] S4Vectors_0.4.0           Biobase_2.26.0            BiocGenerics_0.12.0      
[13] RSQLite_1.0.0             DBI_0.3.1                 knitr_1.7                
[16] edgeR_3.8.2               limma_3.22.1              biomaRt_2.22.0           

loaded via a namespace (and not attached):
 [1] acepack_1.3-3.3          annotate_1.44.0          AnnotationForge_1.8.1   
 [4] base64enc_0.1-2          BatchJobs_1.4            BBmisc_1.7              
 [7] BiocParallel_1.0.0       Biostrings_2.34.0        biovizBase_1.14.0       
[10] bitops_1.0-6             brew_1.0-6               BSgenome_1.34.0         
[13] Category_2.32.0          checkmate_1.5.0          cluster_1.15.3          
[16] codetools_0.2-9          colorspace_1.2-4         dichromat_2.0-0         
[19] digest_0.6.4             evaluate_0.5.5           fail_1.2                
[22] foreach_1.4.2            foreign_0.8-61           formatR_1.0             
[25] Formula_1.1-2            genefilter_1.48.1        geneplotter_1.44.0      
[28] GenomicAlignments_1.2.0  GenomicFeatures_1.18.1   GGally_0.4.8            
[31] ggbio_1.14.0             GO.db_3.0.0              GOstats_2.32.0          
[34] graph_1.44.0             grid_3.1.1               gridExtra_0.9.1         
[37] GSEABase_1.28.0          gtable_0.1.2             Hmisc_3.14-5            
[40] hwriter_1.3.2            iterators_1.0.7          labeling_0.3            
[43] lattice_0.20-29          latticeExtra_0.6-26      locfit_1.5-9.1          
[46] MASS_7.3-35              Matrix_1.1-4             munsell_0.4.2           
[49] nnet_7.3-8               OrganismDbi_1.8.0        PFAM.db_3.0.0           
[52] plyr_1.8.1               proto_0.3-10             R.methodsS3_1.6.1       
[55] R.oo_1.18.0              R.utils_1.34.0           RBGL_1.42.0             
[58] RColorBrewer_1.0-5       RCurl_1.95-4.3           reshape_0.8.5           
[61] reshape2_1.4             rpart_4.1-8              Rsamtools_1.18.1        
[64] rtracklayer_1.26.1       scales_0.2.4             sendmailR_1.2-1         
[67] splines_3.1.1            stringr_0.6.2            survival_2.37-7         
[70] tools_3.1.1              VariantAnnotation_1.12.2 XML_3.98-1.1            
[73] xtable_1.7-4             XVector_0.6.0            zlibbioc_1.12.0    

hi Darya,

You can also use the trick from the edgeR user guide, which I recently posted: DESeq2 paired multifactor test

You would recode the animal variable to:

dds$animal.nested <- factor(c(1,2,3,1,2,4,1,2,3,1,2,4))

and use a design of:

~ stage + stage:animal + stage:region

You can answer then, which show region 2 vs 1 differences in stage E, in stage P, as well as where there are differences in the region 2 vs 1 effect between E and P. See the examples at the linked post.

The E vs P comparisons are problematic if you think there is a strong animal effect, because there are no matched samples across the stages. For example, suppose animals 1,2,3 are very different than animals 4,5,6 which has nothing to do with stage. It makes it problematic to compare E vs P in region 1 if this is the case. You can make the comparison, but you would need to assume that the animal effect is negligible, whereas the comparisons I linked to above actually control for animal differences.