Hello Michael,

So far I have used DESeq 2 with betaPrior = T in the DESeq() function but now I'm testing it with betaPrior = False and using lfcShrink either with the coef or contrast parameter (DESeq2 version 1.18.0). I notice a couple of things which I think need clarifcation

I used a fairly common experimental design with a factor that has 3 levels and DM is initial reference level.

coldata$Cond <- factor(coldata$Cond, levels = c('DM', 'DP', 'Co'))
ddsdm <- DESeqDataSetFromMatrix(countData = counttab, colData = coldata, design = ~ Cond)
ddsdm <- DESeq(ddsdm, parallel=T, betaPrior = F)

With that I ran a couple of comparison and I'm pasting the command calls and the fold changes and was is shown if one does a head on the result obj

1a) Comparison DM vs DP with DP as denumerator

res_dm_dmVdp <- results(ddsdm, contrast = c('Cond', 'DM', 'DP'), alpha = alpha)

Cond DM vs DP and log2 fold change (MLE): Cond DM vs DP
LFC > 0 (up)     : 576
LFC < 0 (down)   : 401

1b) Applying lfcShrink to the result

res_dm_dmVdp_shrink <- lfcShrink(ddsdm, coef = 2, res = res_dm_dmVdp)

Cond DM vs DP and log2 fold change (MAP): Cond DP vs DM
LFC > 0 (up)     : 401
LFC < 0 (down)   : 576

As one can see the genes that were down regulated are now up and vice versa. p-value and padj values stay the same. I went through the documentation (https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html) and tried to find a section which could explain this and didn't find anything. I did however find a post on bioconductor which explains it (DESeq2 lfcShrink() usage of coef vs. usage of contrast) and you explain that this is an intendedbehavior of the function. I think it would be best to describe this better in the manual because many people will just use that program and lfcShrink and may not check it as good as they should and report the wrong results.

2a) This will even more confusion when one has more than two levels in a factor. For example, here I compare Co to DP

res_dm_coVdp <- results(ddsdm, contrast = c('Cond', 'Co', 'DP'), alpha = alpha)

Cond Co vs DP and log2 fold change (MLE): Cond Co vs DP
LFC > 0 (up)     : 100
LFC < 0 (down)   : 371

2b) adding lfcShrink

res_dm_coVdp_shrink <- lfcShrink(ddsdm, coef = 2, res = res_dm_coVdp)

Cond Co vs DP and log2 fold change (MAP): Cond DP vs DM
LFC > 0 (up)     : 222
LFC < 0 (down)   : 249

and I get the foldchanges for DP vs DM but the pvalues and padj values for Co vs DP

2c) One can make this correct by using contrast instead of coef

res_dm_coVdp_shrink_contrast <- lfcShrink(ddsdm, res = res_dm_coVdp, contrast = c('Cond', 'Co', 'DP'))

Cond Co vs DP and log2 fold change (MLE): Cond Co vs DP
LFC > 0 (up)     : 103
LFC < 0 (down)   : 368

and we are back to Co vs DP for fold changes, pvalue and padj. However, 3 genes went from down-regulated to up-regulated. I checked these genes and the highest basemean was 107. So it's lowly expressed but what if the gene has a high expression and it's a good candidate for experimental validation and suddenly the fold change is flipped. I think that's quite hard to explain to biologist. Could you explain why this can happen? I also noticed the the difference in fold change between coef and contrast become higher the lower the expression of a gene is.

Again the different behavior of coef and contrast has already been explained (DESeq2 lfcShrink() usage of coef vs. usage of contrast) but I think it would be good to have it in the documentation too. And maybe show the users a design that relates to DESeq2 version before 1.16.1

3) Lastly my question is how do I use lfcShrink when I have multi level factors and I want to compare Co vs DP, DM vs DP or Co vs DM.  lfcShrink with coef will always give me DP vs DM because that's how it has been set up in the beginning. The only way I can think of is to do a relevel every time

ddsdm_reset <- ddsdm
ddsdm_reset$Cond <- relevel(ddsdm_reset$Cond, ref = "DP")
ddsdm_reset <- DESeq(ddsdm_reset, parallel=T, betaPrior = F)

res_dm_dmVdp_reset <- results(ddsdm_reset, contrast = c('Cond', 'DM', 'DP'), alpha = alpha)

Cond DM vs DP and log2 fold change (MLE): Cond DM vs DP
LFC > 0 (up)     : 576
LFC < 0 (down)   : 401

res_dm_dmVdp_reset_shrink <- lfcShrink(ddsdm_reset, coef = 2, res = res_dm_dmVdp_reset)

Cond DM vs DP and log2 fold change (MLE): Cond DM vs DP
LFC > 0 (up)     : 576
LFC < 0 (down)   : 401

But how do you set it to have Co vs DP for the lfcShrink with the relevel?

You wrote in the two threads (DESeq2 lfcShrink() usage of coef vs. usage of contrast and https://support.bioconductor.org/p/95695/) that the future way is DESeq => results => lfcShrink and coef because that is future proof. So how do I use it for the above shown example.

4) I noticed a result object that has been run with IHW looses the weight columns after passing it to lfcShrink

Thanks for taking the time to answer it and if the above is already explained somewhere, then please can you provide a link.

We are already discussing this at our site, especially now that results are changed and one can go different ways (betaPrior set to True or False but then lfcShrink with coef or contrast) and it would good to have a way to guarantee a backwards comparability to older version (because somehow this has to be explained to the biologists) and have the new improvements too.

Best

Mathias 

Hi Mat,

Thanks for the feedback, your comments have been valuable in the past, and I'm sure we can update the documentation more. I've been spending most of my free time to work on the new methods, and on the primary documentation in the man pages, and now that they are in place, I will have time to add to the DESeq2 vignette.

Regarding, "I simply did not make the connection from the text, maybe due to lack of technical vocabulary. It would be good if it's explained in the html on bioconductor.",

Do you mean that I should put more text about 'coef' and 'contrast' in the section on lfcShrink in the vignette? This is probably a good idea.

If you want to compare all levels of a factor, I'd recommend sticking with the use of 'contrast', when calling lfcShrink(). This will be supported going forward. (Also supported going forward is the use of betaPrior=TRUE, which provides 100% backwards compatibility, avoiding the use of a new function, so that users who don't want to incorporate any additional functionality would not have to change their scripts except for a single argument. I think there hasn't been any change which removes functionality that users had previously.)

Another recommendation from the above, is that, it was possible for a user to overwrite a results table with different LFCs than the p-value, by first using one contrast with results(), then using another with lfcShrink() and providing the previous results table. Two notes about this: the contrast associated with LFC, p-values, etc is always recorded in the mcols(res), and also printed when the results table is printed to console. So at least it can be spotted by eye. But more importantly, in v1.18, the recommended code doesn't provide res anymore, so it is created inside lfcShrink() using the 'contrast' or 'coef'. This is more foolproof, so that users cannot possibly specify different contrasts for the results table and the LFC.