Question: 2 way anova in Bioconductor

0

Dana.Stanley@csiro.au •

**110**wrote: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

ADD COMMENT
• link
•
modified 8.1 years ago
by
Guillaume Meurice •

**210**• written 8.1 years ago by Dana.Stanley@csiro.au •**110**