Background

Bayesian models are essential for studying complex biological systems of interest (SOIs), such as irradiation effects on metabolite concentrations or vaccine impacts on immune gene repertoires. In these models layered parameters describe key biological features of the SOIs. After model fitting, each parameter is characterized by a posterior distribution: a probability distribution representing all plausible effect values given the observed data. With thousands of such posteriors in high-dimensional analyses, identifying large, reliable effects becomes challenging.

Traditional volcano plots address this for frequentist analyses by plotting fold-changes against –log(p-values). We introduce Bayesian volcano plots that instead visualize the posterior mean effect size of parameter \(i\) (\(\theta_i\)) against the probability, \(\pi_i\), where the \(\pi_i\) quantifies the posterior probability that \(\theta_i\) is not equal the null effect. This directly highlights both magnitude and biological relevance. While Sousa et al. (2020) conceptualized Bayesian volcano plots, our R package provides the first practical implementation for automated \(\pi_i\) calculation and visualization of complex biological effects.

Installation

if (!require("BayesVolcano", quietly = TRUE)) {
  install.packages("BayesVolcano")
}

Input data

To create a Bayesian volcano plot, BayesVolcano needs two main inputs:

Posterior samples: A data frame where each column represents a biological effect (e.g., gene expression change) and each row shows posterior values.
Annotation metadata: A data frame mapping parameter names to biological entities.

If you used rstan or brms to fit your model, our extract_fit() function automatically converts your results into the right format. Just install the corresponding package (rstan or brms) first, then run extract_fit(your_model,parameter_name) to prepare your data.

Posterior Samples

data("posterior")
head(posterior[, 1:4])
#>   delta_mu.1.2.3.4 delta_mu.2.2.3.4 delta_mu.3.2.3.4 delta_mu.4.2.3.4
#> 1        1.1527980       0.24442670       -0.5283300       -2.8424959
#> 2        1.5428293       0.41885565       -0.8167788        0.2726311
#> 3        1.4201598      -0.92882009       -0.6627982       -1.8188774
#> 4        0.3277584      -0.40256872        0.3988543       -0.7784557
#> 5        0.8694004      -0.06791806       -0.9230497        0.3559686
#> 6        0.2562036      -2.11318184        1.0385491       -1.1102383

Annotation Metadata

This data frame maps parameter names to biological entities with at least a column named “parameter” and a column named “label”. Additional columns can be provided for later visualization (categorical like “group” or numerical like “value”).

data("annotation_df")
head(annotation_df)
#>          parameter                              label group     value
#> 1 delta_mu.1.2.3.4 1-Aminocyclopropanecarboxylic acid     A  9.533826
#> 2 delta_mu.2.2.3.4                    2-Aminomuconate     B 12.031100
#> 3 delta_mu.3.2.3.4             2-Phosphoglyceric acid     A  3.567905
#> 4 delta_mu.4.2.3.4              3-Hydroxybutyric acid     B  8.407870
#> 5 delta_mu.5.2.3.4                       Acetoacetate     A  9.891394
#> 6 delta_mu.6.2.3.4                  Acetyl coenzyme A     B  9.806375

Workflow

Step 1: Prepare input data

The prepare_volcano_input() function takes as input the posterior and the annotation data frame. Run this now:

d <- prepare_volcano_input(posterior = posterior, annotation = annotation_df)

The output contains

parameter: name of parameter
pi.value: calculated \(\pi\) (see section below)
null.effect: null effect value
parameter.low/high: lower and upper bounds of posterior credible interval
CrI.width: width of posterior credible interval
CrI.level: input credible interval
label: metadata-based name of parameter
group: metadata-based group of parameter
value: metadata-based value of parameter

str(d)
#> 'data.frame':    98 obs. of  11 variables:
#>  $ parameter       : chr  "delta_mu.1.2.3.4" "delta_mu.2.2.3.4" "delta_mu.3.2.3.4" "delta_mu.4.2.3.4" ...
#>  $ pi.value        : num  0.58 0.138 0.292 0.184 0.594 0.916 0.88 0.27 0.23 0.518 ...
#>  $ null.effect     : num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ parameter.median: num  0.643 -0.156 -0.271 -0.213 -0.9 ...
#>  $ parameter.low   : num  -1.23 -2.06 -1.98 -2.11 -3.31 ...
#>  $ parameter.high  : num  2.59 1.96 1.34 1.77 1.5 ...
#>  $ CrI.width       : num  3.82 4.02 3.32 3.88 4.8 ...
#>  $ CrI.level       : num  0.95 0.95 0.95 0.95 0.95 0.95 0.95 0.95 0.95 0.95 ...
#>  $ label           : chr  "1-Aminocyclopropanecarboxylic acid" "2-Aminomuconate" "2-Phosphoglyceric acid" "3-Hydroxybutyric acid" ...
#>  $ group           : chr  "A" "B" "A" "B" ...
#>  $ value           : num  9.53 12.03 3.57 8.41 9.89 ...

Customizing Input Parameters

You can specify the null effect and credible interval level:

