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.

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.

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.

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.

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.