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:

&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


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