Is this design matrix correct (treament factor with three levels)?
1
0
Entering edit mode
Ekarl2 ▴ 80
@ekarl2-7202
Last seen 6.4 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.

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.

edger limma design matrix • 1.5k views
ADD COMMENT
3
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 6 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
ADD COMMENT
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.

ADD REPLY
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?

ADD REPLY
1
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 374 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6