Velocyto/RNA velocity

December 19, 2025 · 6 min read · Tutorials

TL;DR: RNA velocity predicts cell state transitions by comparing nascent (intron-containing) vs mature (spliced) mRNA. Use Velocyto to count introns from BAM files, then scVelo for analysis and visualization. Essential for trajectory inference and understanding dynamic processes.

What is RNA Velocity?

RNA velocity estimates the rate and direction of gene expression changes in individual cells.

The Core Concept

Cells contain two types of mRNA:

Unspliced (with introns) = newly transcribed, nascent RNA
Spliced (no introns) = mature, translated RNA

Key insight: If a gene has more unspliced than spliced mRNA, it's being actively turned on. If the reverse, it's turning off.

RNA Velocity Concept

Why It Matters

Traditional scRNA-seq	RNA Velocity
Static snapshot	Predicts future state
"Cell A is at position X"	"Cell A is moving toward position Y"
Requires pseudotime algorithms	Direct measurement of dynamics
Ambiguous differentiation paths	Resolves branching

When to Use RNA Velocity

✅ Ideal applications:

Differentiation trajectories (stem cells → mature cells)
Cell cycle phase transitions
Immune activation responses
Developmental timecourse studies
Tumor progression dynamics

⚠️ Limitations:

Requires high-quality unspliced reads (not all protocols capture introns)
Assumes RNA transcription/splicing/degradation equilibrium
Less effective in steady-state systems
Computationally intensive

Workflow Overview

BAM file → Velocyto (count introns) → .loom file → scVelo (analysis) → velocity vectors

Part 1: Velocyto (Counting Introns)

Installation

# Create environment with Python 3.8
conda create -n velocyto python=3.8
conda activate velocyto

# Install dependencies
conda install numpy scipy cython numba matplotlib scikit-learn h5py click

# Install via pip
pip install pysam velocyto

Test installation:

velocyto --version

Prepare Required Files

1. Repeat Masking GTF

Download from UCSC Genome Browser:

Go to UCSC Table Browser
Select:
- Clade: Mammal
- Genome: Human (or your organism)
- Assembly: GRCh38 (match your data)
- Group: Repeats
- Track: RepeatMasker
- Output format: GTF
Click "get output"

Save as: GRCh38_rmsk.gtf

2. Gene Annotation GTF

Use the same GTF from your Cell Ranger/STAR alignment:

# Cell Ranger reference
~/refdata-gex-GRCh38-2024-A/genes/genes.gtf

# Or download from Ensembl
wget ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz
gunzip Homo_sapiens.GRCh38.110.gtf.gz

3. BAM File

From Cell Ranger or STAR alignment:

# Cell Ranger
cellranger_output/outs/possorted_genome_bam.bam

# STAR
STAR_output/Aligned.sortedByCoord.out.bam

4. Cell Barcodes

Extract valid cell barcodes:

# From Cell Ranger
cut -d',' -f2 filtered_feature_bc_matrix/barcodes.tsv.gz | gunzip > barcodes.txt

# From STAR
gunzip -c Solo.out/Gene/filtered/barcodes.tsv.gz > barcodes.txt

Running Velocyto

Single Sample

velocyto run \
  -v \  # Verbose
  -m GRCh38_rmsk.gtf \  # Repeat mask
  -b barcodes.txt \  # Cell barcodes
  -o output_dir \  # Output directory
  -@ 12 \  # Threads
  possorted_genome_bam.bam \  # Input BAM
  genes.gtf  # Gene annotation

Output: output_dir/sample.loom
Time: ~30 min for 5K cells

Batch Processing

#!/bin/bash
# process_velocyto.sh

# Samples to process
samples=(sample1 sample2 sample3)

# Reference files
REPEAT_MASK=~/RefData/GRCh38_rmsk.gtf
GTF=~/RefData/Homo_sapiens.GRCh38.110.gtf

for sample in "${samples[@]}"; do
  echo "Processing $sample..."
  
  # Extract barcodes
  gunzip -c ${sample}/filtered_feature_bc_matrix/barcodes.tsv.gz > ${sample}/barcodes.txt
  
  # Run velocyto
  velocyto run \
    -m $REPEAT_MASK \
    -b ${sample}/barcodes.txt \
    -o velocyto_output/${sample} \
    -@ 8 \
    ${sample}/possorted_genome_bam.bam \
    $GTF
    
  echo "$sample complete"
done

Part 2: scVelo (Analysis & Visualization)

Installation

# Install scVelo and dependencies
pip install -U scvelo loompy

# Optional: R integration for enhanced plots
conda install r-base
pip install rpy2

Basic Analysis Pipeline

import scvelo as scv
import scanpy as sc

# Load data
adata = scv.read("sample.loom", cache=True)

# Basic preprocessing
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# Compute velocity
scv.tl.velocity(adata, mode='stochastic')  # or 'dynamical'
scv.tl.velocity_graph(adata)

# Visualize
scv.pl.velocity_embedding_stream(adata, basis='umap', color='celltype')

Integration with Existing Seurat/Scanpy Analysis

