Answers.Body: <p>I found that using the Wald test instead of the Likelihood Ratio Test produced very different results with respect to the p-values.</p>
  
  <pre><code>sleuth_obj = sleuth_prep(
      metadata,
      extra_bootstrap_summary = TRUE,
      read_bootstrap_tpm = TRUE,
      num_cores = 1,
      target_mapping = transcripts_genes,
      aggregation_column = "ens_gene"
  )
  
  # same model of covariates
  sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")
  
  # using Wald test instead of LRT
  sleuth_obj = sleuth_wt(sleuth_obj, "cellLT1", "full")
  
  # extract results
  sleuth_table_genes &lt;- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = TRUE))
  sleuth_table_txs &lt;- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = FALSE))
  </code></pre>
  
  <p><a href="https://i.sstatic.net/8UExE.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8UExE.png" alt="Wald test gene-level"></a>
  <a href="https://i.sstatic.net/8D4je.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8D4je.png" alt="Wald test transcript-level"></a></p>
  
  <p>The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.</p>
  
  <p><a href="https://i.sstatic.net/bht5o.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/bht5o.jpg" alt="LRT vs Wald"></a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-10T14:40:46.627Z
Answers.LastActivityDate: 2019-07-10T14:40:46.627Z
Answers.UserDisplayName: James Hawley
Answers.Body: <p>Try the protocol here for generating a different (smaller) dispersion estimate</p>
  
  <p><a href="https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues" rel="nofollow noreferrer">https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-09T17:22:57.093Z
Answers.Text: I found this previously, but where DESeq2 has its test statistic be normal, Sleuth detects differentially expressed transcripts using a likelihood ratio test. `fdrtool` can create an empirical null model on normal data, but is the test statistic from Sleuth's test normal?
Answers.UserDisplayName: James Hawley
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-08T22:37:11.443Z
Answers.LastActivityDate: 2019-07-08T22:37:11.443Z
Answers.UserDisplayName: swbarnes2
Body: <p>I'm trying to perform pretty standard differential expression analysis using RNA-seq data.</p>
  
  <p>I've used <a href="https://pachterlab.github.io/kallisto/" rel="nofollow noreferrer">Kallisto</a> to perform RNA quantification and am using <a href="https://pachterlab.github.io/sleuth/" rel="nofollow noreferrer">Sleuth</a> to perform the differential expression analysis.</p>
  
  <p>I have 2 conditions that I'd like to compare, which is stored in the <code>metadata</code>:</p>
  
  <pre><code>&gt; metadata
      sample cell replicate batch                   path
  1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
  2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
  3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
  4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
  5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
  6:  Test-3 Test         3   3   /path/to/Batch3/Test-3
  </code></pre>
  
  <p>I'm attempting to test for differentially expressed genes between the <code>Base</code> and <code>Test</code> cells while controlling for <code>batch</code> (the sequencing batch). The R code to do this is adapted from <a href="https://pachterlab.github.io/sleuth_walkthroughs/pval_agg/analysis.html" rel="nofollow noreferrer">Sleuth's manual</a>:</p>
  
  <pre><code>sleuth_obj = sleuth_prep(
      metadata,
      extra_bootstrap_summary = TRUE,
      num_cores = 1,
      target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
      aggregation_column = "ens_gene"
  )
  
  # smooth raw Kallisto abundances
  sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")
  
  # fitting error measurements on reduced model
  sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")
  
  # compare models to calculate differentially expressed transcripts
  sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")
  
  # extract results
  sleuth_table = as.data.table(sleuth_results(
      sleuth_obj,
      'reduced:full',
      'lrt',
      show_all = FALSE
  ))
  sleuth_table_tx = as.data.table(sleuth_results(
      sleuth_obj,
      'reduced:full',
      'lrt',
      show_all = FALSE,
      pval_aggregate = FALSE
  ))
  </code></pre>
  
  <p>To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:</p>
  
  <p><a href="https://i.sstatic.net/FOKOx.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/FOKOx.png" alt="Gene-level p-value histogram"></a>
  <a href="https://i.sstatic.net/TAWCX.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/TAWCX.png" alt="Transcript-level p-value histogram"></a></p>
  
  <p>Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?</p>
  
  <h1>Edit following from @swbarnes2's answer:</h1>
  
  <p>I've used <code>fdrtool</code> to create an empirical null model of Sleuth's likelihood ratio test.</p>
  
  <pre><code>library(fdrtool)
  local_fdr_txs = fdrtool(so<span class="math-container">$tests$</span>lrt<span class="math-container">$`reduced:full`$</span>test_stat, statistic="normal")
  </code></pre>
  
  <p>Which gives the following plots:</p>
  
  <p><a href="https://i.sstatic.net/37rsY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/37rsY.png" alt="fdrtool"></a></p>
  
  <p>Plotting the distribution of these p-values produces:</p>
  
  <p><a href="https://i.sstatic.net/O3hjw.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/O3hjw.png" alt="pvalues after fdrtool"></a></p>
  
  <p>This still has a hump, but looks more anti-conservative than before.</p>
  
  <p>My remaining question is about whether I can use <code>fdrtool</code> to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-01-27T15:01:18.997Z
