[BioC] Desing matrix in limma

[Dario Beraldi]

Hello,

I'm using limma to analyze affymetrix arrays from a time course  
experiment (time points 0, 2, 7, 24). We are interested in comparing 0  
vs 2, 0 vs 7, 0 vs 24.

The experimental design looks like this:

anim_id<- as.factor(rep(c('L1', 'L2', 'L3'), each= 4))
time_point<- rep(c(0, 2, 7, 24), 3)
slide_id<- paste(anim_id, time_point, sep= '.')
(targets<- data.frame(slide_id, anim_id, time_point))

#   slide_id anim_id time_point
#1      L1.0      L1          0
#2      L1.2      L1          2
#3      L1.7      L1          7
#4     L1.24      L1         24
#5      L2.0      L2          0
#6      L2.2      L2          2
#7      L2.7      L2          7
#8     L2.24      L2         24
#9      L3.0      L3          0
#10     L3.2      L3          2
#11     L3.7      L3          7
#12    L3.24      L3         24

I wanted to analyze all the arrays at once to fully take advantage of  
the eBayes function. Also I wanted to capitalize on the fact that the  
same animal is sampled at each time point (0, 2, 7, 24 hours).

So I set up the design matrix like this:

design<- model.matrix(~0+anim_id)
rownames(design)<- paste(targets$anim_id, targets$time_point, sep= '.')
design<- cbind(design, tp.2vs0= c(0,1,0,0, 0,1,0,0, 0,1,0,0))  ## 0 vs 2
design<- cbind(design, tp.7vs0= c(0,0,1,0, 0,0,1,0, 0,0,1,0))  ## 0 vs 7
design<- cbind(design, tp.24vs0= c(0,0,0,1, 0,0,0,1, 0,0,0,1)) ## 0 vs 24
design
#       anim_idL1 anim_idL2 anim_idL3 tp.2vs0 tp.7vs0 tp.24vs0
# L1.0          1         0         0       0       0        0
# L1.2          1         0         0       1       0        0
# L1.7          1         0         0       0       1        0
# L1.24         1         0         0       0       0        1
# L2.0          0         1         0       0       0        0
# L2.2          0         1         0       1       0        0
# L2.7          0         1         0       0       1        0
# L2.24         0         1         0       0       0        1
# L3.0          0         0         1       0       0        0
# L3.2          0         0         1       1       0        0
# L3.7          0         0         1       0       1        0
# L3.24         0         0         1       0       0        1

Now I can procede 'as usual':
fit <- lmFit(limma_affy, design)
fit <- eBayes(fit)

topTable(fit, coef= 4) ## 0 vs 2
topTable(fit, coef= 5) ## 0 vs 7
topTable(fit, coef= 6) ## 0 vs 24

I verified that the logFC I get from topTable is the same as the logFC  
calculated manually, which tells me that the design is correct.

However, I would like to make sure that this procedure is correct,  
Would anyone be so kind to confirm it?

All the best

Dario


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.