Hello List, I've got a bit of a complex design - perhaps you can offer me some advice. Thank you very much in advance for taking the time to read this mail. (I've highlighted my questions with *** if you just want to skip down) I have ~100 two-color spotted cDNA chips, making up a 20-timepoint timecourse, replicated with four different strains. (so ~20 timepoints/strain in a loop, plus a few hybs against a reference RNA). But a few timepoints have only three replicated strains because something went wrong in the lab. Until now, I've been inspired by the end of section 8.2 of limma_2.12.0's limmaUsersGuide(), as follows: ### below, my four strains are named A, D, G and H design <- modelMatrix(targets, ref = "REF") fit <- lmFit(MA, design) cont.matrix = makeContrasts(T_3h_vs0h= (A3h+D3h+G3h+K3h)/4 - (A0h+D0h +G0h+K0h)/4, T_6h_vs0h= (A6h+D6h+G6h)/3 - (A0h+D0h+G0h+K0h)/4, T_12h_vs0h= (A12h+D12h+G12h)/3 - (A0h+D0h+G0h+K0h)/4, T_18h_vs0h= (A18h+D18h+G18h+K18h)/4 - (A0h+D0h+G0h+K0h)/4, T_24h_vs0h= (A24h+D24h+G24h+K24h)/4 - (A0h+D0h+G0h+K0h)/4, T_36h_vs0h= (A36h+D36h+G36h+K36h)/4 - (A0h+D0h+G0h+K0h)/4, T_48h_vs0h= (A48h+D48h+G48h+K48h)/4 - (A0h+D0h+G0h+K0h)/4, T_3_vs0h= (A3+D3+G3+K3)/4 - (A0h+D0h+G0h+K0h)/4, T_4_vs0h= (A4+D4+G4+K4)/4 - (A0h+D0h+G0h+K0h)/4, T_5_vs0h= (A5+G5+K5)/3 - (A0h+D0h+G0h+K0h)/4, T_6_vs0h= (A6+D6+G6+K6)/4 - (A0h+D0h+G0h+K0h)/4, T_7_vs0h= (A7+D7+G7+K7)/4 - (A0h+D0h+G0h+K0h)/4, T_8_vs0h= (A8+D8+G8+K8)/4 - (A0h+D0h+G0h+K0h)/4, T_9_vs0h= (A9+D9+G9+K9)/4 - (A0h+D0h+G0h+K0h)/4, T_10_vs0h= (A10+D10+G10+K10)/4 - (A0h+D0h+G0h+K0h)/4, T_11_vs0h= (A11+D11+G11+K11)/4 - (A0h+D0h+G0h+K0h)/4, T_12_vs0h= (A12+D12+G12+K12)/4 - (A0h+D0h+G0h+K0h)/4, T_13_vs0h= (A13+D13+G13+K13)/4 - (A0h+D0h+G0h+K0h)/4, T_14_vs0h= (A14+D14+G14+K14)/4 - (A0h+D0h+G0h+K0h)/4, T_15_vs0h= (A15+D15+G15+K15)/4 - (A0h+D0h+G0h+K0h)/4, levels=design) fit2 <- contrasts.fit(fit, cont.matrix) Most RNA samples are used twice (eg: the "3hour" sample from strain A is used for both: A0h->A3h and: A3h->A6h), but some are used three times (this would be the case for A3h if I also had a A3h->REF chip). Limma is unable to estimate coefficients for many of my spots, ie fit2 $coefficients give "NA" rows for too many spots. If I understand correctly, this could be because for example both "K0h" chips are have 5000 flagged spots in common (indeed, several chips have thousands of genepix auto-flagged weak spots). NAs are frustrating. One workaround that came up when consulting my local statistician: to treat technical replicates and biological replicates equally. But technical replicates have higher correlation than biological replicates (I calculated a few correlations by hand): - technical replicate (same biological sample, same dye, different chip): median correlation = 0.55 - biological replicate (different sample, but same timepoint and dye, different chip): median correlation = 0.35 Thus I think it would be incorrect to treat my technical replicates as biological replicates. So I'm examining alternative means of creating a design matrix. Ideally, I'd like to use "strain" as a random effect here. *** Is it correct that the only way to add a random effect is with lmFit's block argument?: lmFit(myMA, design=myDesign, block=factor(targets$Colony), correlation=0.55) duplicateCorrelation cannot be used on a complex design. *** Is it correct to use the block argument on such a design, with a correlation I calculated "by hand"? *** Alternatively, I'll just put everything as fixed factors, but it feels wrong: extraFactors = model.matrix(~0+factor(targets$Strain)+factor(targets $Batch)) # several strains, only two possible batches colnames(extraFactors) = c(levels(factor(targets$Strain), "Batch") design_bothFixed = cbind(modelMatrix(targets, ref="REF"), extraFactors) *** Would the design matrix I just created be a correct way of proceeding? ***Since limma is limited to a single random factor, how should I choose which between Strain and Batch to use as fixed and which as random? Thanks very much :o) Yannick -------------------------------------------- yannick . wurm @ unil . ch Ant Genomics, Ecology & Evolution @ Lausanne http://www.unil.ch/dee/page28685_fr.html

