Using DEseq2, the lfcShrink function with the apeglm method requires a single coef value to be specified (i.e. it does not allow a contrast to be specified). Further explanation in the manual.

Therefore, I am having trouble creating a design matrix that represents my comparison of interest. Here is my example code:

cells <- c(rep("A", 6), rep("B", 6))
treat <- c(rep("C", 3), rep("T", 3), rep("C", 3), rep("T", 3))
replicate <- c(rep(c(1:3), 4))
group <- rep(c("W", "X", "Y", "Z"), each = 3)

colData <- data.frame(sample_id = paste0(cells, "_", treat, "_", replicate), cells, treat, replicate, group)

colData$cells <- factor(colData$cells, levels = c("A", "B"))
colData$treat <- factor(colData$treat, levels = c("C", "T"))
colData$group <- factor(colData$group, levels = c("W", "X", "Y", "Z"))

design1 <- model.matrix(~ cells + treat + cells:treat, data = colData)
design2 <- model.matrix(~ group, data = colData)

And the relevant output:

> colData
   sample_id cells treat replicate group
1      A_C_1     A     C         1     W
2      A_C_2     A     C         2     W
3      A_C_3     A     C         3     W
4      A_T_1     A     T         1     X
5      A_T_2     A     T         2     X
6      A_T_3     A     T         3     X
7      B_C_1     B     C         1     Y
8      B_C_2     B     C         2     Y
9      B_C_3     B     C         3     Y
10     B_T_1     B     T         1     Z
11     B_T_2     B     T         2     Z
12     B_T_3     B     T         3     Z

> design1
   (Intercept) cellsB treatT cellsB:treatT
1            1      0      0             0
2            1      0      0             0
3            1      0      0             0
4            1      0      1             0
5            1      0      1             0
6            1      0      1             0
7            1      1      0             0
8            1      1      0             0
9            1      1      0             0
10           1      1      1             1
11           1      1      1             1
12           1      1      1             1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$cells
[1] "contr.treatment"

attr(,"contrasts")$treat
[1] "contr.treatment"

> design2
   (Intercept) groupX groupY groupZ
1            1      0      0      0
2            1      0      0      0
3            1      0      0      0
4            1      1      0      0
5            1      1      0      0
6            1      1      0      0
7            1      0      1      0
8            1      0      1      0
9            1      0      1      0
10           1      0      0      1
11           1      0      0      1
12           1      0      0      1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

I want to compare T vs C treatments within only the A cells. This seems easy with design2, by specifying coef = 2 (i.e. groupX). In that case, the baseline levels are A cells and C treatment. However, I cannot understand how to use an additive (interaction) formula in the model.matrix function (i.e. design1), where explicitly naming multiple factors (e.g. cells, treat) will produce a design matrix to represent this comparison?

Thanks!

Thanks for providing the colData and model.matrix tables, that made it easy to see what you are looking for.

I think this design gives you what you want:

> model.matrix(~cells + cells:treat)
   (Intercept) cellsB cellsA:treatT cellsB:treatT
1            1      0             0             0
2            1      0             0             0
3            1      0             0             0
4            1      0             1             0
5            1      0             1             0
6            1      0             1             0
7            1      1             0             0
8            1      1             0             0
9            1      1             0             0
10           1      1             0             1
11           1      1             0             1
12           1      1             0             1