Applied - Expression Matrices with NumPy

In one line: the canonical NumPy exercise in bioinformatics — normalise, transform, filter and reduce a genes × samples matrix without a single Python loop.

Uses: ndarray · Broadcasting · Axes and Reductions · Boolean Masking · Linear Algebra with NumPy · Vectorization

Setup

import numpy as np, pandas as pd

expr = pd.read_csv("counts.tsv", sep="\t", index_col=0)
counts = expr.to_numpy(dtype=np.float32)          # (n_genes, n_samples)
genes = expr.index.to_numpy()
samples = expr.columns.to_numpy()

counts.shape        # (20000, 500) — genes are ROWS

Hold the orientation in your head. Genes as rows is the bulk RNA-seq convention; scanpy/AnnData use the transpose (cells × genes). Getting it backwards is the most consequential error in this note, because most operations still run and produce plausible numbers.

The habit that prevents it: print the shape of every result and check it against what it should mean.

counts.mean(axis=1).shape       # (20000,)  = one per GENE   ✓
counts.sum(axis=0).shape        # (500,)    = one per SAMPLE ✓

Normalisation

# 1. library size per sample — sum DOWN the genes
lib = counts.sum(axis=0)                    # (500,)

# 2. counts per million
cpm = counts / lib * 1e6                     # (g,s)/(s,) broadcasts ✓

# 3. log transform with a pseudocount
log_cpm = np.log2(cpm + 1)

# 4. centre each gene
gene_mean = log_cpm.mean(axis=1, keepdims=True)      # (20000, 1)
centred = log_cpm - gene_mean

# 5. z-score each gene
gene_sd = log_cpm.std(axis=1, ddof=1, keepdims=True)
z = centred / np.where(gene_sd > 0, gene_sd, 1)      # guard zero-variance genes

Steps 2 and 4 show the Broadcasting asymmetry: (g,s) / (s,) works because shapes right-align, but (g,s) - (g,) does not. keepdims=True removes the asymmetry and is the habit worth forming.

The np.where(gene_sd > 0, gene_sd, 1) guard matters: a gene with zero variance (never detected, or saturated) gives 0/0 = NaN, and one NaN poisons every downstream reduction. See Missing Data in NumPy.

Gene filtering

detected = (counts > 0).sum(axis=1)                   # in how many samples?
expressed = (cpm > 1).sum(axis=1) >= 3                # >1 CPM in ≥3 samples
keep = expressed & (counts.sum(axis=1) >= 10)

counts_f = counts[keep]
genes_f = genes[keep]                                  # ← THE SAME MASK

print(f"kept {keep.sum():,} of {len(keep):,} genes ({keep.mean():.1%})")
assert counts_f.shape[0] == len(genes_f)

Apply the mask to the labels in the same breath. Losing the row–label correspondence is the most invisible bug in matrix-based bioinformatics: every result is misattributed to the wrong gene, and nothing errors. The assertion is cheap insurance.

This is precisely the problem Pandas solves with an index — which is a good argument for doing the filtering in pandas and only dropping to NumPy for the maths.

Highly variable genes

mean = log_cpm.mean(axis=1)
var = log_cpm.var(axis=1, ddof=1)

# variance depends on mean; fit the trend and take the residual
from numpy.polynomial import Polynomial
fit = Polynomial.fit(mean, var, deg=2)
resid = var - fit(mean)

hvg_idx = np.argpartition(resid, -2000)[-2000:]       # top 2000, O(n)
hvg_idx = hvg_idx[np.argsort(resid[hvg_idx])[::-1]]   # sort just those
hvg = genes[hvg_idx]

argpartition for a top-k selection out of 20,000 avoids a full sort. See Fancy Indexing.

The mean–variance trend fit matters biologically: raw variance selects only highly expressed genes, because count noise scales with the mean. The residual selects genes that are more variable than expected for their expression level.

Correlation and PCA

# sample-sample correlation, for QC — NOT gene-gene (that's 20000², 3.2 GB)
sample_corr = np.corrcoef(log_cpm, rowvar=False)      # (500, 500)

# PCA via SVD
X = log_cpm[hvg_idx].T                                 # (samples, genes) ← transpose!
X = X - X.mean(axis=0, keepdims=True)                  # centre each gene

U, S, Vt = np.linalg.svd(X, full_matrices=False)
scores = U * S                                          # (n_samples, n_components)
var_explained = S**2 / np.sum(S**2)

print(f"PC1 {var_explained[0]:.1%}, PC2 {var_explained[1]:.1%}")

Two things people get wrong:

Orientation. PCA convention is observations as rows. Expression matrices are genes as rows. Transpose.
Centring. Without it, PC1 is just the mean expression profile.