Comments.Text: what do mean by null hypothesis as normal ? Anyway, a good question +1.
Comments.UserDisplayName: Subhash C. Davar
Tags: <rna-seq><differential-expression><kallisto><p-values>
Title: Strange p-value histogram for differential gene expression analysis
ThreadText: # Strange p-value histogram for differential gene expression analysis
  
  (Asked by: username_0 on 2019-07-05)
  
  I'm trying to perform pretty standard differential expression analysis using RNA-seq data.
  
  I've used [Kallisto](https://pachterlab.github.io/kallisto/) to perform RNA quantification and am using [Sleuth](https://pachterlab.github.io/sleuth/) to perform the differential expression analysis.
  
  I have 2 conditions that I'd like to compare, which is stored in the `metadata`:
  
  ```
  > metadata
      sample cell replicate batch                   path
  1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
  2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
  3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
  4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
  5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
  6:  Test-3 Test         3   3   /path/to/Batch3/Test-3
  
  ```
  
  I'm attempting to test for differentially expressed genes between the `Base` and `Test` cells while controlling for `batch` (the sequencing batch). The R code to do this is adapted from [Sleuth's manual](https://pachterlab.github.io/sleuth_walkthroughs/pval_agg/analysis.html):
  
  ```
  sleuth_obj = sleuth_prep(
      metadata,
      extra_bootstrap_summary = TRUE,
      num_cores = 1,
      target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
      aggregation_column = "ens_gene"
  )
  
  # smooth raw Kallisto abundances
  sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")
  
  # fitting error measurements on reduced model
  sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")
  
  # compare models to calculate differentially expressed transcripts
  sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")
  
  # extract results
  sleuth_table = as.data.table(sleuth_results(
      sleuth_obj,
      'reduced:full',
      'lrt',
      show_all = FALSE
  ))
  sleuth_table_tx = as.data.table(sleuth_results(
      sleuth_obj,
      'reduced:full',
      'lrt',
      show_all = FALSE,
      pval_aggregate = FALSE
  ))
  
  ```
  
  To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:
  
  [![Gene-level p-value histogram](https://i.sstatic.net/FOKOx.png)](https://i.sstatic.net/FOKOx.png)
  [![Transcript-level p-value histogram](https://i.sstatic.net/TAWCX.png)](https://i.sstatic.net/TAWCX.png)
  
  Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?
  
  # Edit following from @username_2's answer:
  
  I've used `fdrtool` to create an empirical null model of Sleuth's likelihood ratio test.
  
  ```
  library(fdrtool)
  local_fdr_txs = fdrtool(so$tests$lrt$`reduced:full`$test_stat, statistic="normal")
  
  ```
  
  Which gives the following plots:
  
  [![fdrtool](https://i.sstatic.net/37rsY.png)](https://i.sstatic.net/37rsY.png)
  
  Plotting the distribution of these p-values produces:
  
  [![pvalues after fdrtool](https://i.sstatic.net/O3hjw.png)](https://i.sstatic.net/O3hjw.png)
  
  This still has a hump, but looks more anti-conservative than before.
  
  My remaining question is about whether I can use `fdrtool` to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test
  
  **Tags:** rna-seq, differential-expression, kallisto, p-values
  
  > **Comments:**
  > **username_1**: what do mean by null hypothesis as normal ? Anyway, a good question +1.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2019-07-10. Score: 1 | Accepted answer
  
  I found that using the Wald test instead of the Likelihood Ratio Test produced very different results with respect to the p-values.
  
  ```
  sleuth_obj = sleuth_prep(
      metadata,
      extra_bootstrap_summary = TRUE,
      read_bootstrap_tpm = TRUE,
      num_cores = 1,
      target_mapping = transcripts_genes,
      aggregation_column = "ens_gene"
  )
  
  # same model of covariates
  sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")
  
  # using Wald test instead of LRT
  sleuth_obj = sleuth_wt(sleuth_obj, "cellLT1", "full")
  
  # extract results
  sleuth_table_genes <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = TRUE))
  sleuth_table_txs <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = FALSE))
  
  ```
  
  [![Wald test gene-level](https://i.sstatic.net/8UExE.png)](https://i.sstatic.net/8UExE.png)
  [![Wald test transcript-level](https://i.sstatic.net/8D4je.png)](https://i.sstatic.net/8D4je.png)
  
  The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.
  
  [![LRT vs Wald](https://i.sstatic.net/bht5o.jpg)](https://i.sstatic.net/bht5o.jpg)
  
  
  ---
  
  ### Answer by username_2 on 2019-07-08. Score: 1
  
  Try the protocol here for generating a different (smaller) dispersion estimate
  
  <https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues>
  
  > **Comments:**
  > **username_0**: I found this previously, but where DESeq2 has its test statistic be normal, Sleuth detects differentially expressed transcripts using a likelihood ratio test. `fdrtool` can create an empirical null model on normal data, but is the test statistic from Sleuth's test normal?  
  
  ---
  
  