Question

DESeq2 shrunken log fold changes

0

Entering edit mode

dequattro.concetta ▴ 10

@dequattroconcetta-21510

Last seen 2.8 years ago

Italy

Hi! I am doing pairwise comparisons with DESeq2 (version 1.24.0). I would like to shrinkage the log2 Fold Change using the normal approach. I used the following code:

dds <- DESeqDataSet(se_sel, ~ condition)
dds <- DESeq(dds)
resNorm <- lfcShrink(dds, contrast=c("condition", cond1 , cond2 ), type="normal")

resNorm
log2 fold change (MLE): condition cond1 vs cond2
Wald test p-value: condition cond1 vs cond2
                         baseMean     log2FoldChange              lfcSE
                        <numeric>          <numeric>          <numeric>
ENSG00000015592  2252.24746949959   4.93239360176146  0.129430249903588
ENSG00000083857  38573.0261803919   3.96140286121604 0.0734986203310571
ENSG00000090339  2854.89111431782   4.08839129938656   0.10107762978087
ENSG00000171551  12450.6367006511   6.65073841803356  0.132315837227065
ENSG00000076356  8747.08707079265  -3.53308614922331 0.0955398226626009

I was wondering if it is correct to shinkage the LFC with the function lfcShrink or considering that I am using normal approach it is better to run DESeq function with betaPrior=TRUE and then use the function result with contrast to extract the results. In addition, I will considered in my analysis the file with shunkage LFC (resNorm). Is that correct?

Thank you for your help!

Concetta

deseq2 • 2.8k views

ADD COMMENT • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

Hi! After your answers I tested the apeglm and ashr LFC shrinkage methods. I used the following code:

# Analysis
dds <- DESeqDataSetFromMatrix(count=count_data, colData = sampleTable, design = ~ condition)
dds$condition <- relevel(dds$condition, ref = cond2)
dds <- DESeq(dds)
resApe <- lfcShrink(dds, coef=2, type="apeglm")
resApe_2 <- lfcShrink(dds, coef=2, type="apeglm", svalue = TRUE, returnList = TRUE)
############shrinkage ash
resAshr <- lfcShrink(dds, coef=2, type="ashr")
resAshr_2 <- lfcShrink(dds, coef=2, type="ashr", svalue = TRUE, returnList = TRUE)

Setting the option returnList = TRUE I did not understand how to use the output of apeglm or ashr. Can I do just a filtering based on the s-value? In addition, when I draw the MA plot I got this message:

plotMA(resApe_1, ylim=c(-5,5))
thresholding s-values on alpha=0.005 to color points

What does it means? It indicates that it take into account only genes with ajusted p-value < 0.005.

Concerning the LFC shrinkage, I have a couple of questions about the result:

How can I decide which is the best shrinkage method for my data?
How can I exploit the s-value. Should I use it like the adjusted p-value? For example considering only genes with s-value < 0.05.

In addition I would like to ask which is the best approch to filter differentially expressed genes exploiting s-value and/or adjusted p-value:

Should I consider only the p-value adjusted and a LFC (after shrinkage) threshold to discriminate genes up and down regulated?
Considering that the lfcshinkage function provide a s-value, should I consider only the s-value instead of p-value. In this case which s-value value should I set?
Otherwise the best approach is to consider an adjusted p-value < 0.05 and a s-value < 0.05 and then take into account only genes with a LFC >1 or LFC < -1?

Thank you for your help!

Concetta

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

"Can I do just a filtering based on the s-value?"

Yes, this will give a set of genes satisfying an aggregate FSR (see references). You can choose the threshold yourself (plotMA has an argument alpha). See the apeglm vignette for some discussion, but I'd set it a bit lower than the typical 0.05, perhaps 0.01 or 0.005.

Regarding the methods, we recommend apeglm or ashr, given the results from the apeglm publication. We don't recommend type="normal" anymore.

Regarding adjusted p-value or s-value, you're better off picking one or the other, don't combine these two in your results. These are different frameworks. There are good arguments for using s-values, but I admit that most papers continue to use adjusted p-values.

ADD REPLY • link 4.6 years ago Michael Love 43k

0

Entering edit mode

OK thank you so much for your reply. I would like to ask just a clarification about DESeq2. If I understood well the DESeq2 do not perform the shrinkage of log2 Fold Change and so in the results table I have log2 fold change (MLE): dex trt vs untrt that indicates that the LFC is not shrinked. If I want to shrink the log2 fold change I have to use the function lfcShrink. Subsequently, for the downstream analysis if I want to filter genes based on the effect of the treatment it is better to consider the shrinked fold change. Is that correct?

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

