text: user240: I'm trying to perform pretty standard differential expression analysis using RNA-seq data.
  I've used Kallisto to perform RNA quantification and am using 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:
  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:
  
  
  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 @swbarnes2'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:
  
  Plotting the distribution of these p-values produces:
  
  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
  
  user2388: 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
  
  user240: 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))
  
  
  
  The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.
  
  
original_text: 240: <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.stack.imgur.com/FOKOx.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/FOKOx.png" alt="Gene-level p-value histogram"></a>
  <a href="https://i.stack.imgur.com/TAWCX.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/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.stack.imgur.com/37rsY.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/37rsY.png" alt="fdrtool"></a></p>
  
  <p>Plotting the distribution of these p-values produces:</p>
  
  <p><a href="https://i.stack.imgur.com/O3hjw.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/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>
   2388: <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>
   240: <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.stack.imgur.com/8UExE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8UExE.png" alt="Wald test gene-level"></a>
  <a href="https://i.stack.imgur.com/8D4je.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/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.stack.imgur.com/bht5o.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bht5o.jpg" alt="LRT vs Wald"></a></p>
  