My design consists of one factor i.e. treatment with three levels

>dds <- phyloseq_to_deseq2(ps, ~treatment)
>levels(dds$treatment)
'Compost''Control''Forest''Groundnut'

>dds <- DESeq(dds, "LRT", reduced = ~ 1)
>resultsNames(dds)
'Intercept''treatment_Control_vs_Compost''treatment_Forest_vs_Compost''treatment_Groundnut_vs_Compost'

>lfcShrink(
  dds,
  res = results(
    dds,
    name = "treatment_Control_vs_Compost",
    test = "Wald"
  ),
  type = "apeglm",
  coef = "treatment_Control_vs_Compost"
) %>%
  data.frame() %>%
  filter(padj <= 0.05 & abs(log2FoldChange) >= 1)

I then ran lfcShrink for my co-efficient of interest i.e treatment_Control_vs_Compost and after filtering for padj <= 0.05 & abs(log2FoldChange) >= 1, I ended up with 169 features. However, when I repeated this after changing the reference level to Control, I get only 23 features for the same co-efficient

>dds$treatment <- relevel(dds$treatment, ref = "Control")
>dds <- nbinomWaldTest(dds)
>resultsNames(dds)
'Intercept''treatment_Compost_vs_Control''treatment_Forest_vs_Control''treatment_Groundnut_vs_Control'

>lfcShrink(
  dds,
  res = results(
    dds,
    name = "treatment_Compost_vs_Control",
    test = "Wald"
  ),
  type = "apeglm",
  coef = "treatment_Compost_vs_Control"
) %>%
  data.frame() %>%
  filter(padj <= 0.05 & abs(log2FoldChange) >= 1)

I expected more or less the same results, but this kind of surprised me

I don’t think the apeglm results should change much (eg only due to numerical differences).

You are showing that the adjusted pvalues are different, which can occur especially when the data are not a good fit to NB. Note that you’re not looking at adjusted pvalues from apeglm because it doesn’t make pvalues, only passes them through from results. If you want to produce svalues, you can either specify svalue=TRUE or use the threshold argument.

I don’t universally recommend DESeq2 for metagenomic/microbiome data as I’m not convinced the NB GLM per gene is always appropriate here.