Foundations Probability

Probability for Bioinformatics

When to Use

• Modeling variant allele read counts (germline vs somatic) with a Binomial distribution.
• Computing a p-value for "is this gene mutated more than background" with Poisson.
• Explaining/simulating RNA-seq count overdispersion with Negative Binomial (the model behind DESeq2/edgeR).
• Fitting or sampling from a Normal distribution for continuous measurements (log2 expression, quality scores).
• Distinguishing PDF vs PMF vs CDF vs survival function when a user asks "what's the probability of X".

Version Compatibility

• Python ≥3.10, scipy ≥1.11 (scipy.stats), NumPy ≥1.24, pandas ≥2.0, statsmodels ≥0.14.
• R ≥4.3 base stats (stats::dbinom, dpois, dnbinom — no extra package needed).

Prerequisites

• pip install numpy pandas scipy matplotlib statsmodels
• Familiarity with basic distribution shapes (mean/variance) helps but is not required.

scipy.stats Unified API

Goal: Model a continuous measurement (e.g. log2 expression of a housekeeping gene) and answer probability questions about it.

Approach: Build a scipy.stats.norm object from loc/scale (not mean/std), then use .sf, .ppf, .cdf directly — never invert the CDF by hand.

import numpy as np
import scipy.stats as stats

def gene_expression_probabilities(mu_expr=8.5, sigma_expr=1.2, threshold=10.0):
    """Model log2(TPM+1) of a housekeeping gene as Normal(mu, sigma) and
    return key tail/quantile probabilities.
    """
    gene_dist = stats.norm(loc=mu_expr, scale=sigma_expr)
    return {
        "prob_above_threshold": gene_dist.sf(threshold),   # P(expression > threshold)
        "percentile_95": gene_dist.ppf(0.95),              # 95th percentile
        "cdf_at_7": gene_dist.cdf(7),                       # P(expression <= 7)
    }

if __name__ == "__main__":
    result = gene_expression_probabilities()
    print(result)

Goal: Distinguish germline heterozygous SNPs from somatic (subclonal) mutations using alt-read counts at fixed coverage.

Approach: Model alt-supporting reads as Binomial(n=coverage, p=VAF); compare the PMF under a germline (p=0.5) vs somatic (p<0.5) hypothesis to get a likelihood ratio.

import numpy as np
import scipy.stats as stats

def variant_likelihood_ratio(alt_reads, coverage=30, germline_vaf=0.5, somatic_vaf=0.2):
    """Compare P(observed alt reads) under a germline vs somatic Binomial model.

    Returns the likelihood ratio (somatic / germline); >>1 favors somatic origin.
    """
    germ_dist = stats.binom(n=coverage, p=germline_vaf)
    somatic_dist = stats.binom(n=coverage, p=somatic_vaf)
    p_germ = germ_dist.pmf(alt_reads)
    p_somatic = somatic_dist.pmf(alt_reads)
    return p_somatic / p_germ

if __name__ == "__main__":
    lr = variant_likelihood_ratio(alt_reads=6)
    print(f"Likelihood ratio (somatic/germline) for 6/30 alt reads: {lr:.1f}x")

Goal: Flag a gene as a likely mutation "driver" by testing observed mutation count against a Poisson background rate.

Approach: lambda_background = gene_length_kb * mutation_rate_per_kb; use .sf(k-1) to get P(X >= k) as a one-sided p-value.

import scipy.stats as stats

def driver_gene_pvalue(observed_mutations, lambda_background):
    """One-sided Poisson test: P(X >= observed_mutations | background rate).

    Small p-value => gene has more mutations than expected by chance (candidate driver).
    """
    bg_dist = stats.poisson(mu=lambda_background)
    return bg_dist.sf(observed_mutations - 1)   # P(X >= k) = sf(k-1)

if __name__ == "__main__":
    pval = driver_gene_pvalue(observed_mutations=5, lambda_background=0.3)
    print(f"P(mutations >= 5 | background): {pval:.2e}")

Goal: Explain/simulate why RNA-seq read counts are overdispersed (variance > mean) and need a Negative Binomial, not Poisson (this is the model behind DESeq2/edgeR).

