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