Dimensionality Reduction

Dimensionality reduction is the process of replacing a measurement matrix with many columns by a smaller set of derived columns that preserve the variance, distances, or class separation you care about. A gene expression matrix with 20,000 rows and 30 samples becomes a 30 × 5 matrix of component scores, and the geometry of the samples survives well enough to plot, cluster, or regress on.

The reason to do it goes beyond plotting. In spaces of tens of thousands of dimensions, distances concentrate, nearly every pair of points becomes roughly equidistant, and correlations appear between variables that share no biology ¹. Reducing dimensions restores a space where nearest neighbors mean something.

How it works

Two broad approaches, and you will use both.

Matrix factorization splits your data matrix X (samples × features) into a small factor matrix and a loading matrix. PCA is the unsupervised case: components are orthogonal directions of maximum variance, obtained from the SVD of the centered matrix. NMF adds nonnegativity, which makes factors additive and often more interpretable for count data. Multi-assay extensions (MOFA, JIVE, MCIA, RGCCA) find factors shared across blocks, so one factor can load on both a proteomic and a transcriptomic feature set ². A benchmark across ten cancer datasets found that MOFA, RGCCA, and MCIA recovered clinically associated and biologically annotated factors more consistently than the rest, with no single method dominating every metric ³.

Supervised methods use labels. LDA maximizes between-class scatter relative to within-class scatter, so if you already have well-defined groups and roughly Gaussian, similar-covariance classes, it will separate them in fewer dimensions than PCA. That is the whole answer to “why is LDA better than PCA”: it is not better in general, it optimizes a different objective, and it overfits badly when p ≫ n because the within-class covariance is singular and must be regularized. Sparse PLS does the supervised job while zeroing out most loadings, which gives you dimension reduction and variable selection in one pass ⁴. Filter-wrapper feature selection is a different route to the same goal, picking actual features instead of linear combinations ⁵.

Nonlinear embeddings (UMAP, t-SNE) are for visualization and neighborhood structure, not for downstream regression. Modern single-cell toolchains run them on top of 20-50 PCs rather than on raw features, which is both faster and less noisy ⁶. Contrastive-learning embeddings are now competitive on genotype data, producing 2D maps that preserve population structure better than PCA on the same variants ⁷.

In your own data

Where it shows up. Your RNA-seq deliverable includes a counts matrix (salmon.merged.gene_counts.tsv or similar). That is the input. Variance-stabilize first: DESeq2::vst() or edgeR::cpm(log=TRUE, prior.count=2). Never run PCA on raw counts, the first component will track library size. Filter to the top 1,000-2,000 most variable genes before prcomp(t(mat), center=TRUE, scale.=FALSE). Scaling is optional for log-expression and usually wrong there, since it inflates low-expressed noisy genes.

For genotypes, PCA runs on the variant matrix from your VCF. Use PLINK: plink2 --vcf sample.vcf.gz --maf 0.01 --indep-pairwise 50 5 0.2 --make-bed then --pca 10. Skipping LD pruning is the single most common mistake, because a long haplotype block will produce a component that describes one chromosomal region rather than ancestry.

What to check. Plot PC1 vs PC2 colored by draw date, sequencing run, RIN, and library prep batch before you color by anything biological. In a longitudinal single-person profile, the strongest axis is frequently the batch, not the state you were sampling. Look at the scree plot and the cumulative variance: if PC1 holds 60% and loads on hemoglobin genes, you sampled hemolysis. Inspect the top 20 loadings of each retained component by name.

Common mistakes: computing PCs on all samples and then testing group differences on those PCs with the same data, imputing missing proteomic values with zeros (which creates a component that encodes missingness), and reading UMAP inter-cluster distances as meaningful. Standardize per-block before joint factorization, otherwise the assay with the largest raw variance dominates every shared factor ⁸.

Limitations

Components are linear combinations, so a factor rarely maps to one pathway. Unsupervised methods find the largest variance, which may be technical. Nonlinear embeddings have hyperparameters (n_neighbors, min_dist, perplexity) that visibly change the picture, so report them. And with n in the dozens, any supervised reduction needs nested cross-validation or the performance estimate is fiction.

None of this is clinical interpretation. A component that separates your samples tells you structure exists, not what it means for you. Take specific findings to a clinician.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.