Hi all,

once in a while I come upon an DESeq analysis which results in a large number of genes with low p.adjust values (~significant genes) but with a somewhat unsatisfying expression pattern: zero counts in nearly all samples. It is very likely I am missing a modelling parameter but I did search for solution and didn't find anything.

For example, in this re-analysis of a published dataset there are two distinct clouds with extreme fold-changes but at the low end of expression:

plotMA(res, ylim=c(min(res$log2FoldChange, na.rm = T), max(res$log2FoldChange, na.rm = T)))

These contain quite a few differentially expressed genes, and when I take the one with the lowest padj we can see that it has zero counts in all of the control samples (scrambled_KCl)

plotCounts(dds, gene=which.min(res$padj), intgroup="siRNA_stimulus")

gene_id <- rownames(res[which.min(res$padj), ])
samples <- subset(sample_info, siRNA_stimulus %in% c("scrambled_KCl", "NEAT1_KCl")) %>%
    dplyr::arrange(siRNA_stimulus) %>%
    .$sample_id

cnts[gene_id, samples] %>%
    knitr::kable()

On the one had it makes sense that DESeq calls this gene as differentially expressed. However, if we take the 20 genes with the lowest 20 padj values we see the pattern repeating and something even more interesting:

res_20 <- res[order(res$padj), ] %>% head(20)

res_20 %>%
    knitr::kable()


cnts[rownames(res_20), samples]  %>%
    knitr::kable()

Several genes have a single sample with counts > 0 for samples in those groups being compared. In this case one would assume that DESeq would not attribute a low p-value two those because there is little evidence for difference expression in those two groups. Showing the counts for all groups, we can see that this gene happens to be highly expressed in the samples of other groups:

samples_ordered <- sample_info %>%
    dplyr::arrange(siRNA_stimulus)

cnts["ENSG00000237260", samples_ordered$sample_id]

plotCounts(dds, gene="ENSG00000237260", intgroup="siRNA_stimulus")

My question, or doubt, is it expected that this gene would have a low p-value? If seems that while it sort of this, the model might consider that there is enough evidence in the the sample set as whole to attribute a low p-value, but IMO as a biologist, it doesn't pass the "eye test".

As a follow question, is there any way to model out these genes? In principle one could do pre-filtering for genes not expressed in less than x% samples of any condition, but for practical reasons I would rather not (merging results downstream), and I was under the assumption that DESeq would handle these cases elegantly.

Cheers!

What practical reason do you have for not wanting to filter these genes? It is very common to filter genes that have zero counts in most samples because you can get spurious results like this.