[BioC] Desing matrix in limma

[James W. MacDonald]

Hi Dario,

On 4/13/2011 6:46 AM, Dario Beraldi wrote:
> 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?

It is correct, in the strict sense that you are doing what you think you 
are doing.

However, one could argue that it would be a better idea to fit the 
animals as a random effect rather than a fixed effect, using 
duplicateCorrelation() to account for the intra-animal covariance structure.

There are multiple examples of using duplicateCorrelation() in the Limma 
User's Guide.

Best,

Jim
>
> All the best
>
> Dario
>
>

-- 
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues