# Foundations Probability

> Model bioinformatics data with probability distributions (Normal, Binomial, Poisson, Negative Binomial) using scipy.stats — compute PMF/PDF/CDF/SF/PPF, simulate variant allele counts, mutation rates, and RNA-seq overdispersion. Use when doing probability calculations, distribution fitting, p-value derivation from a distribution, VAF/allele-count modeling, or explaining why RNA-seq counts need a negative binomial instead of Poisson.

- Skill: `pavel-kravchenko/foundations-probability` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/foundations-probability`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/foundations-probability/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: AI & ML
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Method | Returns | Use when |
|--------|---------|----------|
| `.pdf(x)` / `.pmf(x)` | Density (continuous) or probability (discrete) at x | Plotting distribution shape |
| `.cdf(x)` | P(X ≤ x) | Computing tail probabilities |
| `.sf(x)` | P(X > x) = 1 - CDF | One-tailed p-values directly (more numerically stable in the tail than `1-cdf`) |
| `.ppf(q)` | x such that P(X ≤ x) = q | Finding critical values / quantiles |
| `.rvs(size)` | Random samples | Simulation and bootstrap |
| `.fit(data)` | MLE parameter estimates | Fitting to real data |

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

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

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

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

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

```r
# 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

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