# Python Bio Numpy

> Vectorize bioinformatics math with NumPy — RPKM/CPM/TPM normalization, per-gene z-scores, broadcasting over genes x samples matrices, position weight matrices (PWM/PSSM) for motif scoring, and O(n) sliding-window GC content via cumsum. Use when normalizing count matrices, computing per-row/per-column statistics on expression data, building or scoring a PWM, calculating GC content over a genome window, or replacing slow Python for-loops over arrays with vectorized operations.

- Skill: `pavel-kravchenko/python-bio-numpy` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/python-bio-numpy`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/python-bio-numpy/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/python-bio-numpy

---


# NumPy for Bioinformatics

## When to Use

- Normalizing raw read counts to RPKM/CPM/TPM without looping gene-by-gene.
- Computing per-gene or per-sample statistics (mean, std, z-score) on a genes x samples matrix and need to reason about `axis=0` vs `axis=1`.
- Building a Position Weight Matrix (PWM/PSSM) from aligned motif sequences and scoring candidate sequences against it.
- Computing GC content in a sliding window across a long sequence and a naive per-window loop is too slow.
- Replacing a Python `for` loop over an array with a vectorized NumPy expression for a 10-100x speedup.

## Version Compatibility

NumPy >= 1.26, Python >= 3.10. Examples also use `scipy.stats` (SciPy >= 1.11) for the optional t-test snippet — install separately if needed.

## Prerequisites

- `pip install numpy` (add `scipy` if you use the t-test helper below).
- Comfortable with Python lists/loops; no prior NumPy needed, but understanding shapes/dtypes helps.
- Related: `bio-expression-matrix-counts-ingest` for loading raw count files into arrays/DataFrames first.

## Key Concepts

- **Vectorized operations** run in compiled C — `arr * 2` is much faster than `[x * 2 for x in arr]`.
- **Broadcasting aligns from the right.** `(4,3) - (3,)` works (broadcasts across rows). `(4,3) - (4,)` fails — reshape to `(4,1)` first, e.g. via `keepdims=True` on a reduction.
- **Views vs copies:** `arr[2:5]` is a view (modifies original). `arr[[0,2,4]]` and `arr[arr > 0]` return copies. Call `.copy()` when you need to be explicit.
- **`axis` in aggregations:** `axis=0` collapses rows (per-column stats); `axis=1` collapses columns (per-row stats). In a genes x samples matrix: `axis=0` = per-sample, `axis=1` = per-gene.

## Normalizing Expression Counts

**Goal:** convert raw read counts into RPKM, CPM, or TPM so samples/genes are comparable.
**Approach:** vectorize the per-million and per-kilobase scaling instead of looping over genes; TPM differs from RPKM by normalizing per-kilobase *before* the per-million scaling.

```python
import numpy as np

def compute_rpkm(counts, lengths):
    """RPKM = reads per kilobase per million mapped reads.

    counts: 1D array of raw read counts per gene
    lengths: 1D array of gene lengths in bp
    """
    rpm = counts * (1_000_000 / np.sum(counts))  # per-million scaling
    return rpm * (1_000 / lengths)                # per-kilobase scaling


def compute_tpm(counts, lengths):
    """TPM = transcripts per million (normalize length first, then depth)."""
    rpk = counts / (lengths / 1_000)   # reads per kilobase
    return rpk / rpk.sum() * 1_000_000


# CPM via broadcasting a (genes, samples) matrix against per-sample totals
expr = np.array([[120, 135, 128], [45, 50, 48], [300, 280, 310]], dtype=float)
sample_totals = expr.sum(axis=0, keepdims=True)   # shape (1, 3)
cpm = expr / sample_totals * 1_000_000
```

## Per-Gene Statistics Across Samples

**Goal:** z-score each gene's expression across samples for heatmaps/clustering.
**Approach:** reduce along `axis=1` (across samples) with `keepdims=True` so broadcasting re-aligns against the original `(genes, samples)` shape.

```python
gene_means = expr.mean(axis=1, keepdims=True)  # shape (genes, 1)
gene_stds = expr.std(axis=1, keepdims=True)
z_scores = (expr - gene_means) / gene_stds