Approach: Parameterize NB by mean mu and dispersion alpha (var = mu + alpha*mu²); convert to SciPy's nbinom(n, p) form.

import scipy.stats as stats

def negative_binomial_from_mean_dispersion(mu, alpha):
    """Build a scipy.stats.nbinom distribution from mean (mu) and dispersion (alpha),
    matching the DESeq2/edgeR parameterization: variance = mu + alpha * mu**2.
    """
    n_nb = 1 / alpha
    p_nb = 1 / (1 + alpha * mu)
    return stats.nbinom(n=n_nb, p=p_nb)

if __name__ == "__main__":
    nb_dist = negative_binomial_from_mean_dispersion(mu=50, alpha=0.1)
    pois_dist = stats.poisson(mu=50)
    print(f"NB variance: {nb_dist.var():.1f}  Poisson variance: {pois_dist.var():.1f}")
    print(f"Overdispersion: {nb_dist.var() / nb_dist.mean():.1f}x")

R Equivalent

R's base stats package uses d/p/q/r prefixes (density, CDF, quantile, random) instead of scipy's .pmf/.cdf/.ppf/.rvs — the same distributions, different naming convention.

# Binomial: exact one-sided test for allele frequency (analogous to a VAF check)
# P(X > 54) for Binomial(100, 0.25), i.e. more alt reads than expected under H0
pbinom(54, size = 100, prob = 0.25, lower.tail = FALSE)

# Poisson: probability of >=5 mutations given background rate lambda = 0.3
ppois(4, lambda = 0.3, lower.tail = FALSE)   # P(X >= 5) = P(X > 4)

# Negative Binomial: RNA-seq count model, mean mu and dispersion size (~1/alpha)
mu <- 50; size <- 10   # size = 1/alpha in the mu-dispersion parameterization
dnbinom(x = 40:60, size = size, mu = mu)

# Exact binomial hypothesis test (e.g. allele-balance test at a heterozygous site)
binom.test(x = 18, n = 30, p = 0.5, alternative = "two.sided")

Descriptive Statistics Across Groups

import numpy as np
import pandas as pd

def simulate_tissue_expression(n_samples_per_tissue=20, n_genes=500, seed=42):
    """Simulate log-normal-ish expression matrices for 3 tissues with
    tissue-specific mean offsets, and return a per-gene summary DataFrame.
    """
    rng = np.random.default_rng(seed)
    tissue_offsets = {"liver": 0.0, "kidney": 0.3, "brain": -0.2}
    expr_data = {
        tissue: rng.normal(loc=8.0 + offset, scale=1.5, size=(n_samples_per_tissue, n_genes))
        for tissue, offset in tissue_offsets.items()
    }
    gene0 = {tissue: expr_data[tissue][:, 0] for tissue in tissue_offsets}
    return pd.DataFrame(gene0)

if __name__ == "__main__":
    summary_df = simulate_tissue_expression()
    print(summary_df.describe().round(3))

Pitfalls

• PDF vs PMF vs CDF: pdf(x) is a density, not a probability (can exceed 1). pmf(k) IS a probability. CDF P(X ≤ x) is always a probability.
• scipy.stats uses loc/scale, not mean/std: for exponential, scale = 1/rate.
• Use .sf(x) instead of 1 - .cdf(x) for tail p-values — avoids catastrophic cancellation when the tail probability is tiny.
• Negative Binomial parameterization varies: NumPy, SciPy (n,p), and R (size,mu or size,prob) all differ — always convert through mean/dispersion, never copy n/p values across libraries.
• Independence is a strong assumption: gene expression levels are correlated (co-regulation, pathway structure); treating genes as independent inflates false-positive rates in downstream testing.
• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based inclusive — mixing them causes off-by-one errors when computing empirical probabilities over intervals.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when computing p-values for thousands of genes/variants simultaneously — a single-test p-value threshold will overcall hits.

See Also

• bio-differential-expression-deseq2-basics — applies the Negative Binomial model here to real RNA-seq differential expression.
• bio-variant-calling-vcf-basics — VCF genotype/AD/DP fields that feed the Binomial VAF model.
• bio-experimental-design-power-analysis — power/sample-size calculations built on these same distributions.
• statistical-analysis — hypothesis testing and confidence intervals built on top of these distributions.