[BioC] limma question

[Gordon Smyth]

At 01:28 AM 3/04/2004, ivan.borozan at utoronto.ca wrote:
>hi there,
>
>I would like to estimate the effect on gene expression levels of two 
>factors Age
>and activity (each with 4 levels) using Limma.
>
>for my design matrix i have
>
>dataB<-data.frame(samples = arrays,Age=factor(FFAge),activity = 
>factor(FFfibrosis))
>
>design<-model.matrix(~activity*Age, data=dataB)
>
>(samples contain the names of my arrays)
>
>fit <- lm.series(MANormBetween$M,design)
>
>toptable(coef=2,num=10,genelist=gal,fit=fit,adjust="fdr",sort.by="P")
>
>I would like to know to which effect the quoted P.values in toptable()
>correspond to ?

Look at colnames(design). You've asked for the 2nd coefficient.

>Also i would like to know how to extract the P.value for the activity:Age 
>effect.

The activity:effect interaction is on 9 degrees of freedom and I assume 
that you understand that you need an F-statistic rather than a t-test 
statistic to test this composite hypothesis. Your fitted factorial model 
has 16 parameters of which the last 9 are interaction terms.

Using limma 1.5.2 or later, you can use

fit <- lmFit(MANormBetween at M, design)
cont.matrix <- rbind( matrix(0,7,9), diag(9) ) # pick out last 9 coefficients
fit <- contrasts.fit(fit, cont.matrix)
fit <- eBayes(fit)
F.stat <- FStat(fit)
P.value <- pf(F.stat, df1=attr(F.stat,"df1"), df2=attr(F.stat,"df2"), 
lower.tail=FALSE)

Gordon

>all the best.