Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Uncertainty in Gene-Level Velocity Estimates

RNA Velocity · Tutorial 397 of 700

Uncertainty in Gene-Level Velocity Estimates

Use parameter variability, fit quality, and changes across cell subsets to judge how stable a gene-level velocity estimate is.

Advanced 12 min read

What You'll Learn

  • Distinguish uncertainty in kinetic rates from uncertainty in a cell’s estimated velocity.
  • Use the inverse observed information to approximate parameter uncertainty.
  • Propagate rate uncertainty to velocity with the delta method.
  • Use cell-level bootstrap refits to assess variability without relying only on local curvature.
  • Check how conclusions change when the cells used for fitting are varied.
  • Recognize when fit quality, weak identifiability, or cell selection makes uncertainty summaries unreliable.

Why Uncertainty Matters for a Gene-Level Velocity

A fitted velocity is an estimate, not a directly observed quantity. For a gene, it depends on kinetic rates estimated from cell counts and on the cell’s observed or estimated expression state. Two fits can produce similar rates but meaningfully different velocities, especially when the terms describing production and degradation nearly cancel. A numerical estimate without an uncertainty assessment can therefore suggest more precision than the data support.

The focus here is uncertainty in the kinetic parameters and in the resulting gene-level velocity. This is distinct from deciding which model to use, and from interpreting or combining velocity vectors. The introductory material on likelihood-based kinetic model fitting provides the setting: parameters are estimated by fitting a kinetic model to the observed counts.

Definition: For one gene, let \(u\) denote unspliced abundance and \(s\) spliced abundance in a cell. In the two-compartment model, the instantaneous changes are $$ \frac{du}{dt}=\alpha-\beta u,\qquad \frac{ds}{dt}=\beta u-\gamma s, $$ where \(\alpha\) is the transcription rate, \(\beta\) the splicing rate, and \(\gamma\) the spliced-RNA degradation rate. The gene-level spliced velocity at state \((u,s)\) is \(v=\beta u-\gamma s\). Uncertainty in \(v\) means uncertainty in this estimated rate of change, conditional on the model and the stated cell state.

Rates and velocities have units tied to the time scale used by the model. For example, if abundance is measured in normalized counts and time in hours, \(\beta\) and \(\gamma\) have units of inverse hours, and \(v\) has units of normalized counts per hour. A velocity is not itself a probability, nor is it a guarantee that a cell will follow a particular future path.

What an Uncertainty Assessment Assumes

Any standard uncertainty calculation is conditional on assumptions. For likelihood-based fitting, suppose the observations used for a gene are described by a correctly specified likelihood \(L(\theta)\), where \(\theta\) collects the fitted kinetic parameters. The likelihood must account adequately for the observation process; observations are treated as independent conditional on the model, unless their dependence is explicitly represented. The fit should converge to a meaningful interior solution, and the parameters should be identifiable from the data. In particular, different parameter values should not produce essentially indistinguishable predictions over the observed cells.

These assumptions are substantive. A narrow uncertainty estimate from a misspecified model can be misleading: it measures precision under the assumed likelihood, not whether the kinetic model is biologically correct. Likewise, a parameter at a boundary, a nearly flat likelihood, or a highly correlated set of rates can make ordinary local approximations unreliable.

Three related questions are useful, but they are not interchangeable:

  • Fit quality: do fitted values describe the observations adequately, without conspicuous systematic discrepancies?
  • Parameter variability: how much might the fitted kinetic rates vary under repeated data collection or resampling?
  • Cell-inclusion sensitivity: does the estimate change when the cells contributing to the fit are changed?

Good fit quality does not prove low parameter uncertainty. Conversely, a parameter estimate can be numerically stable while the model fits poorly. Each question needs its own evidence.

Local Uncertainty from Likelihood Curvature

When the likelihood has a well-defined, approximately quadratic maximum, its curvature gives a local approximation to parameter uncertainty. Let \(\hat{\theta}\) be the maximum-likelihood estimate and let \(\ell(\theta)=\log L(\theta)\). The observed information is the negative matrix of second derivatives of the log-likelihood, evaluated at the estimate. Under regularity conditions and a sufficiently informative sample, its inverse approximates the covariance matrix of \(\hat{\theta}\).

$$ \widehat{\operatorname{Cov}}(\hat{\theta})\approx \mathcal{I}(\hat{\theta})^{-1},\qquad \mathcal{I}(\hat{\theta})= -\left.\frac{\partial^2\ell(\theta)} {\partial\theta\,\partial\theta^\mathsf{T}}\right|_{\theta=\hat{\theta}}. $$

The diagonal entries approximate parameter variances; off-diagonal entries describe estimated covariation between parameters. Those covariances matter when estimating velocity. For fixed state values \(u\) and \(s\), the velocity depends on \(\beta\) and \(\gamma\), so a first-order approximation propagates their joint uncertainty through \(v=\beta u-\gamma s\). This is the delta method: approximate the function locally by its linear expansion around the fitted parameters.

