Foundations Statistics Python

Statistics with Python

When to Use

• Comparing two or more groups of gene expression, protein abundance, or other continuous measurements (t-test, Mann-Whitney U, ANOVA, Kruskal-Wallis).
• Testing categorical/count associations: Hardy-Weinberg genotype counts (chi-square), SNP-disease 2x2 tables (Fisher's exact).
• Modeling variant allele counts, somatic mutation counts, or RNA-seq read counts with binomial/Poisson/negative-binomial distributions.
• Correlating two continuous variables (mRNA vs protein) or fitting linear/multiple regression with statsmodels.
• Correcting p-values across thousands of genes (Bonferroni/BH-FDR) or computing required sample size/power before an experiment.

Version Compatibility

scipy ≥1.11, statsmodels ≥0.14, numpy ≥1.24, pandas ≥2.0, Python ≥3.10.

Prerequisites

pip install scipy statsmodels numpy pandas matplotlib. Assumes familiarity with probability distributions (see foundations-probability) and basic pandas/numpy.

Pitfalls

• ttest_ind defaults to Student's t-test (equal_var=True), which assumes equal variances. Pass equal_var=False explicitly for Welch's t-test — the safer default when variances may differ, as is common comparing tumor vs normal expression.
• mannwhitneyu requires alternative: always pass alternative='two-sided' explicitly — relying on the default risks version-dependent behavior.
• multipletests returns a 4-tuple: (reject, pvals_corrected, alphacSidak, alphacBonf) — unpack all four or index explicitly, don't assume only two values.
• pearsonr/spearmanr return (statistic, pvalue): unpack both — r = stats.pearsonr(x, y) silently gives you a tuple, not a float.
• Coordinate systems: BED is 0-based half-open, VCF/GFF are 1-based inclusive — mixing them in downstream stats causes off-by-one errors.
• Multiple testing: never read raw p-values from a genome-wide scan; always apply FDR (Benjamini-Hochberg) before drawing conclusions.

Core Distributions for Bioinformatics Data

Goal: model gene expression (normal), variant allele counts (binomial), somatic mutation counts (Poisson), and overdispersed RNA-seq counts (negative binomial) using scipy.stats's consistent .pdf/.pmf, .cdf, .sf, .ppf, .rvs interface.

Approach: build the distribution object once with stats.<dist>(params), then call methods on it — no need to re-derive formulas per distribution.

import numpy as np
from scipy import stats

def variant_read_pvalue(alt_reads, coverage, expected_vaf=0.5):
    """P-value for observing >= alt_reads under a binomial null model.

    Use to flag a read count as unlikely under the expected variant
    allele frequency (e.g. 0.5 for a heterozygous germline SNP).
    """
    null_dist = stats.binom(n=coverage, p=expected_vaf)
    return null_dist.sf(alt_reads - 1)  # P(X >= alt_reads)


def mutation_burden_pvalue(observed_muts, background_rate):
    """P-value for a gene's somatic mutation count under a Poisson background model.

    background_rate: expected mutations per gene (e.g. from gene length x
    genome-wide mutation rate). Flags candidate driver genes.
    """
    bg_dist = stats.poisson(mu=background_rate)
    return bg_dist.sf(observed_muts - 1)


def nb_from_mean_dispersion(mu, alpha):
    """Convert bioinformatics (mean, dispersion) to scipy's nbinom(n, p).

    var = mu + alpha * mu**2 (DESeq2/edgeR convention).
    """
    n = 1 / alpha
    p = 1 / (1 + alpha * mu)
    return stats.nbinom(n=n, p=p)


# Example: is 6 alt reads more consistent with germline (p=0.5) or somatic (p=0.2)?
germ_dist = stats.binom(n=30, p=0.5)
somatic_dist = stats.binom(n=30, p=0.2)
likelihood_ratio = somatic_dist.pmf(6) / germ_dist.pmf(6)
print(f"P(alt=6) is {likelihood_ratio:.1f}x more likely under the somatic model")

Choosing and Running the Right Hypothesis Test

Goal: pick the correct test for the study design (parametric vs non-parametric, paired vs independent, two groups vs many) and run it correctly.

Approach: check assumptions first (stats.shapiro for normality, stats.levene for equal variance), then dispatch to the matching scipy/statsmodels function.

import numpy as np
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd

def check_assumptions(group_a, group_b, alpha=0.05):
    """Test normality (Shapiro-Wilk) and equal variance (Levene) for two groups.

    Returns a dict with booleans used to choose the downstream test.
    """
    _, p_norm_a = stats.shapiro(group_a)
    _, p_norm_b = stats.shapiro(group_b)
    _, p_var = stats.levene(group_a, group_b)
    return {
        "normal": (p_norm_a > alpha) and (p_norm_b > alpha),
        "equal_var": p_var > alpha,
    }


def compare_two_groups(group_a, group_b, alpha=0.05):
    """Auto-select t-test (Student/Welch) or Mann-Whitney U based on assumptions.

    Returns (test_name, statistic, p_value).
    """
    assumptions = check_assumptions(group_a, group_b, alpha)
    if assumptions["normal"]:
        stat, p = stats.ttest_ind(group_a, group_b, equal_var=assumptions["equal_var"])
        name = "Student's t-test" if assumptions["equal_var"] else "Welch's t-test"
    else:
        stat, p = stats.mannwhitneyu(group_a, group_b, alternative="two-sided")
        name = "Mann-Whitney U"
    return name, stat, p


# Multi-group comparison: ANOVA -> Tukey HSD post-hoc if significant
liver, kidney, brain = (np.random.normal(m, 1.2, 20) for m in (8.0, 8.3, 7.8))
f_stat, p_anova = stats.f_oneway(liver, kidney, brain)
if p_anova < 0.05:
    labels = ["liver"] * 20 + ["kidney"] * 20 + ["brain"] * 20
    tukey = pairwise_tukeyhsd(np.concatenate([liver, kidney, brain]), labels, alpha=0.05)
    print(tukey.summary())

# Categorical association: chi-square (large n) or Fisher's exact (small n / 2x2)
observed = np.array([[38, 12], [22, 28]])   # [disease, control] x [alt allele, no alt allele]
odds_ratio, p_fisher = stats.fisher_exact(observed, alternative="two-sided")
print(f"Odds ratio: {odds_ratio:.3f}, p={p_fisher:.4f}")

Correlation, Regression, and Multiple Testing Correction

Goal: quantify a relationship between two continuous variables and correct p-values when testing many features (e.g. every gene) at once.

Approach: use pearsonr/spearmanr for a single pair, statsmodels.formula.api.ols for regression with multiple predictors, and multipletests for FDR correction across a p-value array.

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests

def correlate_and_regress(mrna, protein):
    """Pearson + Spearman correlation, then simple OLS: protein ~ mrna.

    Spearman is more robust when protein data contain outliers.
    """
    r_pearson, p_pearson = stats.pearsonr(mrna, protein)
    r_spearman, p_spearman = stats.spearmanr(mrna, protein)
    df = pd.DataFrame({"mrna": mrna, "protein": protein})
    model = smf.ols("protein ~ mrna", data=df).fit()
    return {
        "pearson_r": r_pearson, "pearson_p": p_pearson,
        "spearman_r": r_spearman, "spearman_p": p_spearman,
        "ols_r_squared": model.rsquared, "ols_slope_p": model.pvalues["mrna"],
    }


def fdr_correct(pvalues, alpha=0.05, method="fdr_bh"):
    """Benjamini-Hochberg FDR correction across many p-values (e.g. per-gene DE tests).

    Returns a DataFrame with the original and adjusted p-values plus a
    significance flag — safe to sort/filter directly.
    """
    reject, pvals_corrected, _, _ = multipletests(pvalues, alpha=alpha, method=method)
    return pd.DataFrame({
        "pvalue": pvalues,
        "padj": pvals_corrected,
        "significant": reject,
    })


pvals = np.random.beta(0.5, 5, size=1000)  # simulated per-gene p-values
results = fdr_correct(pvals)
print(f"{results['significant'].sum()} genes significant at FDR 0.05")

See Also

• foundations-probability — probability distributions and the theory behind the tests here.
• foundations-biostatistics-fundamentals — conceptual grounding for hypothesis testing.
• bio-experimental-design-multiple-testing — deeper coverage of FDR/Bonferroni strategies.
• bio-experimental-design-power-analysis — sample-size and power calculations with statsmodels.stats.power.