Creating contrasts for within group effects, interaction and overall effects in DESeq2

Question

[January]

There are two groups, A and B, with a number of individuals in each group. Each individual was treated twice (Test + Control). This is a hierarchical design with repeated measures. I am interested in the following contrasts:

• Between Test and Control in group A
• Between Test and Control in group B
• Interaction between group and treatment
• overall effect of treatment in both groups

Here is an example sample set:

t <- data.frame(group=rep(c("A", "B"), each=4),
                         treat=rep(c("C", "T"), 4),
                         pair=rep(paste0("P", 1:4), each=2))

Now, a naive approach would be something like this:

 t$group.treat <- paste0(t$group, "_", t$treat)
 d <- model.matrix(~ 0 + group.treat + pair, data=t)

And then define the approppriate contrasts (e.g. "(A_C-A_T)-(B_C-B_T)" for the interaction). However, this does not work, because DESeq2 throws an error:

m <- matrix(sample(1:1e6, 8 * 1000, replace=T), ncol=8)
foo.ds2 <- DESeqDataSetFromMatrix(countData=m, colData=t, design=d)
Error in checkFullRank(modelMatrix) : 
  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.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')

I have followed the recommendations in the vignette and elsewhere, and did this:

 t$pair2 <- rep(paste0("P", 1:2), each=2)
 d <- model.matrix(~ group + group:pair2 + group:treat, data=t)

This works and produces the following results names in the DESeq2 object:

> resultsNames(foo.ds2)
[1] "Intercept"      "groupB"         "groupA.pair2P2" "groupB.pair2P2" "groupA.treatT"  "groupB.treatT"

To get the within group effects and the interaction, I run the following:

res.withinA <- results(foo.ds2, name="groupA.treatT")
res.withinB <- results(foo.ds2, name="groupB.treatT")
res.interaction <- results(foo.ds2, contrast=list("groupA.treatT", "groupB.treatT"))

However, how can I define the contrast that shows the overall effect of the group. Would that be the correct approach:

res.AandB <- results(foo.ds2, contrast=c(0, 0, 0, 0, 1, 1))