bio-applied-statistics-for-bioinformatics

name	bio-applied-statistics-for-bioinformatics
description	Run t-test/Mann-Whitney/ANOVA with scipy.stats, apply Bonferroni/BH-FDR via statsmodels, compute Cohen's d and power. Use for comparing expression/counts between groups, correcting p-values, or genomics power analysis.
tool_type	python
primary_tool	scipy.stats

Statistics for Bioinformatics

When to Use

• Comparing a measurement (expression, methylation, mutation count) between two or more groups and need to pick t-test vs Mann-Whitney vs ANOVA vs Kruskal-Wallis.
• Testing thousands of genes/features and need to control the false discovery rate instead of reporting raw p-values.
• Testing association between a categorical variable (mutation status, subtype) and an outcome (Fisher's exact / chi-squared).
• Reporting effect size (Cohen's d) alongside a p-value, or estimating required sample size / power before running an experiment.
• Building a volcano plot or diagnosing a p-value histogram to check whether multiple testing is behaving as expected.

Version Compatibility

Python ≥3.10, scipy ≥1.11, statsmodels ≥0.14, numpy ≥1.24, pandas ≥2.0. R ≥4.3 with base stats (no extra package needed for BH correction).

Prerequisites

pip install scipy statsmodels numpy pandas (Python) or base R (no install needed — p.adjust is in base stats). Assumes you already have per-sample values (expression, counts, methylation fraction) in an array or a tidy DataFrame — see bio-expression-matrix-counts-ingest for building that matrix.

Distribution Selection Guide

Distribution	Biology use case	Example
Normal	Continuous measurements	Log-transformed expression
Poisson	Count data (rare events)	Mutations per gene
Negative binomial	Overdispersed counts	RNA-seq reads (DESeq2 model)
Binomial	Success/failure in n trials	Methylated CpGs out of N sites
Exponential	Waiting times	Distance between mutations, survival times
Uniform	p-values under the null	Expected p-value distribution when H0 is true

Test Selection Guide

Situation	Parametric	Non-parametric
Two independent groups	Independent / Welch's t-test	Mann-Whitney U
Two paired groups	Paired t-test	Wilcoxon signed-rank
3+ independent groups	One-way ANOVA	Kruskal-Wallis
Categorical association	—	Chi-squared / Fisher's exact

Use parametric when data is ~normal (or n>30 per group) with roughly equal variances. Use non-parametric when samples are small, skewed, outlier-heavy, or ordinal.

Goal: decide between a t-test and Mann-Whitney U for a two-group comparison, and get an effect size alongside the p-value. Approach: run both tests (Welch's t-test does not assume equal variance — prefer it over Student's t-test by default), compare; if they disagree strongly, trust Mann-Whitney U for skewed/outlier data. Report Cohen's d for magnitude, since p-values alone don't convey effect size.

import numpy as np
from scipy import stats

def compare_two_groups(treated: np.ndarray, control: np.ndarray) -> dict:
    """Compare two independent groups with Welch's t-test, Mann-Whitney U, and Cohen's d.

    treated, control: 1D arrays of a continuous measurement (e.g. log2 expression).
    Returns a dict of test statistics, p-values, and effect size.
    """
    tw_stat, tw_p = stats.ttest_ind(treated, control, equal_var=False)  # Welch's t-test
    u_stat, u_p = stats.mannwhitneyu(treated, control, alternative="two-sided")

    n1, n2 = len(treated), len(control)
    s_pooled = np.sqrt(
        ((n1 - 1) * treated.std(ddof=1) ** 2 + (n2 - 1) * control.std(ddof=1) ** 2)
        / (n1 + n2 - 2)
    )
    cohens_d = (treated.mean() - control.mean()) / s_pooled

    return {
        "welch_t": tw_stat, "welch_p": tw_p,
        "mannwhitney_u": u_stat, "mannwhitney_p": u_p,
        "cohens_d": cohens_d,
    }

if __name__ == "__main__":
    rng = np.random.default_rng(42)
    control = rng.normal(7.5, 1.8, 40)   # normal tissue, log2 expression
    tumor = rng.normal(9.0, 2.0, 35)     # tumor tissue
    result = compare_two_groups(tumor, control)
    assert 0 <= result["mannwhitney_p"] <= 
    (result)

Multiple Testing Correction

Testing 20,000 genes at alpha=0.05 yields ~1,000 false positives by chance alone.

• Bonferroni: alpha_adj = alpha / m, controls family-wise error rate (FWER). Very conservative; misses true positives when m is large.
• Benjamini-Hochberg (BH/FDR): controls the expected proportion of false discoveries among rejections. Sort p-values ascending, find the largest k where p_(k) <= k/m * q. Standard choice for genomics (q=0.05).

Goal: correct thousands of raw p-values for multiple testing and count true/false positives under each method. Approach: use statsmodels.stats.multitest.multipletests for both Bonferroni and BH; a spike near p=0 in the combined p-value histogram (vs. flat/uniform elsewhere) is the diagnostic that real signal is present.

import numpy as np
from statsmodels.stats.multitest import multipletests

def correct_pvalues(p_values: np.ndarray, alpha: float = 0.05) -> dict:
    """Apply Bonferroni and Benjamini-Hochberg FDR correction to an array of p-values.

    Returns boolean reject arrays and adjusted p-values for each method.
    """
    reject_bonf, pvals_bonf, _, _ = multipletests(p_values, alpha=alpha, method="bonferroni")
    reject_bh, pvals_bh, _, _ = multipletests(p_values, alpha=alpha, method="fdr_bh")
    return {
        "bonferroni": {"reject": reject_bonf, "padj": pvals_bonf},
        "fdr_bh": {"reject": reject_bh, "padj": pvals_bh},
    }

if __name__ == "__main__":
    rng = np.random.default_rng(42)
    n_genes, n_true_de = 20000, 500
    p_null = rng.uniform(0, 1, n_genes - n_true_de)
    p_de = rng.beta(0.3, 5, n_true_de)             # small p-values for true DE genes
    p_values = np.concatenate([p_null, p_de])

    out = correct_pvalues(p_values)
    assert out["bonferroni"]["reject"].sum() <= out["fdr_bh"]["reject"].sum()  # BH always more permissive
    print(f"Bonferroni significant: {out['bonferroni']['reject'].()}")
    ()

Equivalent in R (base stats::p.adjust, no extra package):

set.seed(42)
n_genes <- 20000
n_true_de <- 500
p_null <- runif(n_genes - n_true_de)
p_de <- rbeta(n_true_de, 0.3, 5)
p_values <- c(p_null, p_de)

padj_bonf <- p.adjust(p_values, method = "bonferroni")
padj_bh <- p.adjust(p_values, method = "BH")

cat("Bonferroni significant:", sum(padj_bonf < 0.05), "\n")
cat("BH/FDR significant:    ", sum(padj_bh < 0.05), "\n")

Effect Size and Power/Sample Size

Goal: decide how many samples per group are needed to reliably detect a given effect size before running the experiment. Approach: statsmodels.stats.power.TTestIndPower — solve for nobs1 given effect_size (Cohen's d), alpha, and target power (0.8 is the conventional minimum).

from statsmodels.stats.power import TTestIndPower

def sample_size_for_effect(effect_size: float, alpha: float = 0.05, power: float = 0.8) -> float:
    """Return required sample size per group to detect `effect_size` (Cohen's d)
    at the given alpha and power for an independent two-sample t-test."""
    return TTestIndPower().solve_power(effect_size=effect_size, alpha=alpha, power=power)

if __name__ == "__main__":
    for d in [0.2, 0.5, 0.8, 1.0]:
        n = sample_size_for_effect(d)
        print(f"d={d}: n={n:.0f} per group")
    assert sample_size_for_effect(0.2) > sample_size_for_effect(1.0)  # small effects need more samples

Categorical Association

For a 2x2 (or larger) contingency table (e.g. mutation status vs drug response): use Fisher's exact test (scipy.stats.fisher_exact) for small samples/sparse tables, chi-squared (scipy.stats.chi2_contingency) for larger samples where expected counts are all >= 5.

import numpy as np
from scipy import stats

contingency = np.array([[35, 15],   # Mutated:   Responder, Non-responder
                         [20, 30]])  # Wild-type: Responder, Non-responder

odds_ratio, fisher_p = stats.fisher_exact(contingency)
chi2, chi2_p, dof, expected = stats.chi2_contingency(contingency)

Pitfalls

• Multiple testing: always apply FDR correction when testing thousands of features; raw p-values massively overstate significance.
• Skewed data + t-test: outliers distort t-test results; Mann-Whitney U is robust to them.
• P-value histogram diagnostic: flat/uniform = consistent with all-null; spike near 0 = true signal present; a spike near 1 usually indicates a test assumption violation (e.g. discreteness in Fisher's exact, ties).
• Bonferroni vs BH: Bonferroni controls FWER (probability of any false positive) and is appropriate for a handful of confirmatory tests; BH controls FDR (expected proportion of false positives among discoveries) and is the standard for genome-wide screens — using Bonferroni there discards most true positives.
• Chi-squared with small counts: if any expected cell count is < 5, use Fisher's exact instead — the chi-squared approximation breaks down.
• P-value without effect size: a tiny p-value from a huge sample can correspond to a trivial effect; always report Cohen's d (or log2 fold-change) alongside p-values.

See Also

• bio-experimental-design-multiple-testing
• bio-experimental-design-power-analysis
• bio-experimental-design-sample-size
• bio-differential-expression-de-results