GPTomics/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

| Question | Recommended test | Function | |----------|------------------|----------| | 2 unpaired groups, normal, N≥30 | Welch t-test | t.test(), stat_compare_means(method='t.test') | | 2 unpaired groups, non-normal or small N | Mann-Whitney U (Wilcoxon rank-sum) | wilcox.test(), stat_compare_means(method='wilcox.test') | | 2 paired groups | Paired t-test OR Wilcoxon signed-rank | paired = TRUE | | 3+ groups, normal | One-way ANOVA + Tukey HSD post-hoc | aov(), TukeyHSD() | | 3+ groups, non-normal | Kruskal-Wallis + Dunn post-hoc | kruskal.test(), dunn.test() | | Nested data (cells in patients) | Linear mixed model | lme4::lmer() | | Time-course / repeated measures | Repeated-measures ANOVA OR LMM | nlme::lme() | | Two-way factorial | Two-way ANOVA + interaction term | aov(y ~ a*b) | | Survival / time-to-event | Log-rank, NOT t-test | survdiff() | | Categorical outcome | Chi-square OR Fisher exact | chisq.test(), fisher.test() |

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.

Asterisks shown but p-values not reported anywhere

Trigger: label = 'p.signif' exclusively.

Mechanism: Reader cannot recover the actual p-value.

Symptom: Reviewer asks for exact p; not in figure or supplementary.

Fix: Either show numeric p (label = 'p.format') or include test results table in supplementary.

n.s. annotation hidden

Trigger: Showing only significant brackets, omitting non-significant.

Mechanism: Selective reporting biases interpretation.

Symptom: Reader assumes untested pairs were significant.

Fix: Either annotate all comparisons (with n.s. for non-significant) OR pre-specify which pairs are tested in the legend/caption.

Reading effect size from p-value

Trigger: Conclusion "highly significant difference" from p = 1e-10 on a tiny effect.

Mechanism: Large N inflates significance for trivial differences.

Symptom: Effect size negligible despite extreme p.

Fix: Always report effect size (Cohen's d, Cliff's delta, median difference with CI) alongside p. The bracket should convey direction AND magnitude, not just significance.

Reconciliation: When Tests Disagree

| Pattern | Cause | Action | |---------|-------|--------| | t-test significant; Wilcoxon n.s. | Outliers driving t-test; rank test robust | Trust Wilcoxon for non-normal data | | Unpaired n.s.; paired significant | Within-subject correlation matters | Use paired if subjects are matched | | Pairwise all-significant; ANOVA n.s. | Multiple-testing inflation in pairwise | ANOVA / Kruskal-Wallis is the overall test; pairwise post-hoc only after omnibus significant | | Pseudo-replicated p < 1e-50; LMM p = 0.1 | Pseudoreplication | LMM is correct; pseudo-replicated p is meaningless | | Bonferroni-adjusted n.s.; raw p < 0.05 | Adjustment correctly identified borderline | Trust adjusted; document the test family |

Quantitative Thresholds

| Threshold | Value | Source | |-----------|-------|--------| | α for FWER control | 0.05 family-wise | Standard | | α for FDR control | 0.05 expected FDR (BH) | Benjamini-Hochberg 1995 | | Asterisk convention | * <0.05, ** <0.01, *** <0.001 | Common practice | | Bonferroni cutoff | 0.05 / K(K-1)/2 | Standard | | Holm step-down | better than Bonferroni for all K | Holm 1979 | | FDR (BH) | less strict than FWER | Genomics standard |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | Reviewer asks "why t-test?" | Default not justified | Pre-justify test choice | | Many pairwise-significant; doesn't replicate | No multiple-testing adjustment | Holm or BH | | Effect "highly significant" but tiny | Large N inflates p | Report effect size | | Asterisks only; no p-values | label = 'p.signif' exclusively | Show p.format OR provide table | | Pseudoreplication inflated p | Cells treated as independent | LMM or pseudobulk | | n.s. comparisons hidden | Selective reporting | Annotate all pre-specified pairs | | Numeric p truncated to '<2.22e-16' | R default precision | Manual formatting or report as < 2e-16 |

References

• Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 57:289-300.
• Dunn OJ. 1964. Multiple comparisons using rank sums. Technometrics 6(3):241-252.
• Holm S. 1979. A simple sequentially rejective multiple test procedure. Scand J Stat 6(2):65-70.
• Kassambara A. 2020. Practical Statistics in R for Comparing Groups: Numerical Variables. (ggpubr / rstatix tutorial).
• Wasserstein RL, Lazar NA. 2016. The ASA's statement on p-values: context, process, and purpose. Am Stat 70(2):129-133.

Related Skills

• data-visualization/distribution-plots - Underlying box/violin/raincloud
• clinical-biostatistics/categorical-tests - Chi-square / Fisher tests for categorical outcomes
• clinical-biostatistics/effect-measures - Effect size to report alongside p
• experimental-design/multiple-testing - Methods for controlling FWER and FDR