1
$\begingroup$

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:

Gene-level p-value histogram Transcript-level p-value histogram

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:

fdrtool

Plotting the distribution of these p-values produces:

pvalues after fdrtool

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – James Hawley Jul 9 at 17:22
1
$\begingroup$

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

Wald test gene-level Wald test transcript-level

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

LRT vs Wald

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.