d <- prepare_volcano_input(
  posterior = posterior,
  annotation = annotation_df,
  null.effect = 0.5,
  CrI_level = 0.7
)

Step 2: Understanding \(\pi\)

The function prepare_volcano_input() uses the posterior of each parameter from the annotations (e.g., effect size \(\theta_i\) of parameter \(i\)) to calculate \(\pi_i\):

\(\pi_i = 2 \cdot \max\left(\int_{\theta_i = -\infty}^{\bar{\theta}} p(\theta_i)\mathrm{d}\theta_i, \int_{\theta_i = \bar{\theta}}^{\infty} p(\theta_i)\mathrm{d}\theta_i\right) - 1\)

Where \(\bar{\theta}\) is the null effect. This measures the probability that the effect is in the “direction” away from the null.

\(\pi \approx 1\): Strong evidence for an effect
\(\pi \approx 0\): Negligible evidence for an effect

Important Note: A \(\pi\)-value near 0 doesn’t prove the absence of an effect - it indicates the posterior is widely distributed around the null.

Step 3: Creating the volcano plot

The plot_volcano() function creates a ggplot visual based on the formatted input:

plot_volcano(d)

This plot shows:

X-axis: Posterior median effect size \(\theta\) of parameter \(i\)
Y-axis: \(\pi\)-value
Vertical line: Null effect reference

Where each point is a parameter with median effect size (x-axis) and \(\pi\) (y-axis). The null effect is shown by the vertical line.

Adding Credible Intervals

To visualize effect size uncertainty, add credible intervals as error bars:

plot_volcano(d,
  CrI = TRUE
)

This adds:

Grey error bars showing 95% credible intervals
Subtitle describing the interval level

Visualizing Uncertainty with Point Size

You can also represent credible interval width through point size:

plot_volcano(d,
  CrI = TRUE,
  CrI_width = TRUE
)

Wide intervals = Small points
Narrow intervals = Large points

Color-Coding Parameters

Use metadata columns for color-coding:

plot_volcano(d,
  color = "group"
)

plot_volcano(d,
  color = "value"
)

Combine with other visualizations:

plot_volcano(d,
  CrI = TRUE,
  CrI_width = TRUE,
  color = "group"
)

Step 4: Customizing Your Plot

Adding Titles and Labels

# Customization requires the ggplot2 package
if (!require("BayesVolcano", quietly = TRUE)) {
  install.packages("BayesVolcano")
}
library(ggplot2)

p <- plot_volcano(d)

p + xlab("my informative parameter") +
  ggtitle("My amazing plot")

Adding Labels

Use geom_text() for basic labeling or ggrepel to avoid overplotting:

p <- plot_volcano(d)
p + geom_text(aes(label = label))


# ggrepel version
# library(ggrepel)
# p +
# geom_text_repel(aes(label=label))

Conditional Labeling

Show labels only for reliable effects:

p <- plot_volcano(d)

p + geom_text(aes(label = ifelse(abs(parameter.median) > 1.6 & # only show for parameter value > 0.5
  pi > 0.95, # only show for pi > 0.95
label, # which variable contains the label
""
))) # do not display label if outside of set ranges

This highlights:

Parameters with large effect sizes (|median| > 1.6)
Parameters with strong evidence (\(\pi\) > 0.95)

References

Julie de Sousa, Ondřej Vencálek, Karel Hron, Jan Václavík, David Friedecký, Tomáš Adam, Bayesian multiple hypotheses testing in compositional analysis of untargeted metabolomic data, Analytica Chimica Acta, Volume 1097, 2020, Pages 49-61, ISSN 0003-2670, https://doi.org/10.1016/j.aca.2019.11.006

Corresponding GitHub Repository: https://github.com/sousaju/BayesVolcano

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.2 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Berlin
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.2      BayesVolcano_1.0.1
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.2        cli_3.6.5          knitr_1.51         rlang_1.1.7       
#>  [5] xfun_0.56          otel_0.2.0         purrr_1.2.1        generics_0.1.4    
#>  [9] S7_0.2.1           jsonlite_2.0.0     labeling_0.4.3     glue_1.8.0        
#> [13] htmltools_0.5.9    sass_0.4.10        scales_1.4.0       rmarkdown_2.30    
#> [17] grid_4.5.1         tibble_3.3.1       evaluate_1.0.5     jquerylib_0.1.4   
#> [21] fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5    compiler_4.5.1    
#> [25] dplyr_1.2.0        HDInterval_0.2.4   RColorBrewer_1.1-3 pkgconfig_2.0.3   
#> [29] tidyr_1.3.2        rstudioapi_0.18.0  farver_2.1.2       digest_0.6.39     
#> [33] viridisLite_0.4.3  R6_2.6.1           tidyselect_1.2.1   pillar_1.11.1     
#> [37] magrittr_2.0.4     bslib_0.10.0       withr_3.0.2        tools_4.5.1       
#> [41] gtable_0.3.6       cachem_1.1.0

BayesVolcano

Katja Danielzik (katja.danielzik@uni-due.de)