I am performing a DGE analysis on plant tissues under different stresses with or without a separate treatment: moderate (D) and severe drought (S), with or without a separate treatment (T). Each combination has an unstressed/untreated control, e.g. MUU (moderate drought control, unstressed and untreated), MDU (moderate drought stressed, no treatment), MDT etc., for both levels of stress and different tissues. The analysis is simple in the sense that we are only analysing stressed/treated samples versus their control in a particular tissue: MDUL vs. MUUL, MDTL vs. MUUL in leaf, for example (each sample has three biological replicates).

In order to later use lfcShrink() with apeglm it is necessary to use the coefficient name ("group_MDTL_vs_MUUL") as opposed to a contrast ("group","MDTL","MUUL"). Thus, I releveled the data for each comparison as discussed in the vignette and elsewhere:

dds=DESeqDataSetFromTximport(txi,colData=sampletable,design=~group)
rdds=DESeq(dds)
rdds$group=relevel(rdds$group,ref="MUUL")
rdds=nbinomWaldTest(rdds)
resultsNames(rdds)
[1] "Intercept" ... "group_MDTL_vs_MUUL"
res=results(rdds,name="group_MDTL_vs_MUUL",alpha=0.05,lfcThreshold=1,altHypothesis="greaterAbs")
res_shrink=lfcShrink(dds=rdds,coef="group_MDTL_vs_MUUL",res=res,type="apeglm")

When I was analysing the output for the comparisons I noticed several instances where a gene was called DE despite intuitively seeming like the counts/number of samples with non-zero counts was too low. This included several instances where a gene had literally zero counts in both of the compared samples (bearing in mind the 0.5 pseudocount from plotCounts):

$ res["Solyc00g500050.1.1",]
log2 fold change (MLE): group MDTL vs MUUL 
Wald test p-value: group MDTL vs MUUL 
DataFrame with 1 row and 6 columns
                    baseMean log2FoldChange     lfcSE      stat      pvalue
                   <numeric>      <numeric> <numeric> <numeric>   <numeric>
Solyc00g500050.1.1   31.6121       -20.7792   5.10032  -3.87803 0.000105307
                         padj
                    <numeric>
Solyc00g500050.1.1 0.00381324

$ plotCounts(rdds,"Solyc00g500050.1.1",intgroup="group",return=TRUE) %>% filter(group == "MDTL" | group == "MUUL")
      count group
MDTL1   0.5  MDTL
MDTL2   0.5  MDTL
MDTL3   0.5  MDTL
MUUL1   0.5  MUUL
MUUL2   0.5  MUUL
MUUL3   0.5  MUUL

Using contrast rather than name seemed to prevent these genes being given significant p-values, as did using lfcShrink and filtering genes by s-value instead. I suspect the issue may be related to the highly variable expression of genes across the tissues and conditions, since a gene may be highly expressed in some tissue/condition and zero in many others? This also makes it difficult to pre-filter such genes.

1) Are there any major errors in the pipeline above?

2) Is there any way to avoid identifying such genes, other than a post-hoc filter?

Thanks!

The LFC is unreliable when both groups are 0, and likewise the SE (both are undefined). The solution is to either use contrast or use lfcShrink.