limma: model/contrast matrix for 2x3 factorial design
1
0
Entering edit mode
ckeeling • 0
@ckeeling-10042
Last seen 2.8 years ago

I want to analyze single-colour Agilent microarrays.  I have two Treatment groups (Treat & CTRL) and three Genotypes (WildType, MutantA, & MutantB), and each Treatment/Genotype combination has 4 replicates.

i.e.

Genotype <- factor(rep(c("WildType","MutantA","MutantB"),each=8),levels=c("WildType","MutantA", "MutantB"))
Treatment <- factor(rep(c("CTRL","Treat"),each=4,3),levels=c("CTRL","Treat"))

I want to test the following effects and comparisons:

Treatment Effect
Genotype Effect
Treatment x Genotype Interaction
WildType.Treat versus WildType.CTRL
MutantA.Treat versus MutantA.CTRL
MutantB.Treat versus MutantB.CTRL
MutantA.CTRL versus WildType.CTRL
MutantB.CTRL versus WildType.CTRL
MutantB.CTRL versus MutantA.CTRL

It is not clear to me how to set up the model.matrix and/or contrast.matrix to test all of these.

If I use:

design <- model.matrix(~Treatment*Genotype)

I get one coefficient for Treatment, but two for both Genotype and the interactions.  How do I combine coefficients to get a single Genotype effect and a single Treatment x Genotype interaction effect?

If I use:

design <- model.matrix(~0+Treatment*Genotype)

I get two coefficients for Treatment, and again two for both Genotype and the interactions.

I found the thread here that discusses 2x3: 2 way anova in Bioconductor
But, when I try it, the results from the interaction terms don't appear correct. I.e. the:

summary(decideTests(fit,method="global",adjust.method="BH",p.value=0.05,lfc=0))

shows only a handful of significant genes for each of the two Treatment x Genotype coefficients, but pulling out the results with:

c89 <- topTable(fit2, coef = 8:9, sort="none", n=Inf,adjust.method = "BH") # where 8 and 9 or the two Treatment x Genotype coefficients

gives 20,000 significant genes.

I cannot fully wrap my head around how to interact with the models and coefficients to get out what I want.  Any help would be appreciated.  Thanks!

2
Entering edit mode
@gordon-smyth
Last seen 7 hours ago
WEHI, Melbourne, Australia

Don't use the factorial model. It is better to treat your experiment as having 6 groups. See Section 9.5.2 of the limma User's Guide.

Note that the concepts of "Treatment effect" or "Genotype effect" are undefined in a two factor model, because you cannot sensibly summarize what one factor is doing without considering the other factor at the same time. You need to concentrate on particular contrasts instead.

0
Entering edit mode

Thanks Gordon for your prompt response!  So, something like this?

design = cbind(WildType.CTRL = c(1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), WildType.Treat = c(0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), MutantA.CTRL = c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0), MutantA.Treat = c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0), MutantB.CTRL = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0),MutantB.Treat = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1))

fit = lmFit(raw_aver_filtered, design)
contrastMatrix = makeContrasts("WildType.Treat-WildType.CTRL","MutantA.Treat-MutantA.CTRL","MutantB.Treat-MutantB.CTRL","MutantA.CTRL-WildType.CTRL","MutantB.CTRL-WildType.CTRL","MutantB.CTRL-MutantA.CTRL", levels=design)
fit2 = contrasts.fit(fit, contrastMatrix)
fit2 = eBayes(fit2)

summary(decideTests(fit2,method="global",adjust.method="BH",p.value=0.05,lfc=0))

So, it is not possible to have Treatment or Genotype effects, or their interaction?  This goes against my (albeit limited) understanding of multifactorial design.  How does one answer the question, "Which genes respond differently to treatment in the different genotypes?" in a 2x3 situation (compared to the Diff=(Mu.S-Mu.U)-(WT.S-WT.U)) described in section 9.5.2? Can I only get interactions between two at once?  I.e. Diff_A=(MutantA.Treat-MutantA.CTRL) - WildType.Treat-WildType.CTRL) and Diff_B=(MutantB.Treat-MutantB.CTRL) - WildType.Treat-WildType.CTRL)

2
Entering edit mode

Yes, you could use that design matrix.

You can test for any number of contrasts at a time. If you use:

topTable(fit2, coef=c("Diff_A","Diff_B"))

then you will get F-tests for interaction on 2 df. This is exactly the same as the test for interaction in the factorial model.

There are statistical things called treatment effects but they are scientifically meaningless. The whole purpose of your study to learn how the treatment effect depends on genotype. You wouldn't be doing this unless the treatment effect potentially depended on the genotype. So how could you define a "Treatment effect" without regard to genotype? This just wouldn't be meaningful.

If interactions turn out to be absent, then you could refit the simple Genotype+Treatment model, and then you would have unambiguous genotype and treatment effects. Otherwise, in the presence of interactions, you need to specify which genotype you are considering when you test for a treatment effect.

Traffic: 234 users visited in the last hour
FAQ
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.