Question: Cutoffs of padj and fold change in DESeq2
2
gravatar for pmpierry
2.8 years ago by
pmpierry20
pmpierry20 wrote:

Hello everyone,

I have a doubt about the cutoffs in DESeq2 related to padj and fold change. I initially did all my analyses using the default parameters in "results" function. Afterwards, I decided to use cutoffs of padj<=0.05 e FC>=2.0. To do so, I used the results table from those default analysis and filtered based on these cutoffs I wanted.

Later, I read again DESeq2 documentation and saw the arguments "alpha" and "lfcThreshold" in "results" function. My question is if the way I did before (to filter the results table after default analysis) was correct or I need to do all my analysis again using the arguments in results function.

Thank you!

 

 

deseq2 • 3.1k views
ADD COMMENTlink modified 2.8 years ago • written 2.8 years ago by pmpierry20

Hello Michael,

I followed your suggestions to add alpha and lfcThreshold arguments to results():

res <- results(dds, alpha=0.05, lfcThreshold=1, altHypothesis="greaterAbs")

But, when I call "summary(res)", only alpha has changed, but LFC did not:

> summary(res)

out of 2154 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)  : 84, 3.9% 
LFC < 0 (down)  : 37, 1.7% 
outliers [1]     : 0, 0% 
low counts [2]   : 0, 0% 
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

I thought it would show "LFC > 1 (up)" and "LFC < -1 (down)" or "LFC < |1| (down).

Am I missing some other argument or setup in results()?

Thanks again for all attention and help!

ADD REPLYlink written 2.8 years ago by pmpierry20

I see what you mean. summary() is just a helper function to make a table and I don't have it coded to look at the lfcThreshold. I could try to change that, but anyway, you can see that you are in fact testing against LFC of 1 by looking at the MA plot.

And to get numbers of genes you can just do:

table(res$padj < 0.1 & res$log2FoldChange > 0)

table(res$padj < 0.1 & res$log2FoldChange < 0)

 

ADD REPLYlink written 2.8 years ago by Michael Love24k

Ok, I got it! Thank you again for your help!

ADD REPLYlink written 2.8 years ago by pmpierry20

By the way, I just fixed this in devel (so released in April 2018) so that it shows the LFC threshold value:

> dds <- makeExampleDESeqDataSet()
> dds <- DESeq(dds, quiet=TRUE)
> res <- results(dds, lfcThreshold=log2(1.5))
> summary(res)

out of 997 with nonzero total read count
adjusted p-value   < 0.1
LFC > 0.58 (up)    : 0, 0%
LFC < -0.58 (down) : 0, 0%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

> res <- lfcShrink(dds, coef=2, lfcThreshold=log2(1.5))
> summary(res)

out of 997 with nonzero total read count
adjusted p-value   < 0.1
LFC > 0.58 (up)    : 0, 0%
LFC < -0.58 (down) : 0, 0%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

> summary(results(dds))

out of 997 with nonzero total read count
adjusted p-value   < 0.1
LFC > 0 (up)       : 0, 0%
LFC < 0 (down)     : 0, 0%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?result
ADD REPLYlink written 15 months ago by Michael Love24k
Answer: Cutoffs of padj and fold change in DESeq2
0
gravatar for Michael Love
2.8 years ago by
Michael Love24k
United States
Michael Love24k wrote:

If you want to use a LFC filter you should use the lfcThreshold argument. This is preferable because it is actually testing for a LFC greater than the value you choose. It makes less sense IMO to test for a non-zero LFC and then filter a different value. For more on this, see the "Hypothesis tests with thresholds on effect size" section of the DESeq2 paper.

Regarding the adjusted p-value threshold you can choose what you want, but by supplying the threshold to the 'alpha' argument, you will increase power. The independent filtering routing within results() uses alpha to optimize the low count filtering. Choosing a different value for alpha than the one you use to filter is sub-optimal.

ADD COMMENTlink written 2.8 years ago by Michael Love24k

Thanks a lot for the quick answer!!

ADD REPLYlink written 2.8 years ago by pmpierry20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 268 users visited in the last hour