# Strange p-value histogram for differential gene expression analysis

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

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

• 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 Jul 9 '19 at 17:22

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(
extra_bootstrap_summary = 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.