20 Nov 2012 02:26
Equivalent of contrasts.fit & multi-contrast decideTests for edgeR?
Ryan C. Thompson <rct@...>
2012-11-20 01:26:37 GMT
2012-11-20 01:26:37 GMT
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 _______________________________________________ Bioconductor mailing list Bioconductor@... https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor