I am working with a dataset that has three different treatments (control, A, B). Each sample has only one of these. Each of control, A and B, has three replicates. I want to test for the effect of A on the one hand and the effect of B on the other.

Matrix headers:

ctrl1   trtB1   trtA1   ctrl2   trtB2   trtA2   ctrl3   trtB3  trtA3

Treatment factors:

Treat <- factor(substring(colnames(data.set),1,4))
Treat <- relevel(Treat, ref="ctrl")
Batch <- factor(substring(colnames(data.set),5,5))

Results:

> Treat
[1] ctrl trtB trtA ctrl trtB trtA ctrl trtB trtA

Levels: ctrl trtA trtB

> Batch
[1] 1 1 1 2 2 2 3 3 3
Levels: 1 2 3

Design matrix:

design <- model.matrix(~Batch+Treat)
rownames(design) <- colnames(y)

Results:

> design
      (Intercept) Batch2 Batch3 TreattrtA TreattrtB
ctrl1           1      0      0         0         0
trtB1           1      0      0         0         1
trtA1           1      0      0         1         0
ctrl2           1      1      0         0         0
trtB2           1      1      0         0         1
trtA2           1      1      0         1         0
ctrl3           1      0      1         0         0
trtB3           1      0      1         0         1
trtA3           1      0      1         1         0
attr(,"assign")
[1] 0 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$Batch
[1] "contr.treatment"

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

Now, my question is this:

If I run...

lrt.trtA <- glmLRT(fit, coef=4)
lrt.trtB <- glmLRT(fit, coef=5)

...will this produce the genes that are differentially expressed between trtA and control (for coef=4) and between trtB and control (for coef=5)? Will a positive fold change indicate an upregulation in e. g. trtA compared with ctrl, since ctrl was set as ref with relevel()?

The reason I am a bit unsure is because this is the first time that I am using a treatment factor that has more than two levels.

Yep, looks good to me. Dropping each coef will test for DE between the corresponding treatment and the control.

Edit: Hold on, why do samples treatA1-3 have values of 1 in TreattrtB, and samples treatB1-3 have values of 1 in TreattrtA? It should be the other way around with the model.matrix call you've used.

> Treat <- factor(rep(c("ctrl", "trtA", "trtB"), 3))
> Treat <- relevel(Treat, ref="ctrl")
> Batch <- factor(rep(1:3, each=3))
> design <- model.matrix(~Batch+Treat)
> design
  (Intercept) Batch2 Batch3 TreattrtA TreattrtB
1           1      0      0         0         0
2           1      0      0         1         0
3           1      0      0         0         1
4           1      1      0         0         0
5           1      1      0         1         0
6           1      1      0         0         1
7           1      0      1         0         0
8           1      0      1         1         0
9           1      0      1         0         1