Hi!

This question is regarding the creation of a contrast matrix in limma, edgeR or DESeq2, while analyzing RNA-Seq data. The question is: "How to create contrasts which represent comparisons between two factors which are not the reference level (intercept) using a model which contains an interaction term?"

(I do believe the interaction term is the important component here, since the interaction term will lead to a more complicated contrast matrix.)

Let's have a look at an example:

# Create data frame with information.
groups <- data.frame(group = c(rep(1,6), rep(2,6), rep(3,6), rep(4,6)),
                     subject = c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8),
                     treat = c(rep("a",6),rep("b",6),rep("a",6),rep("b",6)),
                     time = c(rep(c(1,2,3), 4), rep(c(1,4,5), 4)),
                     outcome = c(3,2,4,5,3,4,9,8,7,8,8,10,2,3,1,1,3,2,9,9,7,6,8,6))

# Setting levels for fixed effects.
groups$treat <- factor(groups$treat)
groups$time <- factor(groups$time)
groups$group <- factor(groups$group)
groups$subject <- factor(groups$subject)

# Creating design formula.
design <- model.matrix(outcome ~ treat * time, data = groups)
design
#>    (Intercept) treatb time2 time3 time4 time5 treatb:time2 treatb:time3
#> 1            1      0     0     0     0     0            0            0
#> 2            1      0     1     0     0     0            0            0
#> 3            1      0     0     1     0     0            0            0
#> 4            1      0     0     0     0     0            0            0
#> 5            1      0     1     0     0     0            0            0
#> 6            1      0     0     1     0     0            0            0
#> 7            1      1     0     0     0     0            0            0
#> 8            1      1     1     0     0     0            1            0
#> 9            1      1     0     1     0     0            0            1
#> 10           1      1     0     0     0     0            0            0
#> 11           1      1     1     0     0     0            1            0
#> 12           1      1     0     1     0     0            0            1
#> 13           1      0     0     0     0     0            0            0
#> 14           1      0     0     0     1     0            0            0
#> 15           1      0     0     0     0     1            0            0
#> 16           1      0     0     0     0     0            0            0
#> 17           1      0     0     0     1     0            0            0
#> 18           1      0     0     0     0     1            0            0
#> 19           1      1     0     0     0     0            0            0
#> 20           1      1     0     0     1     0            0            0
#> 21           1      1     0     0     0     1            0            0
#> 22           1      1     0     0     0     0            0            0
#> 23           1      1     0     0     1     0            0            0
#> 24           1      1     0     0     0     1            0            0
#>    treatb:time4 treatb:time5
#> 1             0            0
#> 2             0            0
#> 3             0            0
#> 4             0            0
#> 5             0            0
#> 6             0            0
#> 7             0            0
#> 8             0            0
#> 9             0            0
#> 10            0            0
#> 11            0            0
#> 12            0            0
#> 13            0            0
#> 14            0            0
#> 15            0            0
#> 16            0            0
#> 17            0            0
#> 18            0            0
#> 19            0            0
#> 20            1            0
#> 21            0            1
#> 22            0            0
#> 23            1            0
#> 24            0            1
#> attr(,"assign")
#>  [1] 0 1 2 2 2 2 3 3 3 3
#> attr(,"contrasts")
#> attr(,"contrasts")$treat
#> [1] "contr.treatment"
#> 
#> attr(,"contrasts")$time
#> [1] "contr.treatment"

^{Created on 2019-07-08 by the reprex package (v0.3.0)}

This example shows a table containing the factor levels for the design matrix. After that, you can see the design formula and the design matrix printed to the console.

This dummy dataset contains an outcome (Y), which would be not the case for an RNA-Seq experiment, but I used this example for testing where I wanted an outcome (Y) - so bear with me here.

The intercept represents our reference level, which is "treat a" and "time 1". "treat b" (column 2 in the design matrix) represents the difference of "treat b / time 1" compared to "treat a / time 1". "time 2, 3, 4, 5" represent the respective time differences of "treat a".

Columns 7 to 10 in the design matrix represent the additional effects in "treat b".

Goal is to create a contrast which represents the comparison of "time 5" compared to "time 1" of "treat b". Here, "treat b / time 1" should be our reference level.

I do believe this should be the effect of "treat a / time 5" plus the additional effect captured within the coefficient "treatb:time5", so my contrast would look like this:

makeContrast(contrast = time5 + treatb:time5, ...)

This results in a contrast matrix which does not equals zero, which for me looks wrong.

I did find explanations about expressing "differences of differences" using a contrast matrix (e.g. here), but I do not want to express differences of differences. I do want the logFC and p-values for the comparison of "treat b / time 5" compared to "treat b / time 1".

Thanks for your help!

So you are just asking the same question you already asked here, to which I said to combine the treatment and time. Evidently that didn't make any sense since you have just asked the same question again using a different post name (hint: don't do that - it's considered rude to just keep posting the same question in different posts).

Since words evidently didn't make sense, let's try an example.

> trttime <- paste(groups$treat, groups$time, sep = "_")
> design <- model.matrix(~0+trttime)
> colnames(design) <- gsub("trttime", "", colnames(design))
> design
   a_1 a_2 a_3 a_4 a_5 b_1 b_2 b_3 b_4 b_5
1    1   0   0   0   0   0   0   0   0   0
2    0   1   0   0   0   0   0   0   0   0
3    0   0   1   0   0   0   0   0   0   0
4    1   0   0   0   0   0   0   0   0   0
5    0   1   0   0   0   0   0   0   0   0
6    0   0   1   0   0   0   0   0   0   0
7    0   0   0   0   0   1   0   0   0   0
8    0   0   0   0   0   0   1   0   0   0
9    0   0   0   0   0   0   0   1   0   0
10   0   0   0   0   0   1   0   0   0   0
11   0   0   0   0   0   0   1   0   0   0
12   0   0   0   0   0   0   0   1   0   0
13   1   0   0   0   0   0   0   0   0   0
14   0   0   0   1   0   0   0   0   0   0
15   0   0   0   0   1   0   0   0   0   0
16   1   0   0   0   0   0   0   0   0   0
17   0   0   0   1   0   0   0   0   0   0
18   0   0   0   0   1   0   0   0   0   0
19   0   0   0   0   0   1   0   0   0   0
20   0   0   0   0   0   0   0   0   1   0
21   0   0   0   0   0   0   0   0   0   1
22   0   0   0   0   0   1   0   0   0   0
23   0   0   0   0   0   0   0   0   1   0
24   0   0   0   0   0   0   0   0   0   1
attr(,"assign")
 [1] 1 1 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$trttime
[1] "contr.treatment"

Now you can compare any group to any other group, without having to interpret what the different interaction terms mean. This obviously assumes that you don't care that A) you apparently have repeated measures from different subjects, and B) there is this outcome that you included but ignored as well.