[BioC] confusion using model.matrix with two-color arrays in limma

[Gordon K Smyth]

> Date: Mon, 1 Nov 2004 13:55:12 -0800 (PST)
> From: Dick Beyer <dbeyer at u.washington.edu>
> Subject: [BioC] confusion using model.matrix with two-color arrays in
> 	limma
> To: Bioconductor <bioconductor at stat.math.ethz.ch>
> Message-ID:
> 	<Pine.LNX.4.43.0411011355120.12459 at hymn02.u.washington.edu>
> Content-Type: TEXT/PLAIN; charset=ISO-8859-1; format=flowed
>
> I am struggling with trying to use model.matrix, which is straight-forward for
> single color arrays, on two-color dye-swap data. For two-color data, using modelMatrix is simple
> and straight-forward.
>
> I would like to know how to sort of combine these two calls, model.matrix and modelMatrix, so I
> can easily specify my design matrix (with modelMatrix), and easily specify complex contrasts (as
> with model.matrix).
>
> My particular experiment is a 2x2x2 design with factors, T,R,C each having levels 0/1.  I would
> like to specify the interaction T*R*C, without having to explicitly type it out using
> makeContrasts().
>
> How can I use model.matrix() so the dye-swaps are combined correctly (as is done easily with
> modelMatrix), or how can I specify makeContrasts() so I don't have to type out the T*R*C
> interaction?
>
> My targets are (where H signifies the case for parameter T=0, T signifies parameter T=1, etc):
>  	Cy3	Cy5
> 1	REF	HC
> 2	HC	REF
> 3	REF	HRC
> 4	HRC	REF
> 5	REF	TC
> 6	TC	REF
> 7	REF	TRC
> 8	TRC	REF
> 9	REF	T
> 10	T	REF
> 11	REF	TR
> 12	TR	REF
> 13	REF	H
> 14	H	REF
> 15	REF	HR
> 16	HR	REF
> My parameters are:
>  	T	R	C
> 1	0	0	1
> 2	0	0	1
> 3	0	1	1
> 4	0	1	1
> 5	1	0	1
> 6	1	0	1
> 7	1	1	1
> 8	1	1	1
> 9	1	0	0
> 10	1	0	0
> 11	1	1	0
> 12	1	1	0
> 13	0	0	0
> 14	0	0	0
> 15	0	1	0
> 16	0	1	0
>
> If I use the modelMatrix and makeContrast approach, I would use something like the following to
> get the T*R*C interaction:
> design <- modelMatrix(targets,ref="REF")
> fit <- lmFit(RG, design)
> contrast.matrix <- makeContrasts("H-HC-HR+HRC-T+TC+TR-TRC",levels=design)
> fit2 <- contrasts.fit(fit, contrasts.matrix)
>
> If I use the model.matrix approach (ignoring the necessary dye-swapping):
> dd <- data.frame(
> T=factor(c(0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0)),
> R=factor(c(0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1)),
> C=factor(c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0)))
> design.lm <- model.matrix(~  T*R*C, data=dd, contrasts = list(T="contr.sum", R="contr.sum",
> C="contr.sum"))
> fit <- lmFit(RG, design.lm)
> and have no need for the contrasts.fit() call.
>
> If I knew how to correctly specify a dye-flip parameter in the "dd" data.frame, I think I would be
> where I'd like to be.

Just multiply each row of design.lm by -1 whenever Cy5 is REF:

dyeswap <- 1-2*(targets$Cy5=="REF")
design.lm <- dyeswap * design.lm
fit <- lmFit(RG, design.lm)

BTW there's no compelling reason to incorporate dye-swaps into your experiment if you're using a
common reference.  And if REF was always Cy3, then you could use model.matrix() to create the
design matrix.

Gordon

> Any help or suggestions to try from anyone would be greatly appreciated.
> Thanks very much,
> Dick
> *******************************************************************************
> Richard P. Beyer, Ph.D.	University of Washington
> Tel.:(206) 616 7378	Env. & Occ. Health Sci. , Box 354695
> Fax: (206) 685 4696	4225 Roosevelt Way NE, # 100
>  			Seattle, WA 98105-6099
> http://depts.washington.edu/ceeh/ServiceCores/FC5/FC5.html