Hi. I need to generate two GLM tests of a factorial design with RNA-Seq count data. I have 3 factors with 2 levels apiece (2 species X 2 treatments X 2 times), and 4 separate replicates each (i.e., we made a total of 2*2*2*4 = 32 separate libraries). Our main interest is in the interaction of species*treatment, as we think species A will alter gene expression in the treatment “stress” vs. treatment “benign”, whereas species B is expected to show little change. However, we’d like to also do another test of species*treatment*time, because it is possible that the ability of species A to alter gene expression in response to the “stress” treatment may differ at the 1st versus 2nd time point.

I think the way to set this up, is to create a design matrix as follows, with the lrt test with coef 5 giving the differentially expressed genes for the species*treatment test, and coef 8 giving the the differentially expressed gene for the species*treatment*time test (after calling topTags that is). Yet to ensure I have the statistics correct, my questions are: (1) is this thinking correct, as I don’t see many 3x2 factorial models to follow, and (2) do I need to set up a reference somehow (which I assume would be the set of four samples with TreatmentBenign*SpeciesB*Time2, but I’m not fully sure if that is correct or needed).

Many thanks in advance for your insight!
~Hilary

> designFF <- model.matrix(~Treatment*Species*Age)
> colnames(designFF)
[1] "(Intercept)"
[2] " TreatmentStress"
[3] "SpeciesA "
[4] "Time1"
[5] "TreatmentStress:SpeciesA"
[6] "TreatmentStress:Time1"
[7] "SpeciesA:Time1"
[8] "TreatmentStress:SpeciesA:Time1"

And then to run tests with:
> fit <- glmFit(y, designFF)

> lrtInteractionStressSpecies <- glmLRT(fitFF, coef=5)
> lrtInteractionStressSpeciesTime <- glmLRT(fitFF, coef=8)



	[[alternative HTML version deleted]]

