PCA, t-SNE and UMAP on methylation data, done carefully

Almost every methylation project I have worked on has, at some point, a two-dimensional scatter plot of samples colored by tumor type. These plots are useful. They show whether a cohort has structure, whether a new sample lands near a known class, and whether something technical is going on. They are also easy to over-read, and the defaults in most tools are not chosen with methylation arrays in mind.

This post walks through how I do it: preparing the matrix, running PCA first, then t-SNE and UMAP on top of it, and the checks I run before I believe a picture. The motivating example is the Capper et al. (2018) CNS tumor reference cohort (GEO GSE90496), whose overview figure is itself a t-SNE of the reference samples. I won’t quote numbers from that dataset here; the code is written so you can run it and see your own.

Starting point: a clean beta matrix

I assume you already have a probes by samples matrix of beta values after normalization and probe QC (detection p-value filtering, removal of sex chromosome, SNP-affected and cross-reactive probes). If you are starting from IDAT files, minfi (read.metharray.exp, preprocessNoob, getBeta) or sesame (openSesame) will get you there. The sex chromosome filter matters more for embeddings than for classifiers: leave those probes in, and one of your first principal components may simply be patient sex.

You also want a sample sheet with at least the class label and whatever technical variables you have: array type (450k or EPIC), slide barcode (the first part of the IDAT file name), processing date, center, and sample material (frozen or formalin-fixed paraffin-embedded). You will need these later, and they are much harder to reconstruct after the fact.

import numpy as np
import pandas as pd

# Placeholder paths. betas: probes x samples. meta: one row per sample.
betas = pd.read_csv("betas.tsv.gz", sep="\t", index_col=0)
meta = pd.read_csv("sample_sheet.tsv", sep="\t", index_col=0)
meta = meta.loc[betas.columns]   # same order; errors if IDs don't match

Beta values or M-values

Beta values are bounded between 0 and 1 and read naturally as “fraction methylated”. M-values are the logit (base 2) of beta:

M = log2( beta / (1 - beta) )

Du et al. (2010) showed that beta values are heteroscedastic: their variance is compressed near 0 and 1, where many CpGs sit, while M-values have more uniform variance across the range. That is why M-values are usually recommended for statistical testing. For visualization the argument is less clear-cut. M-values stretch differences near the extremes, so a probe moving from 0.01 to 0.05 counts about as much as one moving from 0.3 to 0.7. Depending on what you care about, that is either a feature or noise amplification.

My default is to select probes and run PCA on betas, because the differences are easy to interpret, and to check that the main structure survives on M-values. If the picture changes a lot between the two, that is worth knowing. The conversion needs a small offset to avoid infinite values:

def beta_to_m(b, offset=1e-3):
    b = np.clip(b, offset, 1 - offset)
    return np.log2(b / (1 - b))

Selecting the most variable probes

Most probes on the array are nearly constant across samples (constitutively methylated or unmethylated) and contribute only noise to distances. Keeping the top few thousand to a few tens of thousands by standard deviation is standard practice. The exact number matters less than you might fear, which you can check by trying two or three values.

betas = betas.dropna(axis=0)              # or impute; PCA can't take NaNs
sd = betas.std(axis=1)
top = sd.sort_values(ascending=False).index[:10_000]

B = betas.loc[top].T.to_numpy(dtype=np.float64)   # samples x probes
M = beta_to_m(B)

Two things to keep in mind. First, the ranking depends on the scale: the top probes by beta standard deviation are not the same as the top probes by M-value standard deviation. Second, the selection is computed on the samples in the plot, so adding a large new batch can change which probes are selected. When you want to project new samples onto an existing embedding, keep the original probe list fixed.

Unlike in a classifier, selecting probes on the full dataset is fine here: there is no held-out performance estimate to leak into.

PCA first

I always run PCA before anything nonlinear, for three reasons. It is deterministic (up to sign) and fast, so it is the reference point. The variance explained tells you how much low-dimensional structure there is. And the top principal components make a better input to t-SNE and UMAP than thousands of raw probes: they denoise the data and make neighbor searches much faster.

Centering and scaling

scikit-learn’s PCA centers each feature but does not scale it. For methylation data I usually leave it that way. Standardizing each probe to unit variance gives every probe equal weight, which in practice means upweighting low-variance probes that you just went to the trouble of filtering against. If your features are on wildly different scales, scaling makes sense; betas are all on the same 0 to 1 scale.

