Thank you all for your help, especially Jenny for such a detailed explanation!!! I followed your email and got it to work!!! Thanks again Dana -----Original Message----- From: Natasha Sahgal [mailto:nsahgal@well.ox.ac.uk] Sent: Thursday, 30 June 2011 7:27 PM To: Jenny Drnevich Cc: Stanley, Dana (LI, Geelong AAHL); bioconductor at r-project.org; bioconductor-bounces at r-project.org Subject: Re: [BioC] 2 way anova in Bioconductor Hi Jenny and the BioC list, Please correct me if I am wrong here, but isn't the contrast "Female - Male" rather than "Male - Female" as you say? I was of the understanding that the contrasts, when not specified by contrasts.matrix, by default is the latter - former (e.g. B - A), in a 2 group comparison (say when levels = c(A,B)). BW, Natasha On 29/06/2011 15:31, Jenny Drnevich wrote: > Actually, you can do ANY sort of factorial design in limma. It took me > awhile to figure out how to expand the sum to zero parametrization > (3rd example, Section 8.7 limmaUsersGuide() ) when you have more than > 2 levels of any factor. Here's an example of a 2x3 design: > > #First set up the 2 group factor, just like in the limma vignette: > > > Sex <- factor(rep(c("Female","Male"),each=6),levels=c("Male","Female")) > # Note that the order you specify the groups in the levels argument > determines the direction of the comparison. See below. > > > contrasts(Sex) <- contr.sum(2) > > Sex > [1] Female Female Female Female Female Female Male Male Male > Male Male Male > attr(,"contrasts") > [,1] > Male 1 > Female -1 > Levels: Male Female > # Note that the contrast is male - female because female was listed > last in the levels argument above > > #Now set up the 3 group factor: > > > Time <- factor(rep(1:3,4),levels=3:1) > # Again, the group order specified in the levels determines > the direction of comparison > > > contrasts(Time) <- contr.sum(3) > # You use contr.sum(3) here because you have 3 groups. If you > have 4 groups, you'd use 4, etc. > > > Time > [1] 1 2 3 1 2 3 1 2 3 1 2 3 > attr(,"contrasts") > [,1] [,2] > 3 1 0 > 2 0 1 > 1 -1 -1 > Levels: 3 2 1 > # Note the contrasts are 3 - 1 and 2 - 1, because group 1 was listed > last in the levels argument above > > > design <- model.matrix(~Sex*Time) > > design > (Intercept) Sex1 Time1 Time2 Sex1:Time1 Sex1:Time2 > 1 1 -1 -1 -1 1 1 > 2 1 -1 0 1 0 -1 > 3 1 -1 1 0 -1 0 > 4 1 -1 -1 -1 1 1 > 5 1 -1 0 1 0 -1 > 6 1 -1 1 0 -1 0 > 7 1 1 -1 -1 -1 -1 > 8 1 1 0 1 0 1 > 9 1 1 1 0 1 0 > 10 1 1 -1 -1 -1 -1 > 11 1 1 0 1 0 1 > 12 1 1 1 0 1 0 > attr(,"assign") > [1] 0 1 2 2 3 3 > attr(,"contrasts") > attr(,"contrasts")$Sex > [,1] > Male 1 > Female -1 > > attr(,"contrasts")$Time > [,1] [,2] > 3 1 0 > 2 0 1 > 1 -1 -1 > > # In this design matrix, the (Intercept) coefficient is the grand > mean, the Sex1 coef is the main effect of Sex, > # the Time1 and Time2 coef taken together will give you the main > effect of Time, and the Sex1:Time1 and > # Sex1:Time2 coef taken together will give you the interaction term. > > # Now fit the model with your data > > > fit.2x3 <- eBayes(lmFit(YourData, design)) > > # To get the overall F test for the 2x3 take all coef except the > Intercept: > > > overall.2x3 <- topTable(fit, coef=2:6, number=Inf) > > # To get the F test (equivalent to t test) main effect of Sex, do: > > > mainSex <- topTable(fit, coef=2, number=Inf) > # Note that the logFC values need to be multiplied by 2 to get the > actual Male:Female logFC value! > > # To get the F test for the main effect of Time, do: > > > mainTime <- topTable(fit, coef=3:4, number=Inf) > # I usually don't use the actual logFC values listed for each coef, so > I usually don't care which group is listed last > > # To get the F test for the interaction term, do: > > > Interact <- topTable(fit, coef=5:6, number=Inf) > > > If you have 4 groups in a level, you'll have 3 coef for the main > effect of the level and 3 coef for the interaction term, and so on up > the line. You can also do n-way factorial designs in the same way. > Just make sure to keep track of the coefficients! > > Cheers, > Jenny > > At 03:39 AM 6/29/2011, axel.klenk at actelion.com wrote: >> Hi Dana, >> >> ooops, obviously I shouldn't be posting immediately after lunch... >> re-reading >> your message after some cups of coffee I realize that my brain has >> skipped >> everything after male/female and treatment/control... so yours is >> clearly >> NOT >> a 2x2 design... my apologies... >> >> Cheers, >> >> - axel >> >> >> Axel Klenk >> Research Informatician >> Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil / >> Switzerland >> >> >> >> >> From: >> axel.klenk at actelion.com >> To: >> <dana.stanley at="" csiro.au=""> >> Cc: >> bioconductor at r-project.org, bioconductor-bounces at r-project.org >> Date: >> 28.06.2011 15:13 >> Subject: >> Re: [BioC] 2 way anova in Bioconductor >> Sent by: >> bioconductor-bounces at r-project.org >> >> >> >> Hi Dana, >> >> limma can do that easily. The current version of the user's guide >> contains one chapter and two case studies on 2x2 factorial designs. >> >> Cheers, >> >> - axel >> >> >> Axel Klenk >> Research Informatician >> Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil / >> Switzerland >> >> >> >> >> >> From: >> <dana.stanley at="" csiro.au=""> >> To: >> <bioconductor at="" r-project.org=""> >> Date: >> 28.06.2011 03:27 >> Subject: >> [BioC] 2 way anova in Bioconductor >> Sent by: >> bioconductor-bounces at r-project.org >> >> >> >> Hi Everyone >> >> I use custom designed chicken high-density multiplex one channel >> microarray (NimbleGen 12x135K). My usual pipeline includes mostly oligo, >> arrayQualityMetrics, preprocessCore, genefilter, GeneSelector and limma. >> So far I only worked with simple design, 2 conditions; control/treatment >> style. Now I have a dataset with embryonic development timeline for male >> and female chicks. I still compare different timepoints using limma >> but I >> want to do 2 way anova to identify interaction significant genes, ie >> genes >> >> where gender influences development timeline. I looked at many packages >> and could not find 2 way anova, though I am sure it is there >> somewhere. I >> tested package ABarray that seemed ideal for my me but after days of >> trying I contacted the author to find out that expression sets (mine is >> coming from limma) can no longer be used by ABarray as the package has >> not been updated for a while. Can anyone suggest a package (and an >> example code if possible) for 2 way anova? I am so temp! >> ted to go and do it in GeneSpring... >> >> Thanks all >> >> Dana