# log-transform raw counts (pseudocount avoids log(0))
log_counts = np.log2(expr + 1)

# Pearson correlation between two gene expression profiles
gene_a, gene_b = expr[0], expr[1]
corr = np.corrcoef(gene_a, gene_b)[0, 1]
```

## Position Weight Matrix (PWM) for Motif Scoring

**Goal:** build a PWM from aligned binding-site sequences and score candidate sequences against it.
**Approach:** count nucleotide frequencies per position into a `(4, seq_len)` array, then score a new sequence as the sum of log-odds vs. a uniform 0.25 background.

```python
def build_pwm(sequences):
    """Build a PWM from aligned sequences of equal length.

    Returns a (4, seq_len) frequency matrix (rows: A, T, G, C).
    """
    mapping = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
    seq_len = len(sequences[0])
    counts = np.zeros((4, seq_len))
    for seq in sequences:
        for i, nuc in enumerate(seq):
            counts[mapping[nuc], i] += 1
    return counts / len(sequences)


def score_sequence(pwm, sequence):
    """Score a sequence against a PWM (sum of log-odds vs uniform 0.25 background)."""
    mapping = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
    score = 0.0
    for i, nuc in enumerate(sequence):
        freq = pwm[mapping[nuc], i]
        score += np.log2(freq / 0.25) if freq > 0 else -10  # penalize zero-freq positions
    return score


sp1_sites = ['GGGCGG', 'GGGCGG', 'GGGCGA', 'GGGCGG', 'AGGCGG']
pwm = build_pwm(sp1_sites)
best = max(['GGGCGG', 'AAATTT', 'GGGCGA'], key=lambda s: score_sequence(pwm, s))
```

## Sliding-Window GC Content (O(n) via cumsum)

**Goal:** compute GC content over every fixed-size window of a long sequence without an O(n*window) nested loop.
**Approach:** convert bases to a 0/1 GC indicator array, take a cumulative sum, then subtract offset cumsums to get each window's sum in O(1) per window.

```python
def sliding_gc(sequence, window=50):
    """O(n) sliding-window GC content using the cumsum trick.

    Returns an array of length len(sequence) - window + 1, one value per window (%).
    """
    gc_binary = np.array([1 if n in 'GC' else 0 for n in sequence])
    cumsum = np.insert(np.cumsum(gc_binary), 0, 0)
    window_sums = cumsum[window:] - cumsum[:-window]
    return window_sums / window * 100
```

## Pitfalls

- Forgetting `keepdims=True` on a reduction before broadcasting it back against the original matrix — `expr.mean(axis=1)` returns shape `(genes,)`, which broadcasts wrong against `(genes, samples)`; you need `(genes, 1)`.
- Mixing up `axis=0`/`axis=1` — in a genes x samples matrix, `axis=0` is "per-sample" (collapses genes), `axis=1` is "per-gene" (collapses samples). Getting this backwards silently produces the wrong-shaped, wrong-meaning result (no error).
- Dividing by zero when a gene has zero counts across all samples/lengths — RPKM/TPM/CPM will produce `inf`/`nan`; filter zero-count genes first or add a pseudocount.
- Using `arr[mask] = value` on a view-derived slice expecting it to persist — fancy/boolean indexing on the *read* side returns a copy, but assignment through it still mutates the original array; know which operation you're doing.
- `score_sequence` here assumes uppercase A/T/G/C only — lowercase bases or ambiguity codes (N, R, Y...) will raise a `KeyError`; validate/uppercase input first.

## See Also

- `bio-expression-matrix-counts-ingest` — loading raw count matrices before normalizing them.
- `bio-expression-matrix-sparse-handling` — when the matrix is too large/sparse for dense NumPy arrays.
- `bio-sequence-manipulation-motif-search` — complementary motif-finding without a PWM.
- `pydeseq2` — for statistically rigorous differential expression instead of manual log2FC/t-tests.