Hi,

After browsing through many of the posts on biostars, seqanswers etc, I still am a bit shaky on how best to handle multiple treatments with DESeq2.  I understand that the analysis takes the condition into account during the creation of the deseqDataSet object via the design argument, but how is that handled by the DESeq function?  

For instance, when using the plotMA function following this guide, which conditions are being plotted against the control by default?  I get a nice plot, but which contrast is being plotted?

Below is how I am going about constructing the deseqDataSet object":

countsTable <-data.matrix(df[0:2301,1:10])
samples<-names(df[1:10])
condition<-c(rep("ctrl",2),rep("A",2),rep("B",2),rep("C",2),rep("D",2))
pData = data.frame(cbind(samples, condition))
ddsfm <- DESeqDataSetFromMatrix(countData = countsTable, colData=pData, design=~condition)

dds<-DESeq(ddsfm)

conds<-unique(condition)

I had assumed that a way to index out the individual contrasts would be to use the following results functions:

aRes<-results(dds, contrast=c("condition", "A", "ctrl"))
bRes<-results(dds, contrast=c("condition", "B", "ctrl"))
cRes<-results(dds, contrast=c("condition", "C", "ctrl"))
dRes<-results(dds, contrast=c("condition", "D", "ctrl"))

This doesn't seem to be returning the contrasts that I was anticipating.  Is there a more idiomatic way to approach this process? 

> 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  
[9] base     

other attached packages:
 [1] plyr_1.8.1              Biobase_2.26.0          DESeq2_1.6.3           
 [4] RcppArmadillo_0.4.600.0 Rcpp_0.11.3             GenomicRanges_1.18.4   
 [7] GenomeInfoDb_1.2.4      IRanges_2.0.1           S4Vectors_0.4.0        
[10] BiocGenerics_0.12.1     ggplot2_1.0.0          

Better than passing 'dds' to plotMA() is to follow the example in the DESeq2 vignette (see section 1.4.1 MA-plot), where we pass the object returned by results() to plotMA(). This allows you to specify what comparison you want to plot:

vignette("DESeq2")

If you supply 'dds' to plotMA without any other arguments, then it is calling results() internally without any argument. To see what happens when you do this, we can check the man page for ?results:

"If results is run without specifying contrast or name, it will return the comparison of the last level of the last variable in the design formula over the first level of this variable."

DESeq() uses the design to build a generalized linear model with terms for the different factor levels. When you contrast D vs Ctrl, it is comparing the difference between the D samples and the Ctrl samples.

"This doesn't seem to be returning the contrasts that I was anticipating."

You can look at an individual row of the results table and compare the log2 fold change with a plot of the normalized counts for each group using the plotCounts() function (see vignette or ?plotCounts). Note that the log2FoldChange column is a posterior estimate (see our publication for more details). If you want to also see the raw log2 fold change, you can use:

results(dds, ..., addMLE=TRUE)