[BioC] Multifactor analysis with deseq

[Simon Anders]

Hi Michel

On 22/03/13 20:52, Michael Muratet wrote:
> fit0 <- fitNbinomGLMs(cds, count~genotype+day+culture)
 > fit1 <- fitNbinomGLMs(cds, count~day+culture)
> then
 > pvalsGLM = nbinomGLMTest( fit1, fit0 )
>
> The result was entirely NaN or NAs.

This is easy. :-) It is simply because you supplied the fits in the 
wrong oder. The fit with 'genotype' is the full model and should come as 
first argument to nbinomGLMTest, and the reduced model without 
'genotype' is second.  So, just swap fit0 and fit1.

> Does the formula in the fitNbinomGLMs() function work the same way as
> lm(), e.g., do I have to explicitly include interactions or will it
> calculate them for me?

Yes, both fitNbinomGLMs and lm internally call model.matrix, which, as 
you know, adds an intercept automatically

> Is there a way to use relevel() to set the reference level? I have a
> wild-type genotype that would be natural for comparing two mutant
> genotypes.

Yes, you can use 'relevel' in the same way as you would with 'lm'.

And to answer your question about the 'fixme'. There are two 
possibilities: You can compare a single genotype with wildtype by simply 
subsetting the CountDataSet to exclude the samples for the othe 
rgenotype before doing the fits. Or you can use our new DESeq2, which 
now offers a Wald test instead of a likelihood ratio test and can so 
give you p values for both genotype coefficient in one go.

> I don't see any references on the list, but has anybody tried to use
> normalized values from a CountDataSet outside of deseq?

Depends on what you want to do with them.

> Lastly, I saw on the list while researching this issue that deseq2 is
> out. I'll give it try.

Sure, do so. And tell us what you think.

   Simon