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

[Gordon K Smyth]

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}}