import scvelo as scv
import scanpy as sc

# Load your processed Seurat/Scanpy object
adata_processed = sc.read_h5ad("processed_data.h5ad")

# Load velocyto counts
adata_velo = scv.read("sample.loom", cache=True)

# Merge with processed data
adata = scv.utils.merge(adata_processed, adata_velo)

# Use existing UMAP/clustering
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

# Plot on existing UMAP
scv.pl.velocity_embedding_stream(
    adata, 
    basis='umap', 
    color='celltype',
    figsize=(8, 6),
    arrow_size=3,
    arrow_length=3
)

Velocity Modes

Mode	Description	When to Use
stochastic	Simple steady-state model	Quick exploration, most datasets
dynamical	Full kinetic model	Precise timing, complex dynamics
deterministic	No noise model	Large, high-quality datasets

# Stochastic (fast, good default)
scv.tl.velocity(adata, mode='stochastic')

# Dynamical (slower, more accurate)
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode='dynamical')

Interpretation

# Velocity confidence
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, c='velocity_confidence', cmap='coolwarm')

# Identify driver genes
scv.tl.rank_velocity_genes(adata, groupby='celltype')
scv.pl.velocity_genes(adata, var_names=['TOP1', 'TOP2A', 'MKI67'])

# Velocity pseudotime
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, c='velocity_pseudotime', cmap='gnuplot')

Advanced Workflows

RNA Velocity from STARsolo

Use STARsolo to generate both spliced and unspliced counts directly:

STAR --runThreadN 16 \
  --genomeDir ./STARindex \
  --readFilesIn R2.fastq.gz R1.fastq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloUMIlen 12 \
  --soloFeatures Gene Velocyto \  # Enable velocyto mode
  --outSAMtype BAM SortedByCoordinate

Output: Solo.out/Velocyto/filtered/ contains spliced/unspliced matrices.

Convert to Loom Format

import scvelo as scv
import scanpy as sc

# Load STARsolo velocyto output
spliced = sc.read_mtx("Solo.out/Velocyto/filtered/spliced.mtx.gz").T
unspliced = sc.read_mtx("Solo.out/Velocyto/filtered/unspliced.mtx.gz").T
ambiguous = sc.read_mtx("Solo.out/Velocyto/filtered/ambiguous.mtx.gz").T

# Combine into AnnData
adata = sc.AnnData(spliced)
adata.layers['spliced'] = spliced.X
adata.layers['unspliced'] = unspliced.X
adata.layers['ambiguous'] = ambiguous.X

# Add cell/gene names
barcodes = pd.read_csv("Solo.out/Velocyto/filtered/barcodes.tsv", header=None)
genes = pd.read_csv("Solo.out/Velocyto/filtered/features.tsv", sep="\t", header=None)
adata.obs_names = barcodes[0].values
adata.var_names = genes[0].values

# Save as loom
adata.write_loom("sample.loom")

Visualization Examples

Stream Plot (Recommended)

scv.pl.velocity_embedding_stream(
    adata, 
    basis='umap',
    color='celltype',
    palette='Set2',
    size=80,
    alpha=0.5,
    arrow_size=2,
    linewidth=1.5,
    density=0.8
)

Grid Plot

scv.pl.velocity_embedding_grid(
    adata,
    basis='umap',
    color='celltype',
    arrow_length=3,
    arrow_size=2,
    density=0.5
)

Phase Portrait (Gene-Specific)

# Show specific genes
scv.pl.velocity(adata, ['MKI67', 'TOP2A', 'PCNA'], ncols=3)

Troubleshooting

Low velocity confidence

Symptom: Arrows pointing in random directions
Causes:

Insufficient intron coverage (low sequencing depth)
Steady-state cells (no active transitions)
Poor data quality

Fix:

# Filter low-quality velocity genes
scv.tl.velocity_confidence(adata)
adata = adata[:, adata.var.velocity_confidence > 0.05]

Velocyto crashes with "No valid cells"

Cause: Barcode mismatch between BAM and barcode list
Fix:

# Ensure barcodes match BAM tags exactly
samtools view possorted_genome_bam.bam | grep -o "CB:Z:[ACGT]*" | head
# Compare to barcodes.txt format

Memory errors

Cause: Large dataset
Fix:

# Process in batches
scv.settings.n_jobs = 1  # Reduce parallelization
scv.pp.filter_genes(adata, min_shared_counts=30)  # Stricter filtering

Comparison: Velocyto vs scVelo

Tool	Purpose	Language	Strengths
Velocyto	Count introns from BAM	Python/CLI	Handles any aligner output
scVelo	Analyze velocity	Python	Advanced models, beautiful plots

Workflow: Always use Velocyto → scVelo or STARsolo → scVelo.

Resources

scVelo tutorial: scvelo.readthedocs.io
Velocyto documentation: velocyto.org
scVelo paper: Bergen et al., Nature Biotech 2020
RNA velocity theory: La Manno et al., Nature 2018

See also:

Originally posted on WeChat Official Account: bioinfo_studyroom

RNA Velocity