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.
p.adjust(..method = "BH") 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 <- fdrtool(ddsRes$stat, statistic= "normal", 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 = "BH")
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,
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.