Volcano plot with aggregated p values from Sleuth

Question

[pv342]

Hi,

I'm using Sleuth to perform a differential expression analysis. If I perform p value aggregation, is it appropriate to calculate log2fc from kallisto_output() and then plot the p values versus these log2fc values?

Thank you.

Answer 1

[dsull]

I see nothing wrong with that mathematically (gene-level log fold change and [aggregated] p-value are two different things -- and it's perfectly valid to plot them against one another or to generate DE gene lists from both). The only issue is that it's not convenient to do -- this will hopefully be fixed soon in a future release of sleuth (see: https://github.com/pachterlab/sleuth/issues/233 )

[pv342]

Do you have any recommendations for comparing gene-level p-values and transcript-level fold changes? Is it okay to aggregate (by sum) transcript abundances (tpm) by gene_id and calculate gene-level fold changes from these? Otherwise I have multiple fold changes for each gene p-value. I would appreciate any thoughts. Thank you!

[dsull]

To get gene-level fold changes, you basically simply set gene_mode=TRUE in sleuth_prep and sleuth_results.

No, don't sum TPMs and calculate gene-level fold changes yourself; to get fold changes, you need to compare between samples which necessitates proper normalization (which sleuth will perform).

[pv342]

I'm running Sleuth (0.30.0 in R4.0.4) in gene mode as you recommend. However, I do not see any fold change related outputs. Do I get this using some other function?

so <- sleuth_prep(design_table, extra_bootstrap_summary = TRUE, target_mapping = t2g, aggregation_column = 'ens_gene', gene_mode = TRUE)

so <- sleuth_fit(so, ~condition, 'full')
so <- sleuth_fit(so, ~1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')

sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE, gene_mode = TRUE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)

The dataframe sleuth_significant has the following columns:

• target_id ext_gene
• pval
• qval
• test_stat
• rss
• degrees_free
• mean_obs
• var_obs
• tech_var
• sigma_sq
• smooth_sigma_sq
• final_sigma_sq

so$gene_mode returns TRUE

[dsull]

Only the wald test ('wt'), not the likelihood ratio test ('lrt'), gives you beta (the effect size; an estimate of the log fold change). This is because the wald test actually tests the null hypothesis of beta being 0. The LRT tests the likelihood of a complex model vs. the likelihood of a simpler model and gives you ratio of those two likelihoods -- it can't give you an estimate of the log fold change therefore sleuth doesn't report it with the LRT.