Why Apparent Clusters May Not Be Separate Populations
A clustering algorithm must divide cells into groups, even when the data vary gradually along a trajectory. For example, cells may change smoothly from one transcriptional state to another as a differentiation program progresses. A partition into two or three clusters can summarize regions of that trajectory, but the partition alone does not show that the cells come from distinct populations.
The question here is narrower than whether a clustering is useful or reproducible: does a cluster boundary coincide with an abrupt change that is difficult to explain as part of a smooth trend? Answering it requires specifying what “smooth” means, choosing a representation in which that assumption is plausible, and calibrating the comparison against a continuous-structure model. There is no assumption-free test that can decide whether clusters are biologically “real.”
These are working statistical models, not claims about every biological system. Real data can contain both: a smooth trajectory within populations, discrete transitions between them, or branching and other structures that neither simple model captures. The goal is to test a specified continuous explanation against a specified boundary diagnostic, not to force every data set into a binary choice.
A Boundary Diagnostic Under a Smooth-Trend Model
Suppose \(u_i\) is a latent or estimated coordinate for cell \(i\), and \(x_{ig}\) is a suitably transformed expression value for gene \(g\). A simple continuous model is
where \(s_g\) is a smooth gene-specific trend and \(\epsilon_{ig}\) is residual variation. “Smooth” means that nearby values of \(u\) are expected to have similar fitted expression; it does not require a straight line. The model does not assume that all genes change in the same direction, or that every cell lies exactly on a curve.
Fit the trends and calculate residuals \(e_{ig}=x_{ig}-\hat{s}_g(u_i)\). For each adjacent pair of cells \(i,j\) in the coordinate ordering, define their mean squared residual contrast across a fixed set of genes:
Let \(B\) be adjacent pairs whose cells have different cluster labels, and let \(W\) be nearby adjacent pairs with the same label. A simple diagnostic is the ratio
A large \(R\) indicates that residual differences are elevated at the cluster boundaries relative to local within-cluster differences. A value near one does not establish continuity; it means this particular boundary diagnostic did not find a pronounced excess. If \(B\) or \(W\) has very few pairs, the ratio is unstable and should not be treated as persuasive evidence.
Calibration and Assumptions
A useful simulation test starts with a null model that represents the continuous explanation being considered. In single-cell data, one possible model fits smooth gene-specific trends while accounting for library size and an appropriate count or transformed-expression noise model. Simulated data should preserve important features of the analysis, such as the number of cells, gene-wise variability, and relevant dependence or overdispersion. The method used to generate simulations must be stated; a smooth curve with independent Gaussian errors, for instance, is not automatically an adequate model for raw counts.
For each simulated data set, repeat the steps that could affect the statistic: preprocessing, estimating the coordinate if it is data-derived, fitting trends, clustering, and calculating \(R\). If the observed cluster number or resolution was selected from the data, that selection should also be repeated when feasible. This matters because a boundary selected for looking unusually sharp will tend to look sharp even under a continuous model. If the coordinate and cluster labels are fixed independently of the expression features used to calculate residuals, the diagnostic is easier to interpret; otherwise, simulation should reproduce their dependence.
With \(M\) simulations and \(K\) simulated statistics at least as large as the observed value, the upper-tail Monte Carlo p-value is
The added one prevents a reported value of zero and gives a finite-simulation estimate. This is a calibrated test only to the extent that the simulation model is credible. If model parameters are estimated, the result is generally a parametric-bootstrap assessment, not an exact test with guaranteed error rates under every data-generating process.
Specify the coordinate, expression scale, smooth trends, and noise assumptions. Decide which genes will contribute to the diagnostic.
Fit the smooth trends, calculate \(R\), and inspect whether boundary pairs are sufficiently numerous and comparable to local within-cluster pairs.
Generate data under the null, then repeat the analysis steps that influence cluster labels, the coordinate, and \(R\).
Compare \(R\) with its simulated reference distribution, examine model fit, and report the result as evidence about that null model rather than a definitive biological classification.
Worked Examples
Worked Example: A Cluster Boundary Along a Smooth Gradient
Suppose a hypothetical experiment has cells ordered by an independently specified coordinate \(u\). After fitting smooth expression trends across a fixed gene set, the mean squared residual contrast for the one adjacent cross-cluster pair is \(0.010\). The three local within-cluster contrasts are \(0.040\), \(0.010\), and \(0.040\). Calculate \(R\).
The mean within-cluster contrast is
Therefore,
The boundary residual contrast is smaller than the average local within-cluster contrast. In this hypothetical data set, the statistic does not suggest an abrupt residual jump at the boundary. Suppose 310 of \(M=1{,}000\) continuous-model simulations produce \(R\) at least as large as the observed value. Then
This is not evidence against the specified smooth model at conventional small significance levels. It does not prove that the cells lie on a continuum: the diagnostic could have limited power, and its conclusion depends on how well the simulated null represents the data.
Worked Example: A Boundary With a Large Residual Jump
Now consider a hypothetical data set with a single cross-cluster pair having residual contrast \(0.810\), while four nearby within-cluster pairs have contrasts \(0.010\), \(0.020\), \(0.010\), and \(0.040\). The mean within-cluster contrast is
Thus,
The residual difference across the boundary is much larger than the local within-cluster differences. Suppose 2 of \(1{,}000\) simulations from the fitted continuous model yield \(R\geq40.5\). The Monte Carlo p-value is
This is evidence that the observed boundary is unusual under the specified smooth null. It remains conditional on adequate model fit and faithful simulation. A misspecified trend, an unmodeled technical batch effect, or an expression feature used both to define the cluster and measure its boundary could also produce a large statistic.
Worked Example: Clusters That Divide a Single-Cell Trajectory
Imagine five hypothetical cells ordered by an independently measured developmental coordinate. A marker score, in transformed expression units, is \(1.0, 1.5, 2.0, 2.5, 3.0\) at coordinates \(0, 0.25, 0.50, 0.75, 1.00\), respectively. A clustering assigns the first two cells to A, the next two to B, and the last cell to C.
The score follows the exact linear trend \(\hat{s}(u)=1+2u\). For example, at \(u=0.50\), \(\hat{s}(0.50)=1+2(0.50)=2.0\), matching the observed score. The residual for every cell is zero, so each adjacent residual contrast is \(0\). In this idealized example, the raw score increases by \(0.5\) at every step, including across both label boundaries; there is no boundary-specific change in the score.
The cluster labels therefore summarize successive regions of a smooth trend in this marker. They do not, by themselves, establish three separated populations. In real data, residuals will not be exactly zero, so the analyst would use several genes, estimate smooth trends, and compare boundary residual contrasts with a simulation-calibrated reference. If the developmental coordinate was estimated from the same expression data, simulations should repeat that estimation, or an independent coordinate or held-out gene set should be considered.
Diagnostics, Alternatives, and Limitations
A numerical result is more useful when paired with checks on whether the diagnostic's assumptions are plausible. Plot expression against the coordinate for representative genes, along with fitted smooth trends and residuals. Check whether residual spread changes substantially across the coordinate, whether a few genes dominate the aggregate statistic, and whether boundary pairs are much farther apart in the coordinate than comparison pairs. Inspect the number and location of cross-label pairs: one boundary pair cannot support a broad claim about a population-wide separation. Repeat the calculation with defensible smoothing choices and gene sets, and report sensitivity rather than selecting only the result that favors a preferred interpretation.
A poorly estimated coordinate can hide a real discontinuity or manufacture apparent smoothness. Branching trajectories also cannot generally be represented by a single ordering. Batch effects, donor composition, library-size artifacts, and preprocessing choices can create or obscure apparent boundaries. The earlier material on preprocessing sensitivity and distance geometry is relevant: a diagnostic depends on the representation and distances used, so its inputs should be described and checked rather than treated as neutral.
Other validity tools address different questions. The gap statistic compares within-cluster dispersion with a chosen reference distribution; silhouette summaries compare within- and between-cluster distances; stability analyses ask whether assignments persist under specified perturbations. A partition can be stable and still divide a continuous gradient, while a boundary diagnostic can find no excess jump even when the partition is a useful summary. None of these results alone determines the biological meaning of a cluster. The later comparison of validity metrics considers what to do when such diagnostics point in different directions.
Check Your Understanding
Answer each question using the specified model and the limits of the diagnostic.
- If the mean boundary residual contrast is \(0.06\) and the mean matched within-cluster contrast is \(0.03\), calculate \(R\). What does the value establish?
- In 499 simulations, 14 statistics are at least as large as the observed statistic. Calculate the Monte Carlo p-value.
- Why should a data-derived coordinate be re-estimated in each simulated data set when calibrating the diagnostic?
- Suppose the observed boundary statistic is not unusually large under a continuous-model simulation. Can one conclude that the clusters are biologically identical? Explain.
Solutions. (1) \(R=0.06/0.03=2.00\). The boundary residual contrast is twice the matched within-cluster contrast, but this alone does not establish a discrete population boundary; its reference distribution is needed. (2) \(\hat{p}_{\mathrm{MC}}=(14+1)/(499+1)=15/500=0.0300\). (3) The estimated coordinate affects fitted trends, residuals, and potentially which pairs are adjacent. Holding it fixed would omit that source of variability from the simulated reference distribution. (4) No. The result says only that this statistic did not provide strong evidence against the specified continuous model; low power, model misspecification, or other biological structure remain possible.