[BioC] limma and nested factors

[Arne.Muller at aventis.com]

Hello,

I'm following the limma users' guide to analyze some affy data. The design in
my analysis us a bit more complex than in the example of the manual. I don't
know if I'm doing it right.

I've two factors 'batch' and 'time' listed below:

> batch
 [1] NEW NEW NEW OLD OLD OLD OLD PRG PRG PRG NEW NEW NEW OLD OLD OLD OLD PRG
PRG
[20] PRG
Levels: NEW OLD PRG
> 

> time
 [1] 04h 04h 04h 04h 04h 04h 04h 04h 04h 04h 24h 24h 24h 24h 24h 24h 24h 24h
24h
[20] 24h
Levels: 04h 24h

You see that BATCH/TIME combinations have 3 replicates (except for OLD/04h or
24h which has 4).

I'd like to know how many gene change in expression over time but within each
batch. Therefore I'm constructing the following model matrix (time nested
within batch, i.e. 04h and 24h at each level of batch):

design.time <- model.matrix( ~ -1 + time %in% batch)
colnames(design.time) <- c("time04h.batchNEW",
                           "time24h.batchNEW",
                           "time04h.batchOLD",
                           "time24h.batchOLD",
                           "time04h.batchPRG",
                           "time24h.batchPRG")

> design.time
   time04h.batchNEW time24h.batchNEW time04h.batchOLD time24h.batchOLD
1                 1                0                0                0
2                 1                0                0                0
3                 1                0                0                0
4                 0                0                1                0
5                 0                0                1                0
6                 0                0                1                0
7                 0                0                1                0
8                 0                0                0                0
9                 0                0                0                0
10                0                0                0                0
11                0                1                0                0
12                0                1                0                0
13                0                1                0                0
14                0                0                0                1
15                0                0                0                1
16                0                0                0                1
17                0                0                0                1
18                0                0                0                0
19                0                0                0                0
20                0                0                0                0
   time04h.batchPRG time24h.batchPRG
1                 0                0
2                 0                0
3                 0                0
4                 0                0
5                 0                0
6                 0                0
7                 0                0
8                 1                0
9                 1                0
10                1                0
11                0                0
12                0                0
13                0                0
14                0                0
15                0                0
16                0                0
17                0                0
18                0                1
19                0                1
20                0                1
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$time
[1] "contr.treatment"

attr(,"contrasts")$batch
[1] "contr.treatment"

> fit.time <- lmFit(emat, design.time)

> contr.mat.time <- makeContrasts(time04h.batchNEW-time24h.batchNEW,
                                  time04h.batchOLD-time24h.batchOLD,
                                  time04h.batchPRG-time24h.batchPRG,
                                  levels=design.time)
> contr.mat.time
                 time04h.batchNEW - time24h.batchNEW
time04h.batchNEW                                   1
time24h.batchNEW                                  -1
time04h.batchOLD                                   0
time24h.batchOLD                                   0
time04h.batchPRG                                   0
time24h.batchPRG                                   0
                 time04h.batchOLD - time24h.batchOLD
time04h.batchNEW                                   0
time24h.batchNEW                                   0
time04h.batchOLD                                   1
time24h.batchOLD                                  -1
time04h.batchPRG                                   0
time24h.batchPRG                                   0
                 time04h.batchPRG - time24h.batchPRG
time04h.batchNEW                                   0
time24h.batchNEW                                   0
time04h.batchOLD                                   0
time24h.batchOLD                                   0
time04h.batchPRG                                   1
time24h.batchPRG                                  -1

> 
fit2.time <- contrasts.fit(fit.time, contr.mat.time)
fit2.time <- eBayes(fit2.time)

> colnames(fit2.time$coefficients)
[1] "time04h.batchNEW - time24h.batchNEW" "time04h.batchOLD -
time24h.batchOLD"
[3] "time04h.batchPRG - time24h.batchPRG

> length(which(topTable(fit2.time, number=12488, coef=1, adjust='fdr',
sort.by='P')[,3] <= 0.01))
[1] 184

> length(which(topTable(fit2.time, number=12488, coef=2, adjust='fdr',
sort.by='P')[,3] <= 0.01))
[1] 349

> length(which(topTable(fit2.time, number=12488, coef=3, adjust='fdr',
sort.by='P')[,3] <= 0.01))
[1] 27

What I see is that there are quite dramatic changes over time depending on
the batch. Actually the batches represent different labs performing an
expression analysis of *untreated* mouse cell cultures (MG-U74Av2 chip with
12488 probe sets) after 04h incubation and 24h incubation ... . The aim is to
see if the three labs have actually treated the cells the same ;-)

Does the above desing makes sense fot this kind of analysis?

	kind regards,

	Arne

--
Arne Muller, Ph.D.
Toxicogenomics, Aventis Pharma
arne dot muller domain=aventis com