[BioC] Affymetrix Limma design/contrast matrix

[Gordon Smyth]

At 04:53 AM 23/06/2004, Matthew  Hannah wrote:
>I know this has been asked several times for various designs, and I
>have searched and read the user guide but I'm getting nowhere fast.
>I would be very grateful if someone could help me out with what is
>probably a simple request to someone familar with lm and Limma.
>
>I was following
>8.4 Estrogen Data: A 2x2 Factorial Experiment with Affymetrix Arrays
>but have got a bit confused - especially if
> > cont.matrix <- cbind(E10=c(0,0,1,0),E48=c(0,0,0,1))
>is not a typo and should read
> > cont.matrix <- cbind(E10=c(0,1,0,0),E48=c(0,0,0,1))

No it is not a typo.

>Anyway rather than say more than I'm statistically inept, I would
>appreciate some help on an appropriate design and contrast matrix
>for the list below.
>
>        Exp Genotype Treatment
>MUTA.1   1      MUT         A
>MUTA.2   2      MUT         A
>MUTA.3   3      MUT         A
>MUTA.4   4      MUT         A
>MUTN.1   1      MUT         N
>MUTN.2   2      MUT         N
>MUTN.3   3      MUT         N
>MUTN.4   4      MUT         N
>ConA.1   1      Con         A
>ConA.2   2      Con         A
>ConA.3   3      Con         A
>ConA.4   4      Con         A
>ConN.1   1      Con         N
>ConN.2   2      Con         N
>ConN.3   3      Con         N
>ConN.4   4      Con         N
>
>I already have it as pData (is there an easy way
>to adapt this?). I tried this design (is it correct?) but also want
>it with the experiment included.
>
> >treatments <- factor(c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4),
>labels=c("MUTA","MUTN","ConA","ConN"))
> > contrasts(treatments) <- cbind(Treat=c(1,0,1,0),MUT=c(1,1,0,0),
>Con=c(0,0,1,1))
> >design <- model.matrix(~treatments)
>
>Then I got very confused with the contrasts - in the example they only
>look at the estrogen effect, what if you want to make the same contrasts
>as in the design (eg: also include time in the estrogen example) do you
>need another fit or do you just use the first one?

No you only need one fit.

>Basically I want to compare MUTA vs ConA, MUTN vs ConN, A vs N.

Perhaps the easiest for you is to use the makeContrasts() function:

treatments <- factor(c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4))
design <- model.matrix(~ 0+treatments)
colnames(design) <- c("MUTA","MUTN","ConA","ConN")
fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(MUTA-ConA, MUTN-ConN, 
AvsN=(MUTA+ConA-MUTN-ConN)/2, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

>Getting slightly more complicated the data is paired (eg: MUTA.1 with
>MUTN.1) and was wondering if this pairwise nature could be taken into
>account and compare the MUTA-MUTN changes vs ConA-ConN changes? I ask
>this as I've found that the changes may be more reproducible than the
>absolute values.

Now you are asking something which is a methodological research question, 
and you really should consider taking on a statistician as a full 
collaborator.
Gordon

>Thanks in advance.
>Matt