Resampling Stability of Inferred Pseudotime

Why Assess Resampling Stability?

A trajectory method returns a path and an ordering of cells from a particular data set, but those outputs can depend on which cells were observed. A path may be determined by a small number of cells in a sparse region; a slightly different sample could connect the same broad regions in another way. Likewise, a cell’s position may shift substantially even when the overall path looks similar. Resampling stability asks how much these inferred structures change when the observed cells are perturbed by a specified sampling scheme.

This is a robustness analysis, not a test that establishes a unique biological trajectory. It is useful only after defining what is being compared: the path, the relative ordering of cells, or both. The comparison must also account for arbitrary pseudotime scales and for the possibility that a fitted path changes direction or branch structure.

Specify the Resampling and Refitting Procedure

Let the observed cells be \(1,\ldots,n\), and let a trajectory-fitting pipeline return a path and a pseudotime \(t_i\) for cells it represents. A resampling analysis creates data sets indexed by \(b=1,\ldots,B\), fits the pipeline to each, and compares the resulting outputs. In a subsampling analysis, each data set contains a specified fraction of the cells sampled without replacement. In a cell bootstrap, \(n\) cells are drawn with replacement; some original cells therefore occur multiple times and some do not occur in a particular draw.

State the sampling unit and what is held fixed. If the question concerns sensitivity to which cells were observed, resample cells and repeat the data-dependent steps that define the trajectory, including graph construction, path fitting, and root selection when root selection is part of the intended pipeline. If the root is fixed by design, keep the same root rule in every fit. A different choice answers a different stability question. Similarly, if preprocessing is treated as part of the method, refit it in each resample; if it is deliberately held fixed, report that the result is conditional on that preprocessing.

For subsamples, the fitted trajectory may not assign pseudotime to omitted cells. Compare ordering on the cells shared by the fits, or use a prespecified projection rule to place omitted cells onto each fitted trajectory. Do not silently treat unassigned cells as if they had unchanged positions. For bootstrap fits, track original cell identities and state how duplicate draws are handled; counting a duplicated cell as multiple independent identities can distort comparisons.

Comparing Cell Orderings

Raw pseudotime values are often unsuitable for direct comparison. Their numerical scale can change across fits even when the order is essentially the same. A rank comparison removes this scale difference. For a fit \(b\) on \(m\) shared cells, let \(r_i^{(b)}\) be cell \(i\)'s rank, using average ranks for ties. Define its normalized rank by

For \(m>1\), the normalized ranks run from 0 to 1 when there are no ties. They describe relative position, not time in hours or a rate of biological change. Before comparing ranks, orient each fit consistently—for example, by using the same root rule. If orientation cannot be aligned, a reversed ordering can look maximally different even when it describes the same unoriented path.

One summary is Spearman's rank correlation, the ordinary correlation of the two rank vectors on the shared cells. With no ties, it can be calculated from rank differences \(d_i=r_i^{(b)}-r_i^{(\mathrm{ref})}\):

This formula follows from the correlation of the two rank vectors. Each vector has mean \((m+1)/2\) and sum of squared deviations \(m(m^2-1)/12\). Since both vectors have the same sum of squares, expanding \(\sum_i(r_i^{(b)}-r_i^{(\mathrm{ref})})^2\) gives \(2\sum_i(r_i-\bar r)^2(1-\rho_s)\). Substituting the rank sum of squares and rearranging yields the displayed formula. With ties, calculate the correlation of average ranks directly rather than applying the no-ties shortcut.

Correlation summarizes broad agreement but can conceal important local changes. Also report a measure such as the mean absolute normalized-rank displacement,

A smaller \(D^{(b)}\) means less average movement in relative position. Pairwise order agreement is another interpretable option: among cell pairs, count the proportion ordered the same way in both fits, with a pair tied in one fit contributing one-half. These summaries should be calculated over a clearly specified set of shared cells.

Comparing Inferred Paths

For a cell-level tree or graph, a direct summary is the frequency with which each edge appears across fitted paths. If \(I_e^{(b)}=1\) when edge \(e\) appears in resample \(b\), and \(0\) otherwise, its empirical edge frequency is

A value near 1 indicates that the edge recurs under the chosen resampling and fitting procedure; a value near 0 indicates that it rarely does. This interpretation requires edges to be comparable across fits. For cell-level edges, use consistent cell identities and specify how missing cells are handled. For geometric paths or skeletons represented differently across fits, define a matching rule before calculating frequencies; otherwise, apparent differences may reflect representation rather than a changed trajectory.

A set-level comparison can use Jaccard similarity. For edge sets \(E_a\) and \(E_b\) on a comparable representation, it is the size of their intersection divided by the size of their union. It equals 1 for identical nonempty sets and 0 for disjoint sets. Report how pairwise similarities are aggregated—for example, the mean similarity between each resampled path and a reference path—because a single average can conceal a small group of very different fits.

Interpret the Distribution, Not Just One Number

With many resamples, summarize the distribution of ordering measures and path features. For each cell, the spread of its normalized ranks across fits can identify cells whose relative positions are uncertain under the specified perturbation. For edges, display frequencies or classify edges using a stated threshold chosen for the analysis. There is no universal cutoff that turns an edge into a true or false biological connection. Preserve and report failed, disconnected, or otherwise unusable fits: their frequency is itself evidence that the pipeline is sensitive to the resampling.

Assumptions, Diagnostics, and Limitations

A cell-level bootstrap treats the observed cells as exchangeable draws from a population relevant to the analysis. That assumption may be implausible when cells are clustered within a small number of donors, samples, or batches. Resampling cells within one data set then assesses sensitivity to cell composition conditional on those data; it does not measure uncertainty across biological replicates. To assess between-replicate variation, the resampling unit must reflect the experimental design, such as resampling independent biological samples when that is justified and feasible.

Inspect whether instability is concentrated in particular regions or reflects a global change. Useful diagnostics include root changes, orientation reversals, branch creation or loss, disconnected fits, cells absent from many subsamples, and large rank displacements for particular cells. Check that comparisons use enough shared cells and that any projection or path-matching rule is applied consistently. Sensitivity to resampling fraction can also matter: very small subsamples may break a graph simply because too few cells remain, while very large subsamples can make resampled fits look similar.

Resampling is not the only robustness check. Parameter sensitivity analysis varies choices such as neighborhood or kernel scales while holding the observed cells fixed; it answers how the output responds to modeling choices, not which cells were observed. Comparing principal curves, MSTs, and diffusion maps assesses dependence on trajectory representation or method. These analyses complement one another and should not be conflated. A stable output can still be systematically wrong because of unsuitable measurements, confounding, or an inappropriate trajectory assumption. Conversely, instability may reflect limited sampling rather than absence of an underlying biological process.

Common Mistake

Wrong reasoning: “The edge appeared in 90% of bootstrap fits, so there is a 90% probability that this biological transition is real.” The edge frequency is not a posterior probability, a p-value, or a calibrated probability that a biological claim is true. It is the proportion of the specified resampled fits in which that edge was recovered, under the chosen data and pipeline.

Correction: Report the resampling scheme, refitting choices, comparison rule, and observed stability summaries. Interpret them as conditional robustness evidence. Do not use them alone to claim biological truth, causal direction, or uncertainty across donors.