Applied - Expression Matrices with NumPy

Core Principle: High-performance processing of gene expression matrices leverages 2D NumPy array vectorization, broadcasting, and axis reductions for count normalization, filtering, and scaling.

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

Setup

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.

Normalisation

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

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

float32 gives ~7 significant digits, far beyond any assay's precision. But accumulate in 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

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.

See also

Axes and Reductions · Broadcasting · Linear Algebra with NumPy · Boolean Masking · Applied - Differential Expression Volcano Plot · Applied - Heatmaps and Clustermaps · Tidy Data

Exercise

Run the lesson's normalise-transform-filter-reduce pipeline end to end with NumPy, without a Python loop over genes or samples.

Write pipeline(text), returning a list of output lines. The first input line holds an integer k. The second is the header: the literal gene followed by the sample names, tab separated. Every later line is a gene name followed by its raw counts. Compute each sample's library size first and drop any sample whose library size is zero -- CPM is undefined there, and dividing anyway yields silent NaNs that poison every later reduction. Over the remaining samples compute CPM, then log2(cpm + 1), then keep only genes with CPM above 1 in at least two samples.

Print samples=<names> and dropped=<names> as comma-separated lists, using none for an empty list, then genes_in=<n> genes_kept=<n>. Then print top=<names> holding the k kept genes with the largest variance of log2(cpm + 1) across samples, ordered by variance descending and gene name ascending as a tie-break, or none when no gene is kept. Then print mean_log2cpm=<v>, the mean of the kept genes' transformed values rounded to three decimals, using 0.0 when nothing is kept. Finally print corr=<v>, the Pearson correlation between the first two retained samples over the kept genes rounded to three decimals; print corr=NA when there are fewer than two retained samples, fewer than two kept genes, or either column has zero variance. Use the population variance, np.var with its default ddof.