Re: DESeq2: df error message with nbinomial LRTest

Dear Mike,

Thanks for explanations, I will follow your recommandations.

Best regards
Marie-Agnès



Le 31/08/2014 15:16, Michael Love a écrit :
>
> Hi Marie-Agnès,
> On Aug 31, 2014 12:16 PM, "INRA - coutellec" 
> <[email protected] 
> <mailto:[email protected]>> 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
> > [email protected] <mailto:[email protected]>
> > https://stat.ethz.ch/mailman/listinfo/bioconductor
> > Search the archives: 
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>


-- 
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
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor