use the apeglm shrinkage estimator, and work out model details

[daler]

According to the extended shrinkage estimator section in the DESeq2 vignette, apeglm has better performance especially with lower-count genes where the default "normal" mode shrinks to aggressively.

However, as they describe there, this prevents us from using the contrast=c('factorame', 'treat', 'control') syntax for generating results. It's still possible to generate results, it just takes some rearrangement of the model to get it to work.

[mitraak]

Alternative shrinkage methods, apeglm and ashr:

The options for type are:

• apeglm is the adaptive t prior shrinkage estimator from the apeglm package (Zhu, Ibrahim, and Love 2018). As of version 1.28.0, it is the default estimator.
• ashr is the adaptive shrinkage estimator from the ashr package (Stephens 2016). Here DESeq2 uses the ashr option to fit a mixture of Normal distributions to form the prior, with method="shrinkage".
• normal is the the original DESeq2 shrinkage estimator, an adaptive Normal distribution as prior.

[mitraak]

Notes on apeglm:

• The main difference from normal is that we cannot use the convenient contrast parameter in our lfcShrink call. So, estimating different comparisons and fold-changes will need some refactoring.
  • dds$condition can be releveled
  • design(dds) can be replaced with a new design
  • the dispersions will not need to be estimated again, only the nbinomWaldTest step needs to be run again (not the full DESeq() call)

NOTE: Most of the following is based on DESeq2 docs linked above.

One-factor design w/ three groups

• In an experiment with three groups, A (reference), B, C, the resultsNames(dds) would be: intercept, groupB and groupC
  • groupB is the contrast: B vs A
  • groupC is the contrast: C vs A
• To get the contrast B vs C, we need to jump through a couple of hoops:

  dds$group <- relevel(dds$group, 'B')
  dds <- nbinomWaldTest(dds)

• Now the resultsNames(dds) will be: intercept, groupA and groupC. So, we can now get the contrast B vs C with

  resBvsC <- lfcShrink(dds, coef='groupC', type='apeglm')

Two factor-design: 3 groups, 2 conditions

• grp: 1,2,3
• cnd: A,B
• model: ~ grp + cnd + grp:cnd

Here the third term grp:cnd is the interaction term. There are three condition effects:

1. grp1, cndB vs cndA
2. grp2, cndB vs cndA
3. grp3, cndB vs cndA

and consequently, three interaction contrasts: 2 vs 1, 3 vs 1 and 3 vs 2.

The normal model matrix will have resultsNames(dds):

## [1]  grp2  grp3  cndB  grp2:cndB  grp3:cndB

• grp2:cndB == 2 vs 1
• similarly, grp3:cndB == 3 vs 1
• To obtain the 3 vs 2, we need to relevel the factor grp

  dds$grp <- relevel(dds$grp, '2')

• Now, after redoing the nbinomWaldTest step, resultsNames(dds) will include the following:

  ## [1] grp1  grp3  cndB  grp1:cndB  grp3:cndB

• Now, the coefficient grp3:cndB == 3 vs 2