I'm trying to perform pretty standard differential expression analysis using RNA-seq data.
I have 2 conditions that I'd like to compare, which is stored in the
> 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
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:
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