[BioC] DESeq2: df error message with nbinomial LRTest

[Michael Love]

Hi Marie-Agnès,
On Aug 31, 2014 12:16 PM, "INRA - coutellec" <
marie-agnes.coutellec at rennes.inra.fr> wrote:
>
> Dear all,
>
> My question is about the right design to use to perform LRT tests with
> DESeq2. I have two crossed factors, tra (2 levels) and pop (4 levels). I
> test the interaction term as follows:
>
> >  ddsFullCountTable <- DESeqDataSetFromMatrix(
>
> + countData = countE,
>
> + colData = cda,
>
> + design = ~ pop*tra)
>
>
>  > dds2inter.LRT <- DESeq(ddsFullCountTable, test="LRT", full=~pop*tra,
> reduced=~pop+tra)
>  > res2<-results(dds2.LRT)
>  > sum( res2$padj < 0.1, na.rm=TRUE )
> [1] 389
>
> and get 389 DEGs (out of 45249 contigs). Although this is very fine with
> me, I would like to know what  exactly is listed in column
> 'log2foldchange' for this interaction term (change between what and
what?).

See the help for ?results, where we describe the meaning of the LFC column
of a results object for test="LRT".

The LFC is for just 1 of the 3 interaction terms (the last level of both
factors), while the p-value is for all interaction terms.

As in linear models, an interaction term is an additional term to account
for differences in the tra effect across pop groups. It is an additive term
in the log2 transformed space, so an additional fold change in the space of
counts.

>
> Next, I want to test for the "tra" effect, alone. I tried :
>
>  > dds.LRT <- DESeq(ddsFullCountTable, test="LRT", full=~pop*tra,
> reduced=~pop+pop:tra)
> or:
>  > dds.LRT <- DESeq(ddsFullCountTable, test="LRT",
> full=~pop+tra+pop:tra, reduced=~pop+pop:tra)
>
> and get the same error message:
> Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior =
> betaPrior,  :
>    less than one degree of freedom, perhaps full and reduced models are
> not in the correct order
>
> whereas the command below works, although the two models differs by two
> elements (which shouldn't work in my view):
>  > dds.LRT <- DESeq(ddsFullCountTable, test="LRT",
> full=~pop+tra+tra:pop, reduced=~pop)
>
>

For the main 'tra' effect, you should just use the Wald test:

results(dds, contrast=c("tra","B","A"))

Filling in the two levels of tra.

Or if this doesn't work (i don't know what version you are using), run this
first:

dds = nbinomWaldTest(dds, betaPrior=FALSE)

Mike

>
> I understand it is somehow meaningless to eliminate a factor while
> keeping it through its interaction in a model. Thus, I plan to:
> 1) take out the 389 "interaction" DEGs
> 2) reestimate size factors and dispersion with an additive design
> (design=~pop+tra), from which I would be able to test the effect of
> "tra" and of "pop" very easily.
>
> Can anyone tell me if this is correct, or if I alternatively should keep
> "pop*tra" as design to estimate dispersion and then use "pop+tra" as
> full model and "pop" as reduced model in my LRT.
>
> Many thanks in advance,
> Marie-Agnès
>
> --
> Marie-Agnès Coutellec
> INRA/Agrocampus-Ouest UMR0985 Ecology and Ecosystem Health
> Ecotoxicology and Quality of Aquatic Environments
> 65 rue de Saint-Brieuc
> 35042 Rennes cedex - FRANCE
> phone: +33 223 485 248
> http://www6.rennes.inra.fr/ese_eng/
>
>
>         [[alternative HTML version deleted]]
>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor

	[[alternative HTML version deleted]]