Dear All, I have a few more question regarding the design matrix bellow (mainly How to interpret the coefficient ?) > > design > (Intercept) Sex1 Time1 Time2 Sex1:Time1 Sex1:Time2 > 1 1 -1 -1 -1 1 1 > 2 1 -1 0 1 0 -1 > 3 1 -1 1 0 -1 0 > 4 1 -1 -1 -1 1 1 > 5 1 -1 0 1 0 -1 > 6 1 -1 1 0 -1 0 > 7 1 1 -1 -1 -1 -1 > 8 1 1 0 1 0 1 > 9 1 1 1 0 1 0 > 10 1 1 -1 -1 -1 -1 > 11 1 1 0 1 0 1 > 12 1 1 1 0 1 0 > attr(,"assign") > [1] 0 1 2 2 3 3 > attr(,"contrasts") > attr(,"contrasts")$Sex > [,1] > Male 1 > Female -1 > > attr(,"contrasts")$Time > [,1] [,2] > 3 1 0 > 2 0 1 > 1 -1 -1 > > # In this design matrix, the (Intercept) coefficient is the grand mean, the Sex1 coef is the main effect of Sex, > # the Time1 and Time2 coef taken together will give you the main effect of Time, and the Sex1:Time1 and > # Sex1:Time2 coef taken together will give you the interaction term. Could someone give more explanation of the interaction terms ? Does such terms catch the variability due to both of these effects ? In the example given above, let's imagine the time effect is very difficult to detect because the signal is highly variable due to Sex effect. Setting the interaction terms in the design matrix allow to discard this variability and then to get genes significantly differentially expressed due to time 1 (using the time 1 coefficient) ? And in this case, the logFC obtained are the good one ? > # Now fit the model with your data > > > fit.2x3 <- eBayes(lmFit(YourData, design)) > # To get the overall F test for the 2x3 take all coef except the Intercept: > > overall.2x3 <- topTable(fit, coef=2:6, number=Inf) Does this raise the features which are a significantly differentially expressed in one coefficient at least ? > > mainSex <- topTable(fit, coef=2, number=Inf) > # Note that the logFC values need to be multiplied by 2 to get the actual Male:Female logFC value! I didn't understand why you need to multiply by 2 here ? The coefficient didn't already give the correct logFC values ? many thanks for the help. -- GM