Answers.Body: <p>Something like this <a href="https://web.stanford.edu/class/bios221/book/Chap-Testing.html#the-false-discovery-rate" rel="nofollow noreferrer">book chapter</a> might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.</p>
  <p>Basically, the <code>p.adjust(..method = &quot;BH&quot;)</code> is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.</p>
  <p>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.</p>
  <p>However, I have to note that in your code, you are not using <code>fdrtool</code> to estimate the fdr as explained above, because:</p>
  <pre><code>FDR.ddsRes &lt;- fdrtool(ddsRes$stat, statistic= &quot;normal&quot;, plot = FALSE)
  </code></pre>
  <p>With this, the raw p-values are stored in <code>FDR.ddsRes$pval</code> and the q value, which is the corrected p-value, is stored in <code>FDR.ddsRes$qval</code>. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:</p>
  <pre><code>p.adjust(FDR.ddsRes$pval, method = &quot;BH&quot;)
  </code></pre>
  <p>You should apply a two tailed test. From the wald stats, it should be:</p>
  <pre><code>derived_p = 2*pnorm(-abs(ddsRes$stat))
  </code></pre>
  <p>In your code, you are running it wrongly by passing the test statistic into <code>fdrtools</code>, we can compare the raw p-values:</p>
  <pre><code>head(data.frame(stat = ddsRes<span class="math-container">$stat,derived = derived_p,
  DESeq2=ddsRes$</span>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
  </code></pre>
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-24T09:09:47.190Z
Answers.LastActivityDate: 2021-05-24T09:09:47.190Z
Answers.UserDisplayName: StupidWolf
Body: <p>I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:</p>
  <pre><code>FDR.ddsRes &lt;- fdrtool(ddsRes$stat, statistic= &quot;normal&quot;, plot = T)
  ddsRes[,&quot;padj&quot;]  &lt;- p.adjust(FDR.ddsRes$pval, method = &quot;BH&quot;)
  </code></pre>
  <p>Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?</p>
  <p>Thanks!!</p>
  
Tags: <r><rna-seq><deseq2><rna><p-values>
Title: How does FDRtool work?
ThreadText: # How does FDRtool work?
  
  (Asked by: username_0 on 2021-05-22)
  
  I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:
  
  ```
  FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = T)
  ddsRes[,"padj"]  <- p.adjust(FDR.ddsRes$pval, method = "BH")
  
  ```
  
  Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?
  
  Thanks!!
  
  **Tags:** r, rna-seq, deseq2, rna, p-values
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-05-24. Score: 1
  
  Something like this [book chapter](https://web.stanford.edu/class/bios221/book/Chap-Testing.html#the-false-discovery-rate) 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 = "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,
  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.
  
  
  ---
  
  