This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Volcano plot with aggregated p values from Sleuth

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.

de kallisto log2fc rna-seq sleuth
• 3,379 views • link
updated by dsull dsull • written by pv342 • pv342

1 answer

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 )

• 0 views • link
written by dsull dsull

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!

• 0 views • link
written by pv342 • pv342

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

• 0 views • link
written by dsull dsull

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

• 0 views • link
written by pv342 • pv342

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.

• 0 views • link
written by dsull dsull

Log in to answer this question.

More posts like this