I am doing comparisons of edgeR, DESeq, and limma for RNA-Seq data (yes, I know limma is not really recommended for RNA-Seq, but I'd like to benchmark them). My experimental design has an experimental factor (Timepoint) and a blocking factor (Donor), something like this:
> expdata
Donor Timepoint
L3_E1 4659 T0
L3_E2 4659 T120
L3_E5 5131 T0
L3_E6 5131 T120
L4_E1 4659 T0
L4_E2 4659 T120
L4_E5 5131 T0
L4_E6 5131 T120
L5_E9 5291 T0
L5_E10 5291 T120
L5_E13 5053 T0
L5_E14 5053 T120
L6_E9 5291 T0
L6_E10 5291 T120
L6_E13 5053 T0
L6_E14 5053 T120
Which you can construct using this code (obtained via `deparse
expdata <- structure(list(Donor = structure(c(1L, 1L, 3L, 3L, 1L, 1L, 3L, 3L, 4L, 4L, 2L, 2L, 4L, 4L, 2L, 2L), .Label = c("4659", "5053", "5131", "5291"), class = "factor"), Timepoint = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("T0", "T120"), class = "factor")), .Names = c("Donor", "Timepoint"), row.names = c("L3_E1", "L3_E2", "L3_E5", "L3_E6", "L4_E1", "L4_E2", "L4_E5", "L4_E6", "L5_E9", "L5_E10", "L5_E13", "L5_E14", "L6_E9", "L6_E10", "L6_E13", "L6_E14"), class = "data.frame")
The way to do this analysis in DESeq and edgeR is to fit a model like ~ Donor + Timepoint and then compare it to a fit of the null model ~ Donor:
## DESeq
alt.formula <- ~Donor + Timepoint
null.formula <- ~Donor
cds <- newCountDataSet(feature.counts, expdata)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds, modelFormula=alt.formula)
fit0 <- fitNbinomGLMs(cds, null.formula)
fit1 <- fitNbinomGLMs(cds, alt.formula)
pvalsGLM <- nbinomGLMTest(fit1, fit0)
padjGLM <- p.adjust(pvalsGLM, method="BH")
results <- data.frame(row.names=featureNames(cds),
LR=fits$null$deviance - fits$alt$deviance,
PValue=pvalsGLM, FDR=padjGLM)
## edgeR
design <- model.matrix(~Donor + Timepoint, data=expdata)
dge <- DGEList(counts=feature.counts,
group=groupvar,
genes=feature.annot)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
lrt <- glmLRT(dge, fit, coef="TimepointT120")
results <- topTags(lrt, n=nrow(lrt))
Now, I can also do the same thing in limma, adding "Donor" as a term in my model formula:
## Limma with Donor in the model
design <- model.matrix(~Donor + Timepoint, data=expdata)
nf <- calcNormFactors(feature.counts)
elist <- voom(feature.counts, design, plot=FALSE,
lib.size=colSums(feature.counts)*nf)
fit <- lmFit(elist, design)
fit <- eBayes(fit)
results <- topTable(fit, coef="TimepointT120",
n=nrow(elist), sort.by="none")
which I understand produces a moderated paired t-test. However, in the limma docs and also in various mailing list postings (e.g. this one) , I see recommentations to use the block argument to limma's lmFit function along with the duplicateCorrelation function:
## Limma with Donor as blocking factor
design <- model.matrix(~Timepoint, data=expdata)
nf <- calcNormFactors(feature.counts)
elist <- voom(feature.counts, design, plot=FALSE,
lib.size=colSums(feature.counts)*nf)
corfit <- duplicateCorrelation(elist, design, block=expdata$Donor)
fit <- lmFit(elist, design, block=expdata$Donor, correlation=corfit$consensus.correlation)
results <- topTable(fit, coef="TimepointT120",
n=nrow(elist), sort.by="none")
which is claimed to give better statistical power. Indeed, when I run both versions of the limma code, the second produces roughly the same gene list (sorted by p-value) only with much more significant p-values.
So now, I am confused because people who are using DESeq and edgeR are saying that you should treat "Donor" as a blocking factor by adding it as a term in your model, while for limma doing exactly the same thing is apparently suboptimal with respect to statistical power. And as far as I can tell, DESeq and edgeR don't have any direct equivalent of limma's "block" argument or duplicateCorrelation function. So is it really true that the analogous procedure (i.e. using ~Donor+Timepoint as the model formula) is correct in edgeR and DESeq but suboptimal in limma? If so, that seems like an odd contrast in what are otherwise quite similar analysis procedures.
So what is the difference between the two blocks of limma code above? Is the second block still doing a paired test? Which of the two blocks is really a closer match to the edgeR and DESeq blocks above?
Very interesting question which I would also want to know the answer to. Maybe ask on the Bioconductor mailing list, where edgeR and limma questions are usually answered?