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!
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!
Comment