Is this design matrix correct (treament factor with three levels)?
1
0
Entering edit mode
Ekarl2 ▴ 80
@ekarl2-7202
Last seen 6.5 years ago
Sweden

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.

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.

edger limma design matrix • 1.5k views
3
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 3 hours ago
The city by the bay

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

0
Entering edit mode

Looks like it was a result of me changing the real treatment names to A and B labels too quickly (real name of treatment A is after real name of treatment B alphabetically), so they got switched. Let me fix that right away.

0
Entering edit mode

I have fixed the display problem now; trtB was in column 2, 5 and 8, whereas trtA was in column 3, 6, 9 in the original matrix file. So that order was not alphabetical. How does it look?

1
Entering edit mode

Fine. And yes, a positive log-fold change for either coef indicates upregulation over the control.