[BioC] Equivalent of contrasts.fit & multi-contrast decideTests for edgeR?

[Gordon K Smyth]

Dear Ryan,

The edgeR glmLRT() function now allows the contrast to be a matrix with 
multiple columns.  The behavior of this argument is now analogous to the 
coef argument.  When multiple contrasts (or multiple coefs) are specified, 
the LR tests are conducted for all the contrasts at once, i.e., the LRTs 
are on several degrees of freedom, one for each independent contrast.

This extension allows very general access to anova-type tests using edgeR.

I committed this change to the official release of edgeR yesterday.

Best wishes
Gordon

---------------------------------------------
Professor Gordon K Smyth,
Bioinformatics Division,
Walter and Eliza Hall Institute of Medical Research,
1G Royal Parade, Parkville, Vic 3052, Australia.
http://www.statsci.org/smyth

On Tue, 20 Nov 2012, Gordon K Smyth wrote:

> Dear Ryan,
>
> It is true that glmLRT() does not allow you to specify multiple contrast 
> vectors to get an F-like test.  There is no good reason for this other than 
> that we haven't got around to it yet.  We will by the next Bioconductor 
> release.
>
> The fact that decideTestsDGE doesn't work with multiple coefficients is 
> unavoidable.  Note that decideTests() in limma doesn't give results for 
> F-tests either, only for individual contrasts.  The basic difference between 
> linear models and glms is that with linear models, the F-test and all the 
> individual t-tests that make it up can be conducted as one computation, but 
> for glms each individual 1df test and the overall F-test have to all be 
> separate computations.  This is why we ask you to do multiple calls to 
> glmLRT() to test for lots of different contrasts, whereas limma can test for 
> any number of contrasts in one step.
>
> Best wishes
> Gordon
>
> --------------- original message --------------
>
> [BioC] Equivalent of contrasts.fit & multi-contrast decideTests for edgeR?
> Ryan C. Thompson rct at thompsonclan.org
> Tue Nov 20 02:26:37 CET 2012
>
> Hi all,
>
> I am trying to compare limma (with voom) and edgeR for RNA-seq
> differential expression analysis, and I have noticed that while edgeR's
> glm functionality closely matches the functionality of limma, one
> feature seems to be missing: testing of multiple contrasts. Specifically:
>
> 1. In glmLRT, the contrast argument only takes a single contrast, not a
> matrix of contrasts (as limma's contrasts.fit would);
> 2. If glmLRT is used with a coef argument containing 2 or more coefs,
> then decideTestsDGE cannot handle the resulting object.
>
> To illustrate what I mean with an example, consider the following
> experimental design with 3 replicates each of 3 timepoints, where I fit
> the same data with two equivalent design matrices, one with an intercept
> term and one without:
>
> library(edgeR)
> dge <- DGEList(...) # Imagine data here
> sampledata <- data.frame(timepoint=,
>                          )
> timepoint <- rep(factor(c("T0", "T1", "T2", "T3")), each=3)
> design <- model.matrix(~timepoint)
> design.noint <- model.matrix(~0+timepoint)
> fit <- glmFit(dge, design)
> fit.noint <- glmFit(dge, design.noint)
> ## Test for changes in any timepoint
> lrt.any.changes <- glmLRT(fit, coef=c(2,3,4))
> ## How can this test be performed on fit.noint?
> lrt.any.changes <- glmLRT(fit.noint, ???)
> ## This throws an error because the DGELRT has multiple columns
> ## "logFC.*" instead of just a single "logFC" that the function
> ## expects.
> decideTestsDGE(lrt.any.changes)
>
>
> By contrast, with limma I can always do the test that I want regardless
> of how I choose to parametrize my design matrix (intercept or not):
>
> library(limma)
> ## Equivalent procedure in limma (I think)
> lfit <- lmFit(voom(dge, design))
> lfit.noint <- lmFit(voom(dge, design.noint))
> ## Test for changes in any timepoint (result is in $F.p.value)
> lfit <- eBayes(lfit)
> ## Same test on the version with no intercept term
> contrasts.anychange.noint <- makeContrasts(timepointT1-timepointT0,
>                                            timepointT2-timepointT0,
>                                            timepointT3-timepointT0,
>                                            levels=design.noint)
> lfit.noint <- eBayes(contrasts.fit(lfit.noint))
> ## Should give identical results?
> decideTests(lfit)
> decideTests(lfit.noint)
>
> Basically, I far as I can tell, with edgeR you can test the null
> hypothesis of multiple model coefficients being zero, but not multiple
> contrasts, despite the fact that both procedures should be statistically
> equivalent. Is edgeR missing this functionality or am I missing the
> proper way to do it? Not having this functionality makes things a little
> confusing, because depending on which one of several equivalent
> parametrizations I choose, different tests are available or not
> available, as illustrated by the code above, in which I can only test
> the hypothesis of "any change between any time points" if I include an
> intercept term. If I'm missing something, can someone pleas enlighten
> me? If edgeR really is missing this functionality, is it planned for the
> future or is there some fundamental difference between lms and glms that
> makes it impossible?
>
> Thanks,
>
> -Ryan Thompson
>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}