hello everybody,
I noticed some inconsistency in the DESeq2 results when calling results() under the default and the alternative hypothesis option.
Some genes are called signifcant above a certain threshold (e.g. LFC>0.5) but not significant when no threshold is applied (i.e. just differentially expressed no matter the strength).
The issue is caused by independent filtering, which in my case filters out lowly-expressed genes under the default (i.e. sets their padj to NA) but does not so under the alternative hypothesis.
Is this
(a) just some negligible consequence of independent filtering and perhaps the tradeoff for having a greater number of DEG under the default
(b) possibly a kind of bug and filtering should be the same for both
or
(c) am I doing/understanding something wrong?
pls let me know if you require additional information
btw: I like to add the stereotypical "I am a biologist and don't know much about statistics/informatics"-excuse
bye
I noticed some inconsistency in the DESeq2 results when calling results() under the default and the alternative hypothesis option.
Some genes are called signifcant above a certain threshold (e.g. LFC>0.5) but not significant when no threshold is applied (i.e. just differentially expressed no matter the strength).
The issue is caused by independent filtering, which in my case filters out lowly-expressed genes under the default (i.e. sets their padj to NA) but does not so under the alternative hypothesis.
Is this
(a) just some negligible consequence of independent filtering and perhaps the tradeoff for having a greater number of DEG under the default
(b) possibly a kind of bug and filtering should be the same for both
or
(c) am I doing/understanding something wrong?
pls let me know if you require additional information
btw: I like to add the stereotypical "I am a biologist and don't know much about statistics/informatics"-excuse
bye
Code:
[COLOR="Blue"]> #without independent filtering > #standard method > res <- results(dds, contrast = c("condition", "C", "A"), independentFiltering = FALSE) > summary(res, alpha=0.01)[/COLOR] out of 13821 with nonzero total read count adjusted p-value < 0.01 LFC > 0 (up) : 547, 4% LFC < 0 (down) : 896, 6.5% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 0) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results [COLOR="blue"]> res <- na.omit(res) > res <- res[res[,6] < 0.01,] > #alternative hypothesis > resGA <- results(dds, contrast = c("condition", "C", "A"), independentFiltering = FALSE, + lfcThreshold=0.5, altHypothesis="greaterAbs") > summary(resGA, alpha=0.01)[/COLOR] out of 13821 with nonzero total read count adjusted p-value < 0.01 LFC > 0 (up) : 40, 0.29% LFC < 0 (down) : 9, 0.065% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 0) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results [COLOR="blue"]> resGA <- na.omit(resGA) > resGA <- resGA[resGA[,6] < 0.01,] > #is resGA subset of res? > table(is.element(rownames(resGA), rownames(res)))[/COLOR] TRUE 49 [COLOR="blue"]> > #with independent filtering > #standard method > res <- results(dds, contrast = c("condition", "C", "A"), independentFiltering = TRUE) > summary(res, alpha=0.01)[/COLOR] out of 13821 with nonzero total read count adjusted p-value < 0.01 LFC > 0 (up) : 549, 4% LFC < 0 (down) : 911, 6.6% outliers [1] : 0, 0% low counts [2] : 1382, 10% (mean count < 15.2) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results [COLOR="blue"]> res <- na.omit(res) > res <- res[res[,6] < 0.01,] > #alternative hypothesis > resGA <- results(dds, contrast = c("condition", "C", "A"), independentFiltering = TRUE, + lfcThreshold=0.5, altHypothesis="greaterAbs") > summary(resGA, alpha=0.01) [/COLOR] out of 13821 with nonzero total read count adjusted p-value < 0.01 LFC > 0 (up) : 40, 0.29% LFC < 0 (down) : 9, 0.065% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 0) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results [COLOR="blue"]> resGA <- na.omit(resGA) > resGA <- resGA[resGA[,6] < 0.01,] > #is resGA subset of res? > table(is.element(rownames(resGA), rownames(res)))[/COLOR] FALSE TRUE 2 47
Comment