<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>

CC BY-SA 4.0
2019-07-10T14:40:46.627Z
2019-07-10T14:40:46.627Z
James Hawley
<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>

CC BY-SA 4.0
2019-07-09T17:22:57.093Z
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?
James Hawley
CC BY-SA 4.0
2019-07-08T22:37:11.443Z
2019-07-08T22:37:11.443Z
swbarnes2
<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>

CC BY-SA 4.0
2022-01-27T15:01:18.997Z
what do mean by null hypothesis as normal ? Anyway, a good question +1.
Subhash C. Davar
<rna-seq><differential-expression><kallisto><p-values>
Strange p-value histogram for differential gene expression analysis
# 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 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

**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))



The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.


---

### 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?  

---

