<p>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:</p>

<pre><code>&gt; 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
</code></pre>

<p>Which you can construct using this code (obtained via `deparse</p>

<pre><code>expdata &lt;- 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")
</code></pre>

<p>The way to do this analysis in DESeq and edgeR is to fit a model like <code>~ Donor + Timepoint</code> and then compare it to a fit of the null model <code>~ Donor</code>:</p>

<pre><code>## DESeq
alt.formula &lt;- ~Donor + Timepoint
null.formula &lt;- ~Donor
cds &lt;- newCountDataSet(feature.counts, expdata)
cds &lt;- estimateSizeFactors(cds)
cds &lt;- estimateDispersions(cds, modelFormula=alt.formula)
fit0 &lt;- fitNbinomGLMs(cds, null.formula)
fit1 &lt;- fitNbinomGLMs(cds, alt.formula)
pvalsGLM &lt;- nbinomGLMTest(fit1, fit0)
padjGLM &lt;- p.adjust(pvalsGLM, method="BH")
results &lt;- data.frame(row.names=featureNames(cds),
                      LR=fits$null$deviance - fits$alt$deviance,
                      PValue=pvalsGLM, FDR=padjGLM)

## edgeR
design &lt;- model.matrix(~Donor + Timepoint, data=expdata)
dge &lt;- DGEList(counts=feature.counts,
               group=groupvar,
               genes=feature.annot)
dge &lt;- estimateGLMCommonDisp(dge, design)
dge &lt;- estimateGLMTrendedDisp(dge, design)
dge &lt;- estimateGLMTagwiseDisp(dge, design)
fit &lt;- glmFit(dge, design)
lrt &lt;- glmLRT(dge, fit, coef="TimepointT120")
results &lt;- topTags(lrt, n=nrow(lrt))
</code></pre>

<p>Now, I can also do the same thing in limma, adding "Donor" as a term in my model formula:</p>

<pre><code>## Limma with Donor in the model
design &lt;- model.matrix(~Donor + Timepoint, data=expdata)
nf &lt;- calcNormFactors(feature.counts)
elist &lt;- voom(feature.counts, design, plot=FALSE,
              lib.size=colSums(feature.counts)*nf)
fit &lt;- lmFit(elist, design)
fit &lt;- eBayes(fit)
results &lt;- topTable(fit, coef="TimepointT120",
                    n=nrow(elist), sort.by="none")
</code></pre>

<p>which I understand produces a moderated paired t-test. However, in the limma docs and also in various mailing list postings (e.g. <a href='http://permalink.gmane.org/gmane.science.biology.informatics.conductor/6622'>this one</a>) , I see recommentations to use the block argument to limma's lmFit function along with the duplicateCorrelation function:</p>

<pre><code>## Limma with Donor as blocking factor
design &lt;- model.matrix(~Timepoint, data=expdata)
nf &lt;- calcNormFactors(feature.counts)
elist &lt;- voom(feature.counts, design, plot=FALSE,
              lib.size=colSums(feature.counts)*nf)
corfit &lt;- duplicateCorrelation(elist, design, block=expdata$Donor)
fit &lt;- lmFit(elist, design, block=expdata$Donor, correlation=corfit$consensus.correlation)
results &lt;- topTable(fit, coef="TimepointT120",
                    n=nrow(elist), sort.by="none")
</code></pre>

<p>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.</p>

<p>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 <code>duplicateCorrelation</code> function. So is it really true that the analogous procedure (i.e. using <code>~Donor+Timepoint</code> 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.</p>

<p>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?</p>

<p>Gordon Smyth answered this on the Bioconductor mailing list. See his explanation <a href="https://mailman.stat.ethz.ch/pipermail/bioconductor/2014-February/057887.html">here</a>. Summarizing his answer, the two approaches are similar. `duplicateCorrelation` fits a random effect for the blocking variable. It is especially useful to use when the design is unbalanced.</p>

<p>Also, for the sake of posterity since this is an old post, using limma with voom is a suitable approach to analyzing RNA-seq data. See this <a href="http://genomebiology.com/content/14/9/R95">comparison of differential expression methods</a> and also the <a href="http://genomebiology.com/content/15/2/R29">voom publication</a>.</p>
