Bio Data Visualization Statistical Annotation

Version Compatibility

Reference examples tested with: ggpubr 0.6+, ggsignif 0.6+, rstatix 0.7+, statannotations 0.6+ (Python), seaborn 0.13+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Statistical Annotation

"Add p-values to my plot" -> Render pairwise group comparisons as brackets with the correct statistical test (parametric vs non-parametric, paired vs unpaired, independent vs nested), adjusted for multiple testing, with rendering of significance as either numerical p OR asterisks. The choices that matter: which test is appropriate for the data, what multiple-testing adjustment applies, and whether to show n.s. (non-significant) results.

• R: ggpubr::stat_compare_means, ggsignif::geom_signif, rstatix::t_test/wilcox_test
• Python: statannotations.Annotator, scipy.stats directly

The Single Most Important Modern Insight -- The Test Must Match the Data

Tool defaults are NOT data-appropriate. ggpubr::stat_compare_means(method='t.test') uses Welch's two-sample t-test assuming approximate normality and unequal variances. This is wrong when:

1. Data are non-normal and N is small (<30): use Mann-Whitney U (method='wilcox.test').
2. Data are paired: use paired t-test or paired Wilcoxon (paired = TRUE).
3. Comparing >2 groups: ANOVA / Kruskal-Wallis with post-hoc, not all-pairs t-test (multiple-testing penalty).
4. Data are nested (cells within patients, replicates within samples): linear mixed model, NOT pairwise test.

The bracket-and-asterisk visual is the same; the underlying statistics are not. Choose the test deliberately.

Decision Tree for Test Selection

Multiple Testing Adjustment

For pairwise comparisons among K groups, there are K(K-1)/2 unadjusted p-values. Without adjustment, family-wise error rate inflates rapidly:

• 3 groups: 3 comparisons; α_FW = 14% at nominal 5%
• 4 groups: 6 comparisons; α_FW = 26%
• 6 groups: 15 comparisons; α_FW = 54%

# rstatix supports per-comparison adjustment
library(rstatix)
df %>%
    pairwise_wilcox_test(value ~ group, p.adjust.method = 'bonferroni') %>%
    add_xy_position(x = 'group')

p.adjust.method options:

• 'bonferroni' — strictest; controls FWER
• 'holm' — stepwise Bonferroni; uniformly more powerful than Bonferroni
• 'BH' (Benjamini-Hochberg) — FDR; less strict than FWER; standard for genomics
• 'fdr' — alias for BH

For figure annotations, holm is the modern default — controls FWER and is more powerful than bonferroni. For a small number of pre-planned comparisons (≤3), Bonferroni is fine.

ggpubr -- Standard ggplot2 Workflow

Goal: Add per-comparison p-value brackets between groups on a distribution plot, using a test appropriate to data shape and adjusting for multiple comparisons.

Approach: Build the base plot with ggboxplot(); add stat_compare_means() with explicit method, comparisons, p.adjust.method, and label arguments; render as asterisks (p.signif) for terse display or numeric (p.format) for precise display.

library(ggpubr)

# Default boxplot + p-value bracket(s)
ggboxplot(df, x = 'group', y = 'value', color = 'group',
          add = 'jitter', palette = 'npg') +
    stat_compare_means(method = 'wilcox.test',           # explicit; default is t-test
                       comparisons = list(c('Control', 'Treatment'),
                                          c('Control', 'Vehicle'),
                                          c('Treatment', 'Vehicle')),
                       label = 'p.signif',               # 'p.signif' for asterisks; 'p.format' for numeric
                       p.adjust.method = 'holm',
                       method.args = list(alternative = 'two.sided'))

For an overall test plus pairwise:

# Overall + per-comparison
ggboxplot(df, x = 'group', y = 'value', color = 'group') +
    stat_compare_means(method = 'kruskal.test',          # overall test
                       label.y = 1.05 * max(df$value)) +
    stat_compare_means(comparisons = pairs,
                       method = 'wilcox.test',
                       p.adjust.method = 'holm',
                       label = 'p.signif')

ggsignif -- Lighter Alternative

library(ggsignif)
ggplot(df, aes(group, value, fill = group)) +
    geom_boxplot() +
    geom_signif(comparisons = list(c('Control', 'Treatment')),
                test = 'wilcox.test',
                map_signif_level = TRUE,                  # asterisks vs numeric p
                step_increase = 0.1) +
    scale_fill_manual(values = c('#0072B2', '#D55E00'))

map_signif_level = TRUE converts p-values to asterisks per Wasserstein-Lazar 2016 convention:

• *** p < 0.001
• ** p < 0.01
• * p < 0.05
• ns p ≥ 0.05

For literal p-values, set FALSE.

statannotations (Python)

import seaborn as sns
from statannotations.Annotator import Annotator

ax = sns.boxplot(x='group', y='value', data=df, palette=['#0072B2', '#D55E00', '#009E73'])

pairs = [('Control', 'Treatment'),
         ('Control', 'Vehicle'),
         ('Treatment', 'Vehicle')]

annotator = Annotator(ax, pairs, data=df, x='group', y='value')
annotator.configure(test='Mann-Whitney',                    # 't-test_ind', 't-test_paired', 'Wilcoxon', etc
                    comparisons_correction='holm',
                    text_format='star',                     # 'star', 'simple', 'full'
                    line_height=0.02,
                    text_offset=0.5)
annotator.apply_and_annotate()

Per-Method Failure Modes

Default t-test on non-normal data

Trigger: stat_compare_means(method='t.test') (default) on log-distributed expression.

Mechanism: t-test assumes approximate normality; non-normal data with small N inflates Type-I.

Symptom: Significant p where rank test gives p > 0.05.

Fix: Switch to method='wilcox.test' for non-normal or small-N data. Verify normality with Shapiro-Wilk if borderline.

Pairwise tests without adjustment

Trigger: Multiple bracket annotations with raw p-values.

Mechanism: K(K-1)/2 comparisons inflate FWER without adjustment.

Symptom: All-pairwise significant at nominal 0.05; doesn't replicate.

Fix: p.adjust.method = 'holm' (or 'bonferroni' or 'BH'). Document choice.

Paired data tested as independent

Trigger: Before/after measurements in same subjects, tested with unpaired t-test.

Mechanism: Ignores within-subject correlation; loses power.

Symptom: Non-significant p where paired test gives significant.

Fix: paired = TRUE (R) or t-test_paired (statannotations). Verify subjects are correctly matched.

Nested data tested with pairwise t

Trigger: Hundreds of cells per patient, tested as if each cell is independent.

Mechanism: Pseudoreplication — within-patient correlation ignored; p-values dramatically over-significant.

Symptom: p < 1e-50 from a dataset where the actual N is ~10 patients.

Fix: Linear mixed model (lme4::lmer(value ~ group + (1|patient_id))); aggregate to per-patient median first; or pseudobulk.