$$ \widehat{\operatorname{Var}}(\hat{v})\approx \nabla v(\hat{\theta})^\mathsf{T} \widehat{\operatorname{Cov}}(\hat{\theta}) \nabla v(\hat{\theta}),\qquad \nabla_{(\beta,\gamma)}v=(u,-s). $$

If \(u\) and \(s\) are treated as fixed, the gradient has only the two displayed components. If the cell state is also estimated with uncertainty, its uncertainty and covariance with the fitted rates should be included in a joint calculation; treating it as fixed then understates total uncertainty. A local interval of the form estimate plus or minus 1.96 standard errors is an approximate 95% interval when the sampling distribution is adequately close to normal. It is not reliable automatically, especially near parameter boundaries or with a flat, asymmetric likelihood.

Worked Example: Propagating Rate Uncertainty to Velocity

Suppose a hypothetical gene fit gives \(\hat{\beta}=0.400\) per hour and \(\hat{\gamma}=0.200\) per hour. At a cell state treated as fixed, \(u=2.00\) normalized counts and \(s=6.00\) normalized counts. The estimated covariance matrix for the two rates is $$ \widehat{\operatorname{Cov}}(\hat{\beta},\hat{\gamma})= \begin{pmatrix} 0.000400 & -0.000100\\ -0.000100 & 0.000100 \end{pmatrix}. $$ Assume the likelihood fit is converged, the rates are identifiable in the fitted data, and a local normal approximation is reasonable.

The fitted velocity is $$ \hat{v}=(0.400)(2.00)-(0.200)(6.00)=0.800-1.200=-0.400 $$ normalized counts per hour. The gradient is \((2.00,-6.00)\). Thus $$ \widehat{\operatorname{Var}}(\hat{v}) =(2.00)^2(0.000400)+(-6.00)^2(0.000100) +2(2.00)(-6.00)(-0.000100)=0.007600. $$ The standard error is \(\sqrt{0.007600}=0.08718\) normalized counts per hour. An approximate 95% interval is $$ -0.400\mathbin{\pm}1.96(0.08718) =-0.400\mathbin{\pm}0.1709 =[-0.5709,-0.2291] $$ normalized counts per hour. Under the stated approximation, this interval indicates uncertainty in the fitted velocity conditional on this cell state and model. It does not account for uncertainty in \(u\) and \(s\), or for model misspecification.

The covariance term contributes \(0.002400\) to the variance. Omitting it would give a different standard error, because uncertainty in \(\hat{\beta}\) and \(\hat{\gamma}\) is related rather than independent.

Resampling Cells and Checking Sensitivity

A cell-level bootstrap offers a different view of variability. Repeatedly sample cells with replacement from the cells used for fitting, refit the kinetic parameters in each resample, and calculate the velocity for each fit. The spread of the resulting estimates reflects how sensitive the fit is to the sampled cells, under the assumption that the observed cells reasonably represent the population of interest. Each resample must repeat the fitting procedure; holding the rates fixed would not measure their fitting variability.

The bootstrap is not a cure for a bad likelihood or unrepresentative data. If cells are clustered by donor, batch, or another sampling unit, resampling individual cells as if independent can give an inappropriate uncertainty estimate. The resampling unit should reflect the design; for example, independent donors may need to be resampled as units. If cells were deliberately sampled across known states, stratified resampling may better preserve that design. A small number of bootstrap refits is useful for illustration but not for a stable interval.

Worked Example: Summarizing Bootstrap Rate and Velocity Variability

Suppose eight hypothetical cell-level bootstrap refits for one gene produce the following rates, in inverse hours:

Refit\(\hat{\beta}\)\(\hat{\gamma}\)Velocity at \(u=2.00,\ s=6.00\)
10.360.18-0.36
20.390.20-0.42
30.410.21-0.44
40.400.19-0.34
50.440.22-0.44
60.380.19-0.38
70.420.20-0.36
80.400.21-0.46

For refit 1, for example, \(v=(0.36)(2.00)-(0.18)(6.00)=0.72-1.08=-0.36\). The mean of the eight \(\hat{\beta}\) values is \(3.20/8=0.400\) per hour. Their squared deviations sum to \(0.004200\), so their sample standard deviation is \(\sqrt{0.004200/7}=0.02449\) per hour. The mean \(\hat{\gamma}\) is \(1.60/8=0.200\) per hour; its squared deviations sum to \(0.001200\), giving a sample standard deviation of \(\sqrt{0.001200/7}=0.01309\) per hour.

The eight velocities sum to \(-3.20\), so their mean is \(-0.400\) normalized counts per hour. Their squared deviations from that mean sum to \(0.01440\), giving a sample standard deviation of \(\sqrt{0.01440/7}=0.04536\) normalized counts per hour. All eight estimates are negative, so their estimated direction is stable across these resamples. With only eight refits, however, this is a descriptive illustration, not a dependable 95% bootstrap interval. In an analysis, many more refits would be used, and the full bootstrap distribution would be inspected for skewness or multiple modes.

