Chapter 10: RNA Velocity Dynamics with Velocyto & scVelo
Chapter 10: RNA Velocity Dynamics with Velocyto & scVelo
While standard single-cell analysis captures a static snapshot of gene expression, RNA velocity leverages the kinetic ratio between unspliced pre-mRNA (containing introns) and mature spliced mRNA (exons only) to forecast the future transcriptional state of individual cells over a timescale of hours.
Transcription, splicing, and degradation are modeled by two coupled ordinary differential equations:
rac{du}{dt} = lpha - eta u, \quad rac{ds}{dt} = eta u - \gamma s
- \(lpha\): Transcription rate (induction of nascent pre-mRNA)
- \(eta\): Splicing rate converting unspliced (\(u\)) to spliced (\(s\)) mRNA
- \(\gamma\): Degradation rate of mature spliced mRNA
When a gene is actively induced, unspliced pre-mRNA rises first; when repressed, unspliced transcripts drop before mature spliced mRNA degrades.
1. Quantifying Spliced and Unspliced Counts (velocyto)
Standard BAM alignment files must be processed with velocyto to generate loom matrices containing separate layers for spliced, unspliced, and spanning reads.
# Run velocyto CLI on 10x BAM file
velocyto run10x \
-m hg38_rmsk.gtf \ # RepeatMasker GTF to exclude noisy repetitive introns
/data/pbmc_10x_run \
/ref/genes.gtf
# Produces: /data/pbmc_10x_run/velocyto/pbmc_10x_run.loom2. Dynamical RNA Velocity Modeling with scVelo
The scVelo toolkit in Python extends velocity modeling beyond steady-state approximations to solve full transcriptional kinetic phase portraits.
import scvelo as scv
import scanpy as sc
# Load loom file with spliced/unspliced matrices
adata = scv.read("pbmc_10x_run.loom", cache=True)
# Preprocessing: filter and normalize spliced and unspliced layers
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
# Estimate RNA velocity using the dynamical model
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
# Project velocity stream vectors onto UMAP embedding
scv.pl.velocity_embedding_stream(
adata,
basis='umap',
color='cell_type',
title='RNA Velocity Differentiation Stream'
)3. Interpreting Phase Portraits & Latent Time
A gene's phase portrait plots unspliced (\(u\)) versus spliced (\(s\)) counts across all cells:
- Points above the steady-state equilibrium line (\(u = rac{\gamma}{eta} s\)) represent cells actively upregulating the gene (positive velocity).
- Points below the line indicate cells where transcription has ceased and mature mRNA is degrading (negative velocity).
# Plot phase portrait for lineage-defining transcription factors
scv.pl.velocity(adata, ['GATA1', 'SPI1', 'RUNX1'], ncols=3)
# Compute reconstructed latent time across developmental trajectories
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', cmap='gnuplot')