DEseq2: can filtering low count genes reduce the number of DE genes?

[JonB]

Hi,

I tried to filter out genes with a count less than 2.0, and it seems that the lowest padj value for these genes is 1.350968e-44, and that removing these genes reduces the number of DE genes. Is this possible or have I done something wrong?

> filterTreshold = 2.0
> keep = rowMeans(counts(ddsHTSeq, normalized=TRUE)) > filterTreshold
> table(keep)
keep
FALSE  TRUE 
 8007 13673 

> min(res$padj[!keep], na.rm=TRUE)
[1] 1.350968e-44
> table(p.adjust(res$pvalue, method="BH") < .1)

FALSE  TRUE 
17870  2610 
> table(p.adjust(res$pvalue[keep], method="BH") < .1)

FALSE  TRUE 
11285  1645

[Michael Love]

hi Jon,

Please read over the DESeq2 vignette, where we describe that independent filtering is automatically performed within the results() function. A threshold on the row means of normalized counts is chosen which maximizes the number of adjusted p-values less than a given threshold. This procedure guarantees that the number of DE genes (as defined by adjusted p-value) will stay equal or increase due to filtering.

[Michael Love]

In other words, yes you can reduce the number of DE genes if you choose an arbitrary threshold (imagine setting the threshold to Inf).

This is why we use the genefilter functions internally to choose the optimal threshold.