# Bio Applied Statistics For Bioinformatics

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

- Skill: `pavel-kravchenko/bio-applied-statistics-for-bioinformatics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-statistics-for-bioinformatics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-statistics-for-bioinformatics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-statistics-for-bioinformatics

---


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

```python
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"] <= 1
    print(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.

```python
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'].sum()}")
    print(f"BH/FDR significant:     {out['fdr_bh']['reject'].sum()}")
```

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

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

```python
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.

```python
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`