Question: Creating contrasts for within group effects, interaction and overall effects in DESeq2
0
January30 wrote:

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.

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)
 "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))
``````
model matrix contrast deseq2 • 647 views