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!
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)andDiff_B=(MutantB.Treat-MutantB.CTRL) - WildType.Treat-WildType.CTRL)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.