maanova purports to handle several random effects. --Naomi At 05:01 PM 6/3/2008, Yannick Wurm wrote: >Hello List, > >I've got a bit of a complex design - perhaps you can offer me some >advice. Thank you very much in advance for taking the time to read >this mail. (I've highlighted my questions with *** if you just want >to skip down) > > > >I have ~100 two-color spotted cDNA chips, making up a 20-timepoint >timecourse, replicated with four different strains. (so ~20 >timepoints/strain in a loop, plus a few hybs against a reference >RNA). But a few timepoints have only three replicated strains because >something went wrong in the lab. > >Until now, I've been inspired by the end of section 8.2 of >limma_2.12.0's limmaUsersGuide(), as follows: > ### below, my four strains are named A, D, G and H > design <- modelMatrix(targets, ref = "REF") > fit <- lmFit(MA, design) > cont.matrix = makeContrasts(T_3h_vs0h= (A3h+D3h+G3h+K3h)/4 > - (A0h+D0h +G0h+K0h)/4, > T_6h_vs0h= (A6h+D6h+G6h)/3 - > (A0h+D0h+G0h+K0h)/4, > T_12h_vs0h= (A12h+D12h+G12h)/3 - > (A0h+D0h+G0h+K0h)/4, > T_18h_vs0h= (A18h+D18h+G18h+K18h)/4 > - (A0h+D0h+G0h+K0h)/4, > T_24h_vs0h= (A24h+D24h+G24h+K24h)/4 > - (A0h+D0h+G0h+K0h)/4, > T_36h_vs0h= (A36h+D36h+G36h+K36h)/4 > - (A0h+D0h+G0h+K0h)/4, > T_48h_vs0h= (A48h+D48h+G48h+K48h)/4 > - (A0h+D0h+G0h+K0h)/4, > T_3_vs0h= (A3+D3+G3+K3)/4 - > (A0h+D0h+G0h+K0h)/4, > T_4_vs0h= (A4+D4+G4+K4)/4 - > (A0h+D0h+G0h+K0h)/4, > T_5_vs0h= (A5+G5+K5)/3 - (A0h+D0h+G0h+K0h)/4, > T_6_vs0h= (A6+D6+G6+K6)/4 - > (A0h+D0h+G0h+K0h)/4, > T_7_vs0h= (A7+D7+G7+K7)/4 - > (A0h+D0h+G0h+K0h)/4, > T_8_vs0h= (A8+D8+G8+K8)/4 - > (A0h+D0h+G0h+K0h)/4, > T_9_vs0h= (A9+D9+G9+K9)/4 - > (A0h+D0h+G0h+K0h)/4, > T_10_vs0h= (A10+D10+G10+K10)/4 - > (A0h+D0h+G0h+K0h)/4, > T_11_vs0h= (A11+D11+G11+K11)/4 - > (A0h+D0h+G0h+K0h)/4, > T_12_vs0h= (A12+D12+G12+K12)/4 - > (A0h+D0h+G0h+K0h)/4, > T_13_vs0h= (A13+D13+G13+K13)/4 - > (A0h+D0h+G0h+K0h)/4, > T_14_vs0h= (A14+D14+G14+K14)/4 - > (A0h+D0h+G0h+K0h)/4, > T_15_vs0h= (A15+D15+G15+K15)/4 - > (A0h+D0h+G0h+K0h)/4, > levels=design) > fit2 <- contrasts.fit(fit, cont.matrix) > >Most RNA samples are used twice (eg: the "3hour" sample from strain A >is used for both: A0h->A3h and: A3h->A6h), but some are used three >times (this would be the case for A3h if I also had a A3h->REF chip). > >Limma is unable to estimate coefficients for many of my spots, ie >fit2 $coefficients give "NA" rows for too many spots. If I understand >correctly, this could be because for example both "K0h" chips are >have 5000 flagged spots in common (indeed, several chips have >thousands of genepix auto-flagged weak spots). > >NAs are frustrating. One workaround that came up when consulting my >local statistician: to treat technical replicates and biological >replicates equally. But technical replicates have higher correlation >than biological replicates (I calculated a few correlations by hand): > - technical replicate (same biological sample, same dye, different >chip): median correlation = 0.55 > - biological replicate (different sample, but same timepoint and >dye, different chip): median correlation = 0.35 >Thus I think it would be incorrect to treat my technical replicates >as biological replicates. > >So I'm examining alternative means of creating a design matrix. > >Ideally, I'd like to use "strain" as a random effect here. *** Is it >correct that the only way to add a random effect is with lmFit's >block argument?: > lmFit(myMA, design=myDesign, block=factor(targets$Colony), >correlation=0.55) >duplicateCorrelation cannot be used on a complex design. *** Is it >correct to use the block argument on such a design, with a >correlation I calculated "by hand"? *** > >Alternatively, I'll just put everything as fixed factors, but it >feels wrong: > extraFactors = > model.matrix(~0+factor(targets$Strain)+factor(targets $Batch)) # > several strains, only two possible batches > colnames(extraFactors) = c(levels(factor(targets$Strain), "Batch") > design_bothFixed = cbind(modelMatrix(targets, ref="REF"), > extraFactors) >*** Would the design matrix I just created be a correct way of >proceeding? > >***Since limma is limited to a single random factor, how should I >choose which between Strain and Batch to use as fixed and which as >random? > >Thanks very much :o) > >Yannick > > >-------------------------------------------- > yannick . wurm @ unil . ch >Ant Genomics, Ecology & Evolution @ Lausanne > http://www.unil.ch/dee/page28685_fr.html > >_______________________________________________ >Bioconductor mailing list >Bioconductor at stat.math.ethz.ch >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: >http://news.gmane.org/gmane.science.biology.informatics.conductor Naomi S. Altman 814-865-3791 (voice) Associate Professor Dept. of Statistics 814-863-7114 (fax) Penn State University 814-865-1348 (Statistics) University Park, PA 16802-2111