A complementary check is cell-inclusion sensitivity. Refit after excluding a cell, a group of cells, or a specified subset, and compare the rates and resulting velocities. This is an influence analysis, not a confidence interval: the selected subsets may not represent repeated random samples. It is especially informative when a few cells dominate the phase portrait or when estimates depend strongly on a narrow range of observed states.

Worked Example: Sensitivity to Which Cells Are Included

Consider a hypothetical gene evaluated at the same fixed state, \(u=2.00\) and \(s=6.00\). The full-cell fit gives \(\hat{\beta}=0.400\) and \(\hat{\gamma}=0.200\) per hour, so \(\hat{v}=0.800-1.200=-0.400\) normalized counts per hour. Refitting after retaining one prespecified group of cells gives \(\hat{\beta}=0.420\) and \(\hat{\gamma}=0.190\), while another prespecified subset gives \(\hat{\beta}=0.360\) and \(\hat{\gamma}=0.220\).

For the first subset, \(\hat{v}=(0.420)(2.00)-(0.190)(6.00)=0.840-1.140=-0.300\) normalized counts per hour. For the second, \(\hat{v}=(0.360)(2.00)-(0.220)(6.00)=0.720-1.320=-0.600\) normalized counts per hour. The range across these three fits is \([-0.600,-0.300]\); the sign is unchanged, but the magnitude varies by a factor of two from the least to the greatest absolute estimate.

This suggests that direction is more stable than magnitude under these particular changes in cell inclusion. It does not show that the sign is certain in a statistical sense, and it does not identify which fit is correct. The subsets should be defined for a scientific or design-based reason, not chosen after inspecting results to obtain a preferred velocity.

Reading the Evidence and Avoiding Overconfidence

The three approaches answer different questions. Likelihood curvature is computationally efficient and uses local information at the fitted optimum, but depends on a good local quadratic approximation. A bootstrap can reveal asymmetry and instability across resamples, but depends on an appropriate resampling scheme and repeated successful fits. Cell-inclusion sensitivity is direct and easy to interpret, but depends on which alternative subsets are examined and is not a sampling interval.

ApproachWhat it summarizesMain caution
Observed-information approximationLocal parameter and velocity uncertainty near the fitted optimumUnreliable with weak identification, boundaries, or a non-quadratic likelihood
Cell-level bootstrapVariability across resampled cell sets and refitsResampling must reflect dependence and sampling design
Cell-inclusion sensitivityChange under specified alternative sets of included cellsNot a confidence interval; subsets may be selective or unrepresentative

Fit quality should be considered alongside these summaries. Check that optimization converged, that the likelihood surface is not nearly flat in important parameter directions, and that fitted values do not show obvious systematic departures from the observed pattern. Large residual structure or many failed refits weakens the interpretation of a narrow standard error. These checks flag reasons to distrust a precision claim; they do not by themselves repair model misspecification.

Common mistake: “The fitted velocity is negative, so the cell is definitely moving toward lower spliced expression.” This treats a point estimate as certain and interprets an instantaneous model-based rate as a guaranteed future outcome. The correction is to report the estimated velocity with units, assess its uncertainty and sensitivity, and describe it as a model-based instantaneous rate. A confidence interval that includes zero would not establish a reliably negative rate under the approximation used; an interval excluding zero still does not prove the model is correct.

A responsible report gives the fitted rate or velocity, its units, how uncertainty was assessed, and the relevant qualifications: fixed versus uncertain cell state, resampling unit, convergence issues, and sensitivity to cell inclusion. A small standard error is not evidence that these qualifications can be ignored.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

  1. For \(u=3\), \(s=4\), \(\hat{\beta}=0.30\) per hour, and \(\hat{\gamma}=0.15\) per hour, calculate the velocity and its units if abundance is in normalized counts.
  2. Why must the covariance between fitted rates be considered when propagating uncertainty to velocity?
  3. A bootstrap gives velocities that are mostly negative but include several positive values. What does this say, and what does it not say?
  4. How does a cell-inclusion sensitivity analysis differ from a bootstrap confidence interval?

Solutions. (1) \(v=(0.30)(3)-(0.15)(4)=0.90-0.60=0.300\) normalized counts per hour. (2) The velocity depends on both rates, and their estimation errors may covary; the variance formula includes a covariance term, so assuming independence can misstate uncertainty. (3) The sign changes across resampled fits, indicating that the estimated direction is sensitive to resampling. This alone does not give a calibrated probability that the true velocity is positive or negative; the bootstrap design and distribution must be assessed. (4) A bootstrap describes variability across resamples intended to represent repeated sampling under a specified design. Sensitivity analysis compares selected alternative cell sets and is a robustness check, not a confidence interval.

Clinical Biostats