DESeq2 multi-factor designs

[id0]

In DESeq2, you can set up a multi-factor design like:

design(dds) <- formula(~ celltype + treatment)

This would give you the differences attributable to the "treatment" while accounting for the "celltype".

However, what if you are interested in the differences in response to treatment in different cell types? Instead of looking for the effects of treatment regardless of cell type, you want the effects of treatments that are unique to a specific cell type. How would you get those?

[dpryan]

Either include an interaction (~celltype*treatment) or just make a "group" factor like:

group <- factor(sprintf("%s_%s", celltype, treatment))

and then use a model of ~group and use contrasts for the comparisons of interest. You'll get the same results with either method.

[Michael Love]

Devon's right with a small caveat. Because of log fold change shrinkage (betaPrior=TRUE is default), the results are not exactly the same (they would be the same with betaPrior=FALSE). The interaction model only shrinks the interaction term, while the group model shrinks all the group effects. My advice is: if you are interested in testing for cell type specificity, use the interaction model, because the interaction terms allow you to easily test this. If you are only interested in comparing the treatment effect in each cell type, use the group model, and the contrasts are easy to form this way.

[id0]

I didn't realize that the interaction and group models would give the same results (at least with betaPrior=FALSE). That is really interesting.

So if I have two cell types (A and B) and two treatments (Yes and No) and then I combine them for the group model, which contrast is equivalent to the interaction model? I am not sure how I could get the difference in response to treatment using the group model.

[Michael Love]

The contrast is (ctA_trY + ctB_trN) - (ctA_trN + ctB_trY)

This is from doing arithmetic on (ctA_trY - ctA_trN) - (ctB_trY - ctB_trN)

Here's an example in base R. Note the interation term (-2.11)

> y <- rnorm(8)
> ct <- c(0,0,1,1,0,0,1,1)
> tr <- c(0,0,0,0,1,1,1,1)
> lm(y ~ ct*tr)

Call:
lm(formula = y ~ ct * tr)

Coefficients:
(Intercept)           ct           tr        ct:tr
    -0.6356       1.1513       0.7684      -2.1120

Now with group. here, I add + 0 to the design, but you wouldn't do this in DESeq2 with betaPrior=TRUE because we fit a term for each group in addition to an intercept.

> g <- factor(c(1,1,2,2,3,3,4,4))
> lm(y ~ g + 0)

Call:
lm(formula = y ~ g + 0)

Coefficients:
     g1       g2       g3       g4
-0.6356   0.5157   0.1329  -0.8278

> beta <- coef(lm(y ~ g + 0))
> (beta[1] + beta[4]) - (beta[2] + beta[3])
       g1
-2.111991

[id0]

Thank you for that great explanation.

The only part I don't understand is this:

The contrast is (ctA_trY + ctB_trN) - (ctA_trN + ctB_trY)

How would I provide that to the DESeq2 results function contrasts parameter?

[Michael Love]

you would use the list style of specifying contrasts, see ?results for more details

results(dds, contrast=list(c("ctA_trY","ctB_trN"),c("ctA_trN","ctB_trY")))