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

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