DESeq: GLMs for RNA-Seq with interaction terms

[EMeyer]

I have a reasonably large RNA-Seq dataset for a non-model plant, and want to fit a GLM with interaction terms:
expression ~ treatment + time + treatment:time

In theory, I should be able to test the significance of each term by comparing a reduced model (dropping that term) with the full model above.

It seems DESeq can handle this for main effect terms: it produces a vector of p-values with a reasonable distribution. See this example, using data from the pasilla package:
data(pasillaGenes)
design <- pData(pasillaGenes)[, c("condition", "type")]
fullCountsTable <- counts(pasillaGenes)
cdsFull <- newCountDataSet(fullCountsTable[1:1000,], design)
cdsFull <- estimateSizeFactors(cdsFull)
cdsFull <- estimateDispersions(cdsFull, method="pooled")
fit1 <- fitNbinomGLMs(cdsFull, count ~ type + condition)
fit0 <- fitNbinomGLMs(cdsFull, count ~ type)
pvalsGLM <- nbinomGLMTest(fit1, fit0)


However, when I try this with an interaction term included it gives clearly wrong results: every p value is either 0 or 1.
fit1 <- fitNbinomGLMs(cdsFull, count ~ type + condition + type:condition)
fit0 <- fitNbinomGLMs(cdsFull, count ~ type + type:condition)
pvalsGLM <- nbinomGLMTest(fit1, fit0)

So my questions:
(a) Is this even possible in DESeq?
(b) If so, can anyone spot the error in how I coded this?

Thanks for any advice!

[Simon Anders]

Hi

you test for a rather odd contrast:

fit1 <- fitNbinomGLMs(cdsFull, count ~ type + condition + type:condition)
fit0 <- fitNbinomGLMs(cdsFull, count ~ type + type:condition)

The standard ANOVA or ANODEV approach is to test first for main effects without any interaction terms present in the fit, and then for the interaction, with all main effect present in both fits, i.e.,

fit1 <- fitNbinomGLMs(cdsFull, count ~ type + condition + type:condition)
fit0 <- fitNbinomGLMs(cdsFull, count ~ type + condition)

For the pasilla data, this shouldn't give many hits, though. (If there were interactions, this would be rather worrying.)

[raphael123]

Hi Simon !
Usually each coeficient in a GLM is associated with a p value, do you know how to get these p value from the GLM of DEseq ?