Question: Two-way anova 2x2 design with limma
0
2.1 years ago by
damoros0
damoros0 wrote:

Hi,
In the 8.7 section (Two-factor example) in limma User Guide, explains three procedures that give the same result.
In the first procedure:

> cont.matrix <- makeContrasts(
+     SvsUinWT=WT.S-WT.U,
+     SvsUinMu=Mu.S-Mu.U,
+     Diff=(Mu.S-Mu.U)-(WT.S-WT.U),
+     levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)

How can I find the statistical F of two-way anova (genotype, treatment and interaction effect)? I need a factorial model formulas (second and third procedure)?

Thank you

modified 2.1 years ago • written 2.1 years ago by damoros0

You should read the limma User's Guide that comes with the version of the limma package you are using. Or else the latest version of the Guide is always available from:

https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

Your reference to Section 8.7 suggests that you are reading a limma User's Guide from many, many years ago.

Answer: Two-way anova 2x2 design with limma
1
2.1 years ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

There is no Section 8.7 in the user's guide. The example code seems to come from 9.5.2. If you want a two-way ANOVA, you will want to drop all non-intercept terms, so that the null is that all factor combinations have the same mean expression. This can be done by setting coef=2:4 in topTable with either of the nested interaction or in the classic interaction parameterisations of that specific experimental design.

Answer: Two-way anova 2x2 design with limma
0
2.1 years ago by
damoros0
damoros0 wrote:

Hi.

I understand that the procedure in 9.5.2 is the recommended one (single factor), but I am required to perform two way anova.

For example, if I extract nomalized intensities of one probe, I can perform a two way anova for that probe:

>mat
data       Strain    Treatment
1  15087.105     WT         U
2  13930.749     WT         U
3  16495.873     WT         U
4  13829.101     WT         U
5  12628.449     WT         S
6  10274.768     WT         S
7  11306.912     WT         S
8  11842.572     WT         S
9  16296.351     Mu         U
10 17066.216     Mu         U
11 15817.936     Mu         U
12 15001.087     Mu         U
13  9952.837     Mu         S
14 12524.801     Mu         S
15 10785.661     Mu         S
16 13326.560     Mu         S
> res.aov=aov(data~Strain*Treatment,data=mat)
> summary(res.aov)
Df   Sum Sq  Mean Sq F value   Pr(>F)
Strain            1  1806283  1806283   1.273    0.281
Treatment         1 59605572 59605572  41.998 3.02e-05 ***
Strain:Treatment  1  1156486  1156486   0.815    0.384
Residuals        12 17031010  1419251
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

it is right? I need to obtain the associated F and  p-values of strain, associated F and  p-values of treatment and associated F and  p-values of interaction effect in the same way that they are obtained from the aov function, using limma and its moderated F statistics. I understand that with topTable(fit2,coef=2:4) I get the results of global or omnibus F test (I'm not interested in that test).

I understand that I should use the classic interaction model (9.5.4 section):

> design <- model.matrix(~Strain*Treatment)

(Intercept) StrainWT TreatmentU StrainWT:TreatmentU
1            1        1          1                   1
2            1        1          1                   1
3            1        1          1                   1
4            1        1          1                   1
5            1        1          0                   0
6            1        1          0                   0
7            1        1          0                   0
8            1        1          0                   0
9            1        0          1                   0
10           1        0          1                   0
11           1        0          1                   0
12           1        0          1                   0
13           1        0          0                   0
14           1        0          0                   0
15           1        0          0                   0
16           1        0          0                   0

> fit <- lmFit(eset, design)
> fit2 <- eBayes(fit2)

How to continue?

2. F-statistics and p-values can be obtained with appropriate settings of coef in topTable. For example, coef=2 will be equivalent to Strain, coef=3 will be equivalent to Treatment, and coef=4 will equivalent to testing the interaction effect. (I say equivalent, not identical, due to the effect of empirical Bayes shrinkage in limma that isn't available with aov.)
3. Having said that, keep in mind that StrainWT and TreatmentU do not represent the main effects of the strain or treatment when the interaction term is non-zero. When there is an interaction between the factors, StrainWT represents the log-fold change of WT over Mu in treatment S only, while TreatmentU represents the log-fold change of U over S in the Mu strain only. The same applies for the output of the aov call, which also uses model.matrix internally.

Thank you Aaron. I understand that the three t-statistic obtained with topTable in this case is equivalent to the three F statistics of aov function (and p-values).

Two last questions:

a)  In the case of a two-color analysis using the separate channels procedure, would this be the procedure to follow?

> corfit <- intraspotCorrelation(MA, design)

> fit <- lmscFit(MA, design, correlation=corfit\$consensus)

> fit2 <- eBayes(fit)

> topTable(fit2, coef=2, adjust="BH") # Strain

> topTable(fit2, coef=3, adjust="BH") # Treatment

> topTable(fit2, coef=4, adjust="BH") #Interaction

b) Using separate channels procedure, should I include the dye effect in the design matrix? In the user guide it is not included, but I have seen it in other posts and tutorials. Biological replicates have been labeled alternately with cy5 and cy3 (two with cy5 and two with cy3).