[BioC] DESeq: GLMs for RNA-Seq with interaction terms

[Kasper Daniel Hansen]

First remove the interaction, then remove the main effect.  Or do both
at the same time.  Ie. either

  type+condition+type:condition
to
  type+condition
to
  type

Or just compare
  type+condition+type:condition
to
  type

Kasper

On Wed, Jul 13, 2011 at 5:39 PM, Eli Meyer <EliMeyer at mail.utexas.edu> wrote:
> Thanks for the fast response, Kasper.  No, thats not confusing, we were
> thinking the same thing early on... its often a question of figuring out how
> to express the model in the given software package.
> How would you recommend testing these main effect terms within DESeq, in
> that case?
>
> On Wed, Jul 13, 2011 at 3:17 PM, Kasper Daniel Hansen
> <kasperdanielhansen at gmail.com> wrote:
>>
>> On Wed, Jul 13, 2011 at 3:56 PM, Eli Meyer <EliMeyer at mail.utexas.edu>
>> wrote:
>> > 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)
>>
>> The smaller model makes no sense.  You are (in words) asking for an
>> interaction between condition and type without including a condition
>> (main) effect.  If you are confused about this (and I agree it might
>> be confusing at first depending on your training) you should read any
>> introduction to linear models.
>>
>> Kasper
>
>
>
> --
> Eli Meyer
> Postdoctoral Fellow
> Department of Integrative Biology
> University of Texas at Austin
> Austin, TX 78712
> office: (512) 471-3278
> cell: (310) 618-4483
>