Scvelo
K-DenseCommunityCreated by Kuan-lin HuangDistributed by K-Dense
- Skill
- Workflow
- Science and research
- BSD-3-Clause
- Low risk
- Has scripts
AutoGineer summary
Written by AutoGineer from the skill’s own files, not by the author
This skill provides a structured workflow for performing RNA velocity analysis using scVelo. It enables researchers to estimate expression derivatives from spliced and unspliced single-cell RNA counts, fitting deterministic or dynamical models to explore kinetic trajectories. The workflow includes data validation, preprocessing, model fitting, and visualization of phase portraits and latent time. It is designed to be used alongside Scanpy for analyzing directional cell state transitions, while emphasizing that velocity estimates are model-dependent and do not constitute proof of cell fate or causal drivers.
What it does
- Fits deterministic and dynamical kinetic models to spliced/unspliced RNA counts.
- Generates gene-level phase portraits and velocity graphs for trajectory analysis.
- Estimates relative latent time and ranks velocity-associated genes.
- Provides a validated environment configuration to handle specific library compatibility.
- Includes diagnostic tools to assess model fit and velocity coherence.
- Projects high-dimensional velocity estimates onto UMAP embeddings.
When to use it
- When you need to generate directional trajectory hypotheses from single-cell RNA-seq data.
- When analyzing kinetic models of gene expression using spliced and unspliced counts.
- When you require gene-level phase portraits to inspect induction and repression dynamics.
- When you need to estimate relative latent time within a cell population.
Advantages
- Provides a standardized, tested workflow for complex scVelo analysis.
- Includes explicit validation steps to ensure input data quality and model compatibility.
- Supports both deterministic and dynamical modeling approaches.
- Offers clear diagnostic outputs to help identify potential model misspecification.
Usage examples
- “Run a dynamical RNA velocity analysis on my H5AD file using the clusters annotation.”
- The script will validate the input, perform dynamical model fitting, compute latent time, and save the results including velocity-annotated H5AD files and plots.
- “How do I identify the top genes driving the velocity in my dataset?”
- The skill provides a method to rank velocity-associated genes based on group-specific correlations, which can then be inspected via phase portraits.
Install
Download the skill as a ZIP:
https://autogineer.com/agent-skills/k-dense-ai/scvelo/downloadUnzip it into the cross-client skills folder of your project:
mkdir -p .agents/skills && unzip scvelo.zip -d .agents/skills
Documentation by Kuan-lin Huang
scVelo RNA Velocity
When to use
Use for kinetic analysis of aligned spliced/unspliced RNA counts, directional trajectory hypotheses, gene-level phase portraits, and dynamical latent time. Use the model reference for assumptions, transition probabilities, and an optional CellRank handoff.
RNA velocity estimates an expression derivative under a model. Smooth arrows are not observed cell movement, lineage tracing, causal drivers, or proof of future fate. The dynamical estimator is not automatically more accurate than a steady-state estimator on every dataset.
Tested environment and compatibility
The maintained workflow targets scVelo 0.3.4, Scanpy 1.12.4, AnnData 0.13.4, NumPy 2.5.3, pandas 2.3.3, SciPy 1.18.1, Matplotlib 3.11.2, and loompy 3.0.8 on Python 3.13. Use an isolated environment:
uv venv --python 3.13 .venv-velocity
uv pip install --python .venv-velocity/bin/python \
scvelo==0.3.4 scanpy==1.12.4 anndata==0.13.4 numpy==2.5.3 \
pandas==2.3.3 scipy==1.18.1 matplotlib==3.11.2 loompy==3.0.8The packages above were resolved from cache and tested in an isolated uv
environment; shell paths above use POSIX syntax. Installation requires network
access unless packages are cached. No API key is required.
Deterministic and dynamical models: native synthetic fitting, graphs, relative times, plots and H5AD round trips are tested on this stack.
Stochastic model: scVelo 0.3.4's default GLS scalar assignment fails with NumPy 2. Scanpy 1.12 requires NumPy 2, so installing
numpy<2beside current Scanpy is not a solution. A separately validated legacy Scanpy/NumPy environment is required; that legacy stack was not executed in this review. The helper rejects this combination before modifying data, without changing the model.pandas 3: scVelo 0.3.4 dynamical fitting and some plots use operations no longer supported by pandas 3. Keep the explicit pandas pin.
The upstream stable tutorial still contains removed
scv.readand older preprocessing calls. Check the 0.3.4 source when a tutorial disagrees with the installed API.
Input contract
Obtain spliced and unspliced counts from a velocity-aware quantifier. Keep quantifier/version, reference annotation, counting mode, sample IDs and cell barcode mapping. Upstream quantification is outside the bundled script.
Supply cells by genes matrices in
adata.layers['spliced']andadata.layers['unspliced'], with identical cell/gene ordering and unique IDs. Counts must be finite and nonnegative; fractional count estimates are allowed. Never use logged, scaled, residualized, or batch-corrected values as counts. Numerical inspection alone cannot establish that data are raw.Inspect per-library/cluster coverage, doublets, ambient RNA, zero-count cells, and annotation compatibility before fitting. Preserve the original full-gene count file. Retained-gene backups do not preserve filtered-out genes.
Resolve barcode prefixes and gene identifiers explicitly when combining files. Do not silently strip library IDs, intersect away most cells, or make duplicated biological IDs unique without understanding why they repeat.
Read H5AD through AnnData. scVelo 0.3.4 has no scv.read or scv.DataFrame:
import anndata as ad
adata = ad.read_h5ad("velocity_counts.h5ad")For legacy loom input use ad.io.read_loom('counts.loom', X_name='spliced', sparse=True). AnnData 0.13 deprecates loom; convert to H5AD for further work.
Its layers[None] aliases X: do not delete it or iterate all layers as if every
key were a string. The helper accesses only named velocity layers. Native loom
reading and exact, reordered metadata alignment are covered by the tests;
AnnData 0.13 write_loom is not used.
Run the maintained workflow
The bundled script accepts local files and never downloads a demonstration dataset on startup:
MPLBACKEND=Agg .venv-velocity/bin/python scripts/rna_velocity_workflow.py \
velocity_counts.h5ad --mode dynamical --groupby clusters \
--n-top-genes 2000 --n-neighbors 30 --n-jobs 1 --output-dir velocity_results
# Import a loom and align its raw counts to an existing annotation file:
MPLBACKEND=Agg .venv-velocity/bin/python scripts/rna_velocity_workflow.py \
counts.loom --processed-h5ad annotated.h5ad --groupby clusters \
--mode deterministic --output-dir velocity_checkRun these from the skill directory. Input filenames and biological annotations
are illustrative; the same CLI and functions are tested on small generated
kinetic fixtures. Omit --groupby when annotations are unavailable. Use
--no-plots for analysis without computing UMAP.
The helper performs these steps:
Validate layers, IDs, model compatibility and grouping; refuse previously generated moments/velocity so preprocessing cannot silently run twice.
Back up retained raw layers as
spliced_countsandunspliced_countsand rebuildXfrom raw spliced counts. Filter genes and normalizeXand the two count layers on a linear scale.Apply
sc.pp.log1ponly to X, select highly variable genes with Scanpy, and rebuild PCA and neighbors after subsetting. Reusing stale PCA/neighbors from an unrelated feature set can silently change the velocity model.Compute
MsandMuwithscv.pp.moments(adata, n_neighbors=None)from the explicit Scanpy graph. Moments are dense: budget memory for multiple cells-by-genes arrays, not only the sparse input.For dynamical mode, call
recover_dynamics(var_names='all')on the selected genes beforevelocity(mode='dynamical'); otherwise fit the selected model. Require usable genes and a nonempty velocity graph.Compute velocity coherence and velocity pseudotime; add latent time only for dynamical fits. Rank velocity-associated genes only with at least two groups and at least two cells per group. Record/warn about nonfinite or constant time/coherence outputs in
uns['velocity_workflow']['diagnostics']; a constant pseudotime does not support a trajectory ordering. Ranking is exploratory.For plots, recompute UMAP from the rebuilt graph, project velocities, save PNGs directly to the requested directory, then save H5AD with package versions and selected parameters. Existing labels remain annotations, not validated cell identities.
The Python function mutates its argument in place; pass adata.copy() to keep
the original object. The CLI writes to output_dir/adata_velocity.h5ad.
Version 2 changes the former helper's behavior deliberately: raw layer geometry
is rebuilt, missing requested labels are errors, and plots/parallelism can be
controlled explicitly. Runtime depends on cells, genes and fit difficulty;
there is no universal 10–30 minute expectation.
Inspect the evidence before interpretation
After a dynamical run (dataset-specific gene selection is illustrative):
import pandas as pd
import scvelo as scv
# Candidate kinetic genes, not experimentally established drivers.
candidates = adata.var['fit_likelihood'].dropna().nlargest(6).index.tolist()
scv.pl.velocity(adata, var_names=candidates, basis='umap', show=False)
# This ranks group-associated velocities, not causal influence or condition DE.
scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=0.3)
ranked = pd.DataFrame(adata.uns['rank_velocity_genes']['names'])
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, color=['velocity_length', 'velocity_confidence'], show=False)Inspect spliced-versus-unspliced phase portraits, coverage across induction and repression, fitted parameters, failed/NaN fits, and branch-specific kinetics. Check sensitivity to gene set, neighbors, subsampling and model assumptions. Use time-course labels, perturbations, lineage tracing or labeling experiments as independent directional evidence where available. A plausible UMAP alone cannot validate a fit, and tuning until arrows match a desired story is not a validation strategy.
There is no general minimum of 2,000 cells, universal unspliced fraction, or rule that root cells must have the highest unspliced/spliced ratio. Coverage of relevant kinetic states and measurement quality matter. Negative velocity can represent repression or model misspecification; it is not by itself evidence that layers were swapped.
Output interpretation
| Field | Interpretation |
|---|---|
layers['velocity'] | Model-estimated derivative in processed gene-expression space |
var['velocity_genes'] | Genes selected for the velocity graph; distinct from all HVGs |
layers['Ms'], layers['Mu'] | Neighbor-averaged linear-scale spliced/unspliced expression |
layers['fit_t'] | Gene-specific fitted time coordinates, dynamical model only |
var['fit_alpha/beta/gamma'] | Fitted rates on the model's inferred scale; not calibrated physical rates |
var['fit_likelihood'] | Relative model-fit diagnostic; not a posterior probability of biological truth |
uns['velocity_graph'] | Sparse positive cosine correlations for candidate transitions; not row-stochastic |
obsm['velocity_umap'] | Projected vectors, only after embedding computation/plotting |
obs['velocity_pseudotime'] | Graph-based relative ordering |
obs['latent_time'] | Coupled dynamical ordering; normally scaled 0–1, not elapsed hours |
obs['velocity_length'] | Processed-space vector magnitude; not physical cell speed |
obs['velocity_confidence'] | Neighbor velocity coherence, not calibrated uncertainty |
For fate probabilities use an explicitly normalized transition kernel and a validated terminal-state definition; see the optional CellRank example in the reference. PAGA requires optional igraph and compatible Scanpy internals; it is not part of the tested core workflow or a substitute for fate inference.
Troubleshooting
| Problem | Action |
|---|---|
| Missing or mismatched layers | Revisit quantification and explicit ID alignment; X cannot substitute for unspliced counts |
| Very few velocity genes | Inspect depth, state coverage and phase portraits before altering thresholds |
| Smooth but implausible arrows | Check model assumptions, batch geometry and individual genes; compare independent evidence |
| Nonfinite fit or empty graph | Stop interpretation and investigate degenerate features/coverage |
| Excessive memory or fitting time | Use a justified gene set and n_jobs=1; account for dense moments and fit arrays |
| Stochastic failure on NumPy 2 | Use a separately validated compatible legacy stack or explicitly reconsider the model |
pandas unique error in fit/plot | Use the tested pandas 2.3.3 pin |
Sources and verification scope
Reviewed 2026-10-01 against the released source, API, kinetic-model caveats, and AnnData loom reader. The main paper is Bergen et al., 2020. Synthetic execution tests verify software contracts and file/figure production; they do not establish biological accuracy, identifiability or dataset-specific parameter recovery. Optional CellRank/PAGA analyses remain source-reviewed, illustrative extensions. No authenticated remote service is used.
Files
- ATTRIBUTION.mdAdded by AutoGineer899 B
- SKILL.md12.2 KiB
- references/velocity_models.md9.2 KiB
- scripts/rna_velocity_workflow.pyScript12.9 KiB
references/velocity_models.md
Velocity models and interpretation
Reviewed 2026-10-01 for scVelo 0.3.4. API examples below assume the
preprocessed adata from the main workflow, including Ms, Mu, a neighbors
graph and, for dynamical operations, fitted dynamics. The core runtime stack
and its stochastic/pandas limitations are specified in SKILL.md.
Kinetic framework
For unspliced abundance u and spliced abundance s:
du/dt = alpha(t) - beta*u
ds/dt = beta*u - gamma*s
velocity = ds/dtA positive velocity means increasing spliced abundance under the fitted model, not guaranteed induction of an entire cell program. Rates, relative scaling of the two modalities and inferred time can trade off. Snapshot data alone do not identify an absolute biological clock.
| Model | Assumption and appropriate diagnostic |
|---|---|
deterministic (steady_state alias) | Fits a steady-state spliced/unspliced relationship from expression extremes; requires informative steady-state coverage |
stochastic | Adds second-order moments to the steady-state treatment; depends on neighborhood estimates and the same biological assumptions |
dynamical | Fits induction/repression kinetics and cell-specific times, relaxing observed steady-state requirements but retaining the specified kinetic model |
The dynamical model usually costs more. None of these names guarantees accuracy or publication readiness. Incomplete cycles, changing rates across branches, transcriptional bursts, mature populations without true transitions and measurement biases can make seemingly coherent estimates wrong.
# Tested on the current core stack:
scv.tl.velocity(adata, mode='deterministic')
# Dynamical fitting must precede dynamical velocity:
scv.tl.recover_dynamics(adata, var_names='all', n_jobs=1, show_progress_bar=False)
scv.tl.velocity(adata, mode='dynamical')These are alternative analyses on separately prepared objects. Default
recover_dynamics(var_names='velocity_genes') can select genes using an initial
steady-state estimate; explicit var_names='all' fits all genes in the already
filtered analysis object. It does not restore excluded genes.
The default stochastic GLS path is not supported by the tested NumPy 2 stack. Do not enable offsets or substitute another estimator solely to bypass an API failure: those choices change model assumptions. A legacy NumPy<2 environment also needs an older compatible Scanpy/SciPy/Python combination, which was not executed in this review.
Graphs and projections
scv.tl.velocity_graph(adata, n_jobs=1, show_progress_bar=False)
T = scv.utils.get_transition_matrix(adata)uns['velocity_graph'] stores positive cosine correlations between expression
changes and velocity; velocity_graph_neg stores the negative part. They are
not normalized probabilities. The transition utility applies a kernel and
normalization; check finite/nonnegative entries, each row sum, isolated cells,
and self-transition choices before using T for a Markov calculation. Changing
the kernel, connectivity blend or terminal states changes inferred fates.
For velocity_graph, approx=True performs cosine calculations in a reduced
PCA representation (first 30 PCs); it does not select an approximate nearest
neighbor backend. sqrt_transform=None is model-dependent in the released
implementation. Set it explicitly only with a scientific reason.
UMAP streams are interpolated projections of high-dimensional estimates, with additional grid smoothing. Arrows can bridge populations that are not real transitions. Check gene-space phase portraits and sensitivity to embedding, neighbors and selected genes before interpreting projected topology.
Relative time and fit diagnostics
layers['fit_t'] contains gene-specific fitted coordinates; it does not impose
universal landmarks such as peak induction at 0.5. tl.latent_time filters genes
by likelihood, aligns/root-adjusts gene times, builds a shared time and smooths
it using neighborhood information. It can incorporate velocity-pseudotime
information. It is not simply a likelihood-weighted mean.
scv.tl.latent_time(adata) # requires fitted dynamics and a velocity graph
scv.pl.scatter(adata, color='latent_time', basis='umap', show=False)The default output is scaled to 0–1. t_max rescales that ordering; supplying a
number does not measure hours. root_key and end_key encode root/end priors;
record external evidence for those priors and assess direction sensitivity.
There is no t_max_rank argument in the current latent-time API.
fit_likelihoodis a fit score, not causal evidence or a calibrated probability that a gene is biologically correct. Threshold 0.1 is an implementation default in some downstream methods, not a universal scientific cutoff.fit_r2reflects the steady-state regression used during fitting/selection; it is not a complete validation of the nonlinear kinetic model.velocity_lengthis processed-space vector magnitude, dependent on feature scaling and gene selection. It is not comparable physical cell speed.scVelo 0.3.4 clips negative neighbor correlations to zero for
velocity_confidence; ordinary finite scores lie around 0–1. Zero-vector degeneracy can produce NaNs. Report nonfinite cells, rather than describing a coherence cutoff as a probability of correctness.rank_velocity_genescompares group-associated velocities. It does not show perturbational causality or replace biological-replicate-aware condition DE.
Diagnostic plotting uses existing functions:
scv.pl.proportions(adata, groupby='clusters', show=False)
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, color='velocity_confidence', show=False)There is no scv.pl.velocity_confidence function in 0.3.4. Use pandas
DataFrame for the structured gene-ranking table, not removed
scv.DataFrame aliases.
Preprocessing with Scanpy
In 0.3.4, filter_and_normalize filters genes and normalizes per cell; its stale
docstring still mentions log transformation. It no longer performs logging/HVG
selection. Apply sc.pp.log1p to X and sc.pp.highly_variable_genes explicitly.
Spliced/unspliced layers remain linear. A preprocessed Scanpy X is not evidence
that those layers have been normalized correctly.
Rebuild PCA and neighbors after selecting the velocity feature set. Reusing an integration-derived graph is a separate analysis choice: batch correction changes local state geometry and can invent links or remove temporal signal. This helper deliberately uses uncorrected log-spliced PCA. Any separately chosen integration must follow an independently validated integration workflow.
pp.moments writes dense first moments (Ms, Mu). The stochastic estimator
computes second-order moments when needed; do not describe the two stored
first-moment arrays as means and variances.
Optional CellRank 2.3.3 handoff (illustrative; not executed here)
Install and validate CellRank in a compatible isolated environment before use. This extension is not a core dependency of the scVelo helper. The current GPCCA API requires terminal-state assignment before fate probabilities.
import cellrank as cr
from cellrank.kernels import VelocityKernel, ConnectivityKernel
vk = VelocityKernel(adata).compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
combined = 0.8 * vk + 0.2 * ck # illustrative weight; perform sensitivity analysis
estimator = cr.estimators.GPCCA(combined)
estimator.compute_macrostates(n_states=4, cluster_key='clusters')
estimator.plot_macrostates(which='all')
# Propose terminal states, then inspect and validate them biologically.
estimator.predict_terminal_states()
estimator.plot_macrostates(which='terminal')
# Only proceed after checking terminal states and connectivity:
estimator.compute_fate_probabilities()
estimator.plot_fate_probabilities()Four macrostates and the 0.8 velocity weight are examples, not defaults or
recommendations. Explicitly set validated states with set_terminal_states
when prior evidence supports them. Probabilities remain conditional on the
kernel, graph, observed cells and selected terminal states.
PAGA is a separate coarse cluster-connectivity summary. scv.tl.paga and
scv.pl.paga remain documented, but need optional igraph and rely on Scanpy's
private PAGA implementation. They were source-reviewed, not executed in this
core test environment. A PAGA edge does not establish a cell fate probability.
Primary references
Requirements and compatibility
Works with: Claude Code, Claude apps, OpenAI Codex, GitHub Copilot, Cursor, Gemini CLI, Kiro, OpenCode, Goose, Amp, Roo Code, JetBrains Junie, Factory Droid, Devin (ex-Windsurf), Mistral Vibe, OpenHands, Letta
“Works with” means the tool reads the open Agent Skills format; a skill written for one tool may use features others lack.
Needs: Needs an API key, Needs Python
The author says: Requires Python 3.13 for the tested stack with scvelo 0.3.4, scanpy 1.12.4, anndata 0.13.4, numpy 2.5.3 and pandas 2.3.3. Loom import requires loompy. Local H5AD analysis needs no network or credentials. The default stochastic solver is incompatible with NumPy 2.
Security
What our scanner noticed· Low risk
- Instructions that use the network
- Scripts that your agent may run
This skill includes scripts that your agent may run. Read them before you use it.
A static scan is not a guarantee. Skills are third-party content: review a skill before you let an agent use it.
Attribution and licence
- Created by
- Kuan-lin Huang
- Publisher
- K-Dense (Community)
- Distributed by
- K-Dense
- Source
- K-Dense-AI/scientific-agent-skills · skills/scvelo · commit a25cf9e
- Licence
- BSD 3-Clause License (redistribution allowed)
- ATTRIBUTION.md (written by AutoGineer)
- Synced
AutoGineer is not affiliated with or endorsed by Kuan-lin Huang or K-Dense.
Download ZIP
The ZIP holds the skill’s files, its licence files and an ATTRIBUTION.md. Download scvelo.zip