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

[Ryan C. Thompson]

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