You can filter on the "un-shrunken" or MLE LFC, just use lfcThreshold. Search the vignette for this term.

ADD REPLY • link 4.6 years ago Michael Love 43k

0

Entering edit mode

I tested lfcThreshold on my data in the us-shrunken" MLE LFC. However, why did I obtain a different results from the result that I obtained with the following command?

print(length(rownames(res[res$padj < 0.05 & !is.na(res$padj) & abs(res$log2FoldChange) > 1,]))

In addition, do you suggest to work on us-shrunken" MLE LFC or on shrunken apeglm LFC?

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

That is answered in the 2014 paper, where we describe what the methods are doing.

If you just want to do what most DESeq2 users are doing in papers then it’s adjusted pvalues and maybe lfcThreshold. But it’s up to you.

ADD REPLY • link 4.6 years ago Michael Love 43k

score 0 · Answer 1 · 2020-04-20

0

Entering edit mode

Michael Love 43k

@mikelove

Last seen 10 hours ago

United States

Here is what lfcShrink prints in the latest version of DESeq2:

> packageVersion("DESeq2")
[1] ‘1.26.0’
> dds <- makeExampleDESeqDataSet()
> dds <- DESeq(dds, quiet=TRUE)
> res <- lfcShrink(dds, coef=2, type="normal")
using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).

Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
Reference: https://doi.org/10.1093/bioinformatics/bty895

betaPrior=TRUE is quasi-deprecated. It's non-default, but still allowed because I didn't want to break anyone's downstream scripts or pipelines unnecessarily. lfcShrink with apeglm or ashr is now recommended.

ADD COMMENT • link 4.6 years ago Michael Love 43k

0

Entering edit mode

Thank you for your reply. Is it right to use contrast instead of coef in the lfcShrink function? Does it affect the result?

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

Check out the vignette — we have a lot of information about this in the extended section on shrinkage estimators.

Then if after reading over that section you still have a question feel free to post here again.

ADD REPLY • link 4.6 years ago Michael Love 43k

0

Entering edit mode

Thank you. I have read the vignette. Now it is clear, however I tested the example code in the manual and I observed that genes log2FC is different in res.shr and res.shr_2 (obtained with contrast). I did not understand why the LFC is different because if I understood well the two approaches, contrast and coef, should be the same.

set.seed(1)
dds <- makeExampleDESeqDataSet(n=500,betaSD=1)
dds <- DESeq(dds)
res <- results(dds)
# these are the coefficients from the model
# we can specify them using 'coef' by name or number below
resultsNames(dds)
res.shr <- lfcShrink(dds=dds, coef=2)
res.shr_2 <- lfcShrink(dds=dds, contrast=c("condition","B","A"))

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

In ?lfcShrink, we explain that, for type="normal" across these two variations, the software will use different design matrices, and so you will get different results, "coef will make use of standard model matrices, while contrast will make use of expanded model matrices".

Nevertheless, type="normal" is no longer recommended, as demonstrated in the apeglm paper.

ADD REPLY • link 4.6 years ago Michael Love 43k

0

Entering edit mode

OK thank you. However which approach do you suggest to use, coef or contrast? In addition in my DE analysis I have two conditions and I compared the two condition, building the dds and extracting the result using contrast function as follow.

dds <- DESeqDataSet(se, ~ condition)
dds <- DESeq(dds)
resNorm <- lfcShrink(dds, contrast=c("condition", cond1 , cond2 ), type="normal")

Is it OK for two conditions to use contrast function or it is better to specify the reference level? Does the two approaches change the results? However I am using type="normal", because I performed a prevoius analysis with earlier version of DESeq2.

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10

0

Entering edit mode

which approach do you suggest to use

I suggest not to use type="normal", and don't really have opinions on the options within the unrecommended type. We've kept in the older method for backwards compatibility, but not as a recommended pipeline. Are you trying to exactly recreate your previous analysis, or perform a best practice new analysis? Is there a reason that you have to stick with the old results/methods?

I don't think re-leveling will affect contrast and type="normal".

ADD REPLY • link 4.6 years ago Michael Love 43k

0

Entering edit mode

OK thank you. I want also to use lfcShrink with apeglm.
Now I am doing the analysis also with type="normal" because it has been performed a DE analysis on other samples with a version of DESeq2 < 1.16 and I want to compare the results.

ADD REPLY • link 4.6 years ago dequattro.concetta ▴ 10