from sklearn.decomposition import PCA

pca = PCA(n_components=50, svd_solver="randomized", random_state=0)
pcs = pca.fit_transform(B)

evr = pca.explained_variance_ratio_
print(np.round(evr[:10], 3))
print("cumulative, first 10 PCs:", evr[:10].sum().round(3))

The randomized solver is much faster on wide matrices and is accurate for the top components; random_state makes it reproducible.

How many components

A scree plot (variance explained per component) usually shows a few large components and then a long tail. For a cohort with dozens of classes, the useful signal is spread over more components than an “elbow” suggests, because each pair of related classes may need its own axis. A simple sanity check is a permutation null: shuffle each probe independently, which destroys correlation between probes while keeping each probe’s distribution, and see how much variance the top components explain by chance.

def parallel_analysis(X, n_components=50, n_perm=5, seed=0):
    rng = np.random.default_rng(seed)
    null = np.zeros((n_perm, n_components))
    for i in range(n_perm):
        Xp = np.column_stack([rng.permutation(col) for col in X.T])
        null[i] = PCA(n_components=n_components, svd_solver="randomized",
                      random_state=i).fit(Xp).explained_variance_ratio_
    return null.max(axis=0)

null_evr = parallel_analysis(B)
above = evr > null_evr
n_keep = int(np.argmin(above)) if not above.all() else len(above)
print("PCs above the permutation null:", n_keep)

For t-SNE and UMAP input, 30 to 50 components is a common and reasonable default for a cohort of this size. The exact number usually changes the picture very little, and you can check that too.

Looking for technical structure in the PCs

Before making any pretty plots, I check which variables the top components track. For a categorical variable, the fraction of a component’s variance explained by group means (an R² from a one-way ANOVA) is a quick summary:

def pc_covariate_r2(pcs, covariate, n=10):
    g = np.asarray(covariate).astype(str)
    out = {}
    for j in range(n):
        v = pd.Series(pcs[:, j])
        fitted = v.groupby(g).transform("mean")
        out[f"PC{j + 1}"] = 1 - ((v - fitted) ** 2).sum() / ((v - v.mean()) ** 2).sum()
    return pd.Series(out)

r2 = pd.DataFrame({
    col: pc_covariate_r2(pcs, meta[col])
    for col in ["methylation_class", "array_type", "slide", "material"]   # your column names
})
print(r2.round(2))

A caution on reading this: a variable with many levels (like slide barcode) will soak up some variance by chance, and biology and batch are often confounded. If all samples of one class were run on two slides, the slide R² will be high for a component that is really about tumor type. The table tells you where to look, not what is true. A component driven mostly by array type or material, with class explaining little, is a red flag.

t-SNE

t-SNE (van der Maaten and Hinton, 2008) converts distances into neighbor probabilities in the high-dimensional space and then arranges points in 2D so that the 2D neighbor probabilities match. It is very good at keeping local neighborhoods together, which is why it separates methylation classes so nicely. It makes no promise about anything else.

from sklearn.manifold import TSNE

tsne = TSNE(
    n_components=2,
    perplexity=30,
    init="pca",
    learning_rate="auto",
    random_state=0,
)
emb_tsne = tsne.fit_transform(pcs)

Perplexity

Perplexity is roughly the effective number of neighbors each point pays attention to. Low values (5 to 10) emphasize very local structure and can break real classes into fragments. High values (50 to 100 or more) pull in more global structure and can merge small classes into their neighbors. It must be smaller than the number of samples, and it interacts with class size: if your smallest class has 8 samples, a perplexity of 100 asks each of those samples to consider many neighbors from other classes.

The honest approach is to run a small grid and look at all of them:

runs = {}
for perp in [5, 30, 100]:
    for seed in [0, 1, 2]:
        runs[(perp, seed)] = TSNE(
            perplexity=perp, init="pca", learning_rate="auto", random_state=seed
        ).fit_transform(pcs)

Features that appear across perplexities and seeds are probably real. Features that appear in only one run are probably not.

Seeds and initialization