Always print the variance explained and put it on the axis labels. A PC1 explaining 45% and one explaining 4% look identical on a plot and mean entirely different things. See Linear Algebra with NumPy.

Distance matrices

from scipy.spatial.distance import pdist, squareform

d = squareform(pdist(log_cpm[hvg_idx].T, metric="correlation"))    # (500, 500)

The pure-broadcasting version materialises an (n, n, d) intermediate:

diff = X[:, None, :] - X[None, :, :]        # (500, 500, 2000) = 4 GB in float32

Correct, elegant, and a memory disaster. Use scipy.spatial.distance for real distance work. This is the standard example of broadcasting being free for the operands and expensive for the result. See Broadcasting.

Differential expression, the simple version

group = np.array([s.startswith("treated") for s in samples])

mean_t = log_cpm[:, group].mean(axis=1)
mean_c = log_cpm[:, ~group].mean(axis=1)
log2fc = mean_t - mean_c                              # already in log space

from scipy import stats
t, p = stats.ttest_ind(log_cpm[:, group], log_cpm[:, ~group], axis=1)

from statsmodels.stats.multitest import multipletests
_, padj, _, _ = multipletests(p, method="fdr_bh")

de = pd.DataFrame({"gene": genes, "log2fc": log2fc, "pvalue": p, "padj": padj})

stats.ttest_ind(..., axis=1) runs 20,000 tests in one vectorised call. That is the whole point of this note.

But use DESeq2, edgeR or limma-voom for real work. A t-test on log-CPM ignores the count nature of the data, the mean–variance relationship, and the small-sample variance instability that shrinkage estimators exist to fix. With n=3 per group, a t-test's variance estimate is essentially noise. Use pydeseq2 if you want to stay in Python.

Batch correction, conceptually

# naive: centre within each batch
for b in np.unique(batch):
    m = batch == b
    log_cpm[:, m] -= log_cpm[:, m].mean(axis=1, keepdims=True)

This removes the batch mean and also removes any real biological difference that is confounded with batch. If all your treated samples were run in batch 2, this deletes your effect.

Check for confounding before correcting anything:

pd.crosstab(batch, condition)      # ← if this is diagonal, you cannot correct

A perfectly confounded design is not fixable by any method. Use pycombat / limma::removeBatchEffect with the biological covariate protected, and only when the design permits it.

Memory

counts.nbytes / 1e6                       # float64: 80 MB per 20000×500
counts.astype(np.float32).nbytes / 1e6    # 40 MB

np.log2(cpm, out=cpm)                     # transform in place — no second copy

float32 gives ~7 significant digits, far beyond any assay's precision. But accumulate in float64:

counts.sum(axis=0, dtype=np.float64)

For single-cell (100k cells × 30k genes, >90% zeros), a dense array is 12 GB and mostly zeros. Use scipy.sparse.

The full pipeline

def process(counts, genes, samples, min_cpm=1, min_samples=3, n_hvg=2000):
    lib = counts.sum(axis=0)
    cpm = counts / lib * 1e6
    keep = (cpm > min_cpm).sum(axis=1) >= min_samples
    counts, cpm, genes = counts[keep], cpm[keep], genes[keep]
    print(f"kept {keep.sum():,}/{len(keep):,} genes")

    log_cpm = np.log2(cpm + 1)

    mean, var = log_cpm.mean(axis=1), log_cpm.var(axis=1, ddof=1)
    resid = var - np.poly1d(np.polyfit(mean, var, 2))(mean)
    hvg = np.argpartition(resid, -n_hvg)[-n_hvg:]

    X = log_cpm[hvg].T
    X = X - X.mean(axis=0, keepdims=True)
    U, S, Vt = np.linalg.svd(X, full_matrices=False)

    return {"log_cpm": log_cpm, "genes": genes, "samples": samples,
            "hvg": genes[hvg], "pcs": U * S,
            "var_explained": S**2 / (S**2).sum()}

Common mistakes

axis=0 vs axis=1 reversed. Plausible, silent, wrong.
Filtering the matrix but not the gene labels.
PCA without transposing or centring.
Not reporting variance explained.
Zero-variance genes producing NaN.
Materialising an (n,n,d) distance intermediate.
A t-test where a proper DE method belongs.
Batch-correcting a confounded design.
int16 for counts — highly expressed genes overflow.
Dense arrays for single-cell data.

Applied - Expression Matrices with NumPy

Applied - Expression Matrices with NumPy

Setup

Normalisation

Gene filtering

Highly variable genes

Correlation and PCA

Distance matrices

Differential expression, the simple version

Batch correction, conceptually

Memory

The full pipeline

Common mistakes

See also