Due to the increase in functionality in edgeR, I've become slightly confused about the correct order of calling the different functions, since there seems to be multiple ways of doing the same thing. Specifically, there are currently 3-4 ways of estimating dispersion:

1) 3-step method:

• estimateGLMCommonDisp
• estimateGLMTrendedDisp
• estimateGLMTagwiseDisp

2) All-in-one method:

• estimateDisp

3) Robust method

• estimateGLMRobustDisp

4) QL

• estimateDisp
• glmQLFit

The estimateDisp function help file says it's similar to calling estimateGLMCommonDisp > estimateGLMTrendedDisp > estimateGLMTagwiseDisp, but also that it gives slightly different results. What's the difference between these two approaches, and is estimateDisp always the preferred?

The estimateDisp also has an argument called robust, which is suggested to be set to TRUE in the edgeR user guide. Is running estimateDisp(robust=TRUE) similar to running estimateGLMRobustDisp? 

To make the confusion total, there's the separate approach of using the QL approach instead. The manual suggests using estimateDisp(robust=TRUE), even glmQLFit re-estimates the dispersions. Again, glmQLFit also has an argument called robust. In addition, it has an argument called abundance.trend (which I have not seen mentioned in any manuals), should this ever be set to TRUE?

What's considered a good rule of thumb for running the different implementations? Something like below:

edgeR:

1. estimateDisp(robust=TRUE)
2. glmFit
3. glmLRT

edgeR robust:

1. estimateGLMRobustDisp
2. glmFit
3. glmLRT

edgeR QL:

1. estimateDisp(robust=TRUE)
2. glmQLFit(robust=TRUE)
3. glmQLFTest

Any advise is much appreciated!

Your call sequences are correct for all three approaches. A couple of comments, though:

1. As for why estimateDisp is different from the common/trended/tagwise trio of methods, see C: edger, trended or common dispersion.
2. No, estimateGLMRobustDisp is distinct from robust=TRUE in estimateDisp. The former uses observation weights to reduce the impact of outlier counts within each gene. The latter doesn't protect each gene from outlier observations, but instead protects the empirical Bayes shrinkage from genes with outlier dispersions (that are caused by outlier counts). This difference in behaviour has some practical consequences. For example, a gene that is DE but has a couple of outlier counts in one group will (hopefully) still be detected as DE with estimateGLMRobustDisp. This is because those outlier counts should be downweighted, thus avoiding increases to the dispersion estimate. In contrast, with estimateDisp, the dispersion gets inflated by the outliers so there won't be enough power to detect DE - however, any deleterious effect of the inflated dispersion on EB shrinkage of all other genes is prevented with robust=TRUE. I use estimateDisp as it's faster (no need to iteratively compute weights), but there can be some benefit from robustifying against outlier observations - read the NAR paper.
3. glmQLFit estimates QL dispersions, it doesn't re-estimate the NB dispersions; keep in mind that they are two separate sets of values. Technically, you don't need to set robust=TRUE in estimateDisp because that only affects the tagwise NB dispersions that are never used in the QL framework - only the trended NB dispersions are used. That said, robustifying doesn't do any harm so you might as well do it if you plan to use the tagwise NB dispersions for diagnostics later, e.g., with plotBCV. Note that the EB shrinkage is repeated within glmQLFit using the QL dispersions. Here, robust=TRUE in glmQLFit is analogous to that in estimateDisp, as it protects the EB shrinkage from outlier dispersions (QL, not NB). Finally, abundance.trend=TRUE is already the default, and just ensures that EB shrinkage of the QL dispersions is performed towards a mean-dependent trend based on the QL dispersions. This trend is distinct from that generated by estimateDisp, which concerns the NB dispersions.

So, all in all, there's a lot of options that can be very confusing, but all of them work reasonably well so it's not too critical what you pick. I would recommend using the QL pipeline because it provides more accurate error control as well as robustness to outliers. (The fact that I'm first author on the associated report is, of course, just a coincidence.)

Fair question. We do offer a lot of alternative pipelines and it could be confusing.

Aaron has already covered most things, but to understand what robust=TRUE does in estimateDisp(), glmQLFit() and eBayes() see our recently accepted manuscript:

  http://arxiv.org/abs/1602.08678

The option has the same meaning in all three functions.

Here is quick overview of the different robustness concepts:

1. One or more samples are outliers: use arrayWeights() or voomWithQualityWeights() in limma.
2. One or more genes are more variable than the others: use robust=TRUE in the Bayes step in either edgeR or limma.
3. There are outlier observations, but not associated with any particular sample or any particular gene: use lmFit(,method="robust") in limma or estimateGLMRobustDisp() in edgeR.

The third option is what statisticians usually think of as robustness, but the first two options are more efficient because they involve global models. Option 1 assumes that outlier observations tend to occur for particular samples. Option 2 assumes that outlier observations tend to be associated with particular genes or with a group of genes. This is more common than you might think.