t-SNE optimizes a non-convex objective, so different random seeds give different layouts. With init="pca" (the default in current scikit-learn) the global arrangement is more stable between runs than with random initialization, but it is not fixed. Always set random_state for anything you publish, and always look at more than one seed before you interpret anything.

What the distances mean

This is where most over-reading happens:

Cluster sizes mean nothing. t-SNE adapts to local density, so a tight class and a diffuse class can come out the same size. You cannot read within-class heterogeneity from cluster area.
Distances between clusters mean little. Two clusters being far apart in a t-SNE plot does not mean they are more different than two clusters that happen to sit close together. Global arrangement is weakly constrained.
Gaps are not evidence of discrete classes. t-SNE can make continuous variation look clumpy, especially at low perplexity.
What you can trust: points that are close together are usually close in the input space. A sample sitting in the middle of a known class is a reasonable first hint about what it resembles.

If a question depends on distances between classes, answer it in PCA space or on the original features, not on the t-SNE plot.

UMAP

UMAP (McInnes, Healy and Melville, 2018) builds a weighted k-nearest-neighbor graph in the input space and optimizes a 2D layout of that graph. In practice it looks similar to t-SNE, is faster on large cohorts, and tends to keep somewhat more of the global arrangement, although “somewhat more” is not the same as “faithful”.

import umap

reducer = umap.UMAP(
    n_neighbors=15,
    min_dist=0.1,
    metric="euclidean",
    random_state=0,
)
emb_umap = reducer.fit_transform(pcs)

The three parameters that matter:

n_neighbors plays a role similar to perplexity: small values emphasize local structure, large values emphasize the broader layout. The same caution about small classes applies.
min_dist controls how tightly points are packed in the embedding. Small values give dense, well-separated clumps; larger values spread points out. This is a display parameter: it changes how the picture looks, not which samples are neighbors.
metric sets the input distance. Euclidean on PCA scores is my default. Running UMAP directly on the top-probe beta matrix with metric="correlation" is a reasonable alternative and a useful comparison, because correlation distance focuses on the pattern across probes rather than the overall level.

reducer_corr = umap.UMAP(n_neighbors=15, min_dist=0.1,
                         metric="correlation", random_state=0)
emb_umap_corr = reducer_corr.fit_transform(B)

On reproducibility: setting random_state in umap-learn makes results reproducible but disables most of its parallelism, so it runs slower. umap-learn will warn you about this. Accept the warning for anything you plan to show someone.

UMAP has one practical advantage for methylation work: a fitted reducer has a transform method, so you can place new samples onto an existing embedding (reducer.transform(new_pcs)) as long as you apply the same probe list and the same fitted PCA. That is handy for showing where a new case lands relative to a reference cohort. It is still a visualization, not a classifier.

Checking stability

I like at least one number to go with my eyes. Two simple ones:

Neighbor overlap between runs. For each sample, compare its k nearest neighbors in one embedding with its k nearest neighbors in another (a different seed, or the input PCA space). The average overlap tells you how consistent the local structure is.

from sklearn.neighbors import NearestNeighbors
from sklearn.manifold import trustworthiness

def knn_index(Z, k=15):
    nn = NearestNeighbors(n_neighbors=k + 1).fit(Z)
    return nn.kneighbors(Z, return_distance=False)[:, 1:]   # drop self

def knn_overlap(Z1, Z2, k=15):
    a, b = knn_index(Z1, k), knn_index(Z2, k)
    return np.mean([len(set(r1) & set(r2)) / k for r1, r2 in zip(a, b)])

print("t-SNE seed 0 vs seed 1:", round(knn_overlap(runs[(30, 0)], runs[(30, 1)]), 2))
print("t-SNE vs PCA input:    ", round(knn_overlap(runs[(30, 0)], pcs), 2))

Trustworthiness. scikit-learn’s trustworthiness measures how many of the neighbors in the embedding are also neighbors in the input space, penalizing “false” neighbors. Values close to 1 mean the embedding rarely puts unrelated samples next to each other.

print("trustworthiness:", round(trustworthiness(pcs, emb_tsne, n_neighbors=15), 3))

Neither number says the plot is “correct”. They tell you whether the local structure you are looking at is consistent and grounded in the input, which is the part of the plot that deserves your trust anyway.

Color by everything, not just class

The single most useful habit is to make the same embedding several times, colored by different variables: class, then array type, slide, center, material, sex, and anything else in the sample sheet. If samples group by slide within a class, or EPIC samples sit on one side of every cluster, you have technical structure that will also affect any downstream analysis.

import matplotlib.pyplot as plt

def scatter_by(emb, labels, ax, title, max_legend=20):
    labels = pd.Series(labels).astype(str).to_numpy()
    cats = pd.unique(labels)
    cmap = plt.get_cmap("tab20")
    for i, c in enumerate(cats):
        sel = labels == c
        ax.scatter(emb[sel, 0], emb[sel, 1], s=4, color=cmap(i % 20),
                   label=c if len(cats) <= max_legend else None, linewidths=0)
    ax.set_title(title)
    ax.set_xticks([])
    ax.set_yticks([])
    if len(cats) <= max_legend:
        ax.legend(markerscale=3, fontsize=6, frameon=False)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
scatter_by(emb_umap, meta["methylation_class"], axes[0], "UMAP, colored by class")
scatter_by(emb_umap, meta["array_type"], axes[1], "UMAP, colored by array type")
scatter_by(emb_umap, meta["slide"], axes[2], "UMAP, colored by slide")
fig.tight_layout()
fig.savefig("umap_checks.png", dpi=150)

With 80 or more classes, no color palette will make every class distinguishable, and the function above simply drops the legend. Two workarounds: color by class family (a coarser grouping) and label cluster centers with text, or highlight one class at a time in color against grey for everything else. The second is surprisingly effective for checking whether a specific class is coherent.

I also drop the axis ticks on t-SNE and UMAP plots. The axis values carry no meaning, and leaving them off discourages people from reading them. (PCA axes do mean something, so keep them there, ideally with the variance explained in the axis label.)

The same thing in R

If your preprocessing is in R anyway, prcomp, Rtsne and uwot cover the same ground. Here betas is the probes by samples matrix from minfi or sesame.

library(Rtsne)
library(uwot)

betas <- betas[complete.cases(betas), ]
sds   <- apply(betas, 1, sd)
top   <- names(sort(sds, decreasing = TRUE))[1:10000]
x     <- t(betas[top, ])                      # samples x probes

pc <- prcomp(x, center = TRUE, scale. = FALSE, rank. = 50)
var_explained <- pc$sdev^2 / sum(pc$sdev^2)
round(var_explained[1:10], 3)

set.seed(1)
ts <- Rtsne(pc$x, pca = FALSE, perplexity = 30, theta = 0.5,
            check_duplicates = FALSE)
emb_tsne <- ts$Y

set.seed(1)
emb_umap <- uwot::umap(pc$x, n_neighbors = 15, min_dist = 0.1,
                       metric = "euclidean")

pca = FALSE in Rtsne matters because the input is already principal components; by default Rtsne would run its own PCA. matrixStats::rowSds is a much faster replacement for the apply call on a full array. For uwot, set.seed gives reproducible results as long as the optimization runs single-threaded, which is the default for the stochastic gradient descent step; check the n_sgd_threads documentation if you change thread settings. I call uwot::umap explicitly because the separate umap package exports a function with the same name.

Probe QC done, sex chromosome probes removed, probe selection documented (how many, by which statistic, on betas or M-values).
PCA run first; variance explained reported; top components checked against batch, array type, slide and material.
t-SNE and UMAP run on PCA scores with parameters stated in the figure legend or methods.
At least three seeds and two or three perplexity or n_neighbors values looked at; the conclusions hold across them.
Same embedding colored by class and by each technical variable.
No claims based on cluster size or distance between clusters in t-SNE, and only cautious ones in UMAP.
Seeds, package versions, probe list and PCA loadings saved, so the figure can be regenerated and new samples projected onto it.

If you want to go from “this sample lands near class X” to an actual prediction, that is a classification problem with its own rules about validation (feature selection inside cross-validation, a locked test set, calibrated scores), and a topic for its own post.

References

Capper D, et al. DNA methylation-based classification of central nervous system tumours. Nature 555, 469–474 (2018).
Du P, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics 11, 587 (2010).
van der Maaten L, Hinton G. Visualizing data using t-SNE. Journal of Machine Learning Research 9, 2579–2605 (2008).
McInnes L, Healy J, Melville J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03426 (2018).
Aryee MJ, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363–1369 (2014).

← All writing