Question: Creating contrasts for within group effects, interaction and overall effects in DESeq2
gravatar for January
15 months ago by
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.

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


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))
model matrix contrast deseq2 • 647 views
ADD COMMENTlink modified 13 months ago by Biostar ♦♦ 20 • written 15 months ago by January30
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1111 users visited in the last hour