[BioC] 2 way anova in Bioconductor
Jenny Drnevich
drnevich at illinois.edu
Wed Jun 29 16:31:49 CEST 2011
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
>
>
> [[alternative HTML version deleted]]
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>
>
>The information of this email and in any file transmitted with it is
>strictly confidential and may be legally privileged.
>It is intended solely for the addressee. If you are not the intended
>recipient, any copying, distribution or any other use of this email is
>prohibited and may be unlawful. In such case, you should please notify the
>sender immediately and destroy this email.
>The content of this email is not legally binding unless confirmed by
>letter.
>Any views expressed in this message are those of the individual sender,
>except where the message states otherwise and the sender is authorised to
>state them to be the views of the sender's company. For further
>information about Actelion please see our website at
>http://www.actelion.com
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>
>
>The information of this email and in any file transmitted with it is
>strictly confidential and may be legally privileged.
>It is intended solely for the addressee. If you are not the intended
>recipient, any copying, distribution or any other use of this email
>is prohibited and may be unlawful. In such case, you should please
>notify the sender immediately and destroy this email.
>The content of this email is not legally binding unless confirmed by letter.
>Any views expressed in this message are those of the individual
>sender, except where the message states otherwise and the sender is
>authorised to state them to be the views of the sender's company.
>For further information about Actelion please see our website at
>http://www.actelion.com
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list