How does FDRtool work?
Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.
Basically, the p.adjust(..method = &quot;BH&quot;) is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.
When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.
However, I have to note that in your code, you are not using fdrtool to estimate the fdr as explained above, because:
FDR.ddsRes &lt;- fdrtool(ddsRes$stat, statistic= &quot;normal&quot;, plot = FALSE)

With this, the raw p-values are stored in FDR.ddsRes$pval and the q value, which is the corrected p-value, is stored in FDR.ddsRes$qval. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:
p.adjust(FDR.ddsRes$pval, method = &quot;BH&quot;)

You should apply a two tailed test. From the wald stats, it should be:
derived_p = 2*pnorm(-abs(ddsRes$stat))

In your code, you are running it wrongly by passing the test statistic into fdrtools, we can compare the raw p-values:
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.