[BioC] limma design: a basics questions for factorial analysis

[Marcelo Luiz de Laia]

Dear Biocondctors Users and Developers,

In a few days a go I send to this list a question about a design experiment.

I receive suggestion about that design in two way: the first suggestion 
indicated that design was a 2-way factorial. The second suggestion 
advice me to analyse it in limma and send me a script.

I start to learn about both suggestions. But, I see that my design have 
more one factor - time -, and this arrive a more few doubts.

I choose to analyse it in limma and, for that, I edit the previous 
script send me by a friend.

The background of our experiment:

Varietys: (V) 4
Treatments: 2 (treated [Tr] and not treated [Ct])
Time points: (T) 3
Biological replicates: 3
Genes: 3,575 printed 2 times
one color cDNA array

                        | - R1 - 7,150 spots
            | - Tr -- | - R2 - 7,150 spots
            |           | - R3 - 7,150 spots
     | T1-|
     |      |           | - R1  - 7,150 spots
     |      | -  Ct -- | - R2  - 7,150 spots
     |      |           | - R3  - 7,150 spots
V1 |      |
     |      |
   (......)

Then, I have a raw data from these situations:

V1T1TrR1- Var 1, Time 1, Treated, Rep 1
V1T1TrR2- Var 1, Time 1, Treated, Rep 2
V1T1TrR3- Var 1, Time 1, Treated, Rep 3
V1T1CtR1- Var 1, Time 1, Control, Rep 1
V1T1CtR2- Var 1, Time 1, Control, Rep 2
V1T1CtR3- Var 1, Time 1, Control, Rep 3
V1T2TrR1- Var 1, Time 2, Treated, Rep 1
(...)
V1T2CtR3- Var 1, Time 2, Control, Rep 3
V1T3TrR1- Var 1, Time 3, Treated, Rep 1
(...)
V1T3CtR3- Var 1, Time 3, Control, Rep 3
(...)
(...)
V4T3CtR3- Var 4, Time 3, Control, Rep 3

My data are in a matrix [7,150x72]
           [,1]              [,1]              [,1]              [,1]   
            (...)             [,72]
V1T1CtR1   V1T1TrR1   V1T1CtR2   V1T1TrR2              (...)   V4T3TrR3
1           10                  5                  7             100    
          (...)                90
1           10                  5                  7             100    
          (...)                90
2             7                  6                  9               98 
             (...)                80
(...)

-------------here start the script-----------------
# I edited his script for include the time factor
# because this changes, that I made, I was with doubts
# my doubts is not related with original script

library(limma)
y <- matrix(rnorm(7200)+6, ncol=72)
rownames(y) <- paste( "Gene", rep(1:50, each=2))
var <- factor(rep(LETTERS[1:4], each=18))
treat <- factor(rep(1:2, 18))
times <- factor(rep(1:3, each=6))
mVar <- matrix(0, 72, 4)
mTreat <- matrix(0, 72, 2)
mTimes <- matrix(0, 72, 3)
iV <- cbind(1:72, var)
iTr <- cbind(1:72, treat)
iTe <- cbind(1:72, times)
mVar[iV] <- 1
mTreat[iTr] <- 1
mTimes[iTe] <- 1
design <- cbind(geral=1, mVar[,-4], mTreat[,-2], mTimes[,-3])
correlacao <- duplicateCorrelation(y, design, ndups=2)
fit <- lmFit(y, design, ndups=2, correlation=correlacao$consensus)
fit2 <- eBayes(fit)
## Top 10: Var 4 vs. Var 1
topTable(fit2, coef=2)
## Top 10: Var 4 vs. Var 2
topTable(fit2, coef=3)
## Top 10: Treat 2 vs. Treat 1
topTable(fit2, coef=5)
## Top 10: Time 3 vs Time 1
topTable(fit2, coef=6)
volcanoplot( fit2, coef=5, highlight=4)
-----------script end here------------

My doubts:

With this design, is possible to compare Var 4 vs others ones. But, if I 
need to compare (Var 4 and Var 2) vs (Var 3 and Var 1) in all times 
together, only in time 1, only in time 2, only in time 3? Because Var 4 
and Var 2 are resistant and Var 3 and Var 1 are susceptible. What I need 
to change in design expression?

If I change the design in this way (mVar[,4] to mVar[,2],

design <- cbind(geral=1, mVar[,-2], mTreat[,-2], mTimes[,-3])

I will have:

coef=2 => Var 2 vs Var 1
coef=3 => Var 2 vs Var 3
coef=4 => Var 2 vs Var 4

This is correct?

In the same way I will can change time 3 by 2 or 1 in design expression 
end get comparison between 2 vs 1, and 2 vs 3?

I appreciate any commentaries about this one and/or about the best way 
to analyse this data set. I will love any suggestion.

I am very thanks to you!

Best wishes.

-- 
Marcelo Luiz de Laia
Ph.D Candidate
São Paulo State University (http://www.unesp.br/eng/)
School of Agricultural and Veterinary Sciences
Department of  Technology
Via de Acesso Prof.Paulo Donato Castellane s/n
14884-900   Jaboticabal - SP - Brazil
Fone: +55-016-3209-2675
Cell: +55-016-97098526