RegVelo: RNA Velocity Estimation with Variational Inference

RegVelo is a deep generative model that infers splicing kinetics, transcription rates, and latent cellular time by integrating Gene Regulatory Network (GRN) priors.

Unlike previous methods that use decoupled one-dimensional ODEs, RegVelo uses a single, high-dimensional ODE system where transcription is coupled through gene regulation. The model follows these dynamics:

1. Cell Representation: Encodes unspliced (u) and spliced (s) RNA readouts into a low-dimensional latent variable using a neural network.
2. Latent Time: An additional neural network takes samples of the cell representation to parameterize gene-wise latent time.
3. Splicing Dynamics: Modeled via ODEs where transcription ($\alpha_g$) is informed by a GRN weight matrix ($W$) and base transcription ($b$): $$\alpha_g = h \left( \left [ W s(t) +b \right ] _{g} \right)$$ where $h$ is a non-linear activation function.
4. Inference: Uses a parallelizable ODE solver to predict spliced and unspliced abundances. The model assumes predicted abundances are the expected value of a Gaussian likelihood and uses gradient-based optimization to update parameters.
5. Velocity: Cell-gene-specific velocities are defined as splicing velocities based on estimated splicing/degradation rates and predicted abundances.

This approach allows for sampling predicted readouts and velocities from the learned posterior distribution.

Because RegVelo is a generative model that couples cellular dynamics with regulatory networks, it can be used for in silico counterfactual inference to test cellular responses to unseen transcription factor (TF) perturbations.

Perturbation Workflow

1. Deplete the Regulon: For a trained model, simulate a perturbation (e.g., a TF depletion) by removing all its downstream targets from the GRN.
2. Generate Perturbed Field: Generate the perturbed velocity vector field based on the modified GRN.
3. Quantify Effect: Calculate the perturbation effect score by measuring the dissimilarity between the original and perturbed cell velocities. This is typically quantified using cosine dissimilarity.

Connecting to Cell Fate (CellRank Integration)

RegVelo can be combined with CellRank to connect gene regulation to long-term cell fate decisions:

• Local Dynamics: RegVelo provides the vector field.
• Long-term Fate: CellRank uses the vector field to model state transitions via Markov chains and calculate fate probabilities for terminal states.
• Perturbation Analysis: By comparing CellRank's fate probability predictions for the original vs. perturbed vector fields, you can identify enrichment (increased probability) or depletion (decreased probability) effects toward specific terminal states.

After training the RegVelo model, you can use CellRank to compute cell fates and commitment scores.

Workflow:

1. Compute a velocity kernel using cr.kernels.VelocityKernel(adata).compute_transition_matrix() and write it to adata.
2. Use cr.estimators.GPCCA(vk) to compute macrostates.
3. Set terminal states using .set_terminal_states(TERMINAL_STATES).
4. Compute fate probabilities with .compute_fate_probabilities(solver="direct").
5. Visualize commitment scores using rgv.pl.commitment_score.

import cellrank as cr
import regvelo as rgv

# Compute transition matrix
vk = cr.kernels.VelocityKernel(adata).compute_transition_matrix()
vk.write_to_adata()

# Compute macrostates
estimator = cr.estimators.GPCCA(vk)
estimator.compute_macrostates(n_states=9, cluster_key="cell_type")

# Set terminal states and compute fate
TERMINAL_STATES = ["terminal_state_1", "terminal_state_2"]
estimator.set_terminal_states(TERMINAL_STATES)
estimator.compute_fate_probabilities(solver="direct")

# Visualize
estimator.plot_fate_probabilities(same_plot=False, basis="umap")
rgv.pl.commitment_score(
    adata=adata, lineage_key="lineages_fwd", frameon=False,
    s=40, cmap="coolwarm", title="Commitment score",
)

RegVelo allows you to predict the effect of a transcription factor (TF) knockout on cell fate decisions. This is done by removing the TF's outgoing edges from the GRN and recomputing the velocity field.

Workflow:

1. Use rgv.tl.in_silico_block_simulation to generate a perturbed AnnData object and a perturbed model.
2. Recompute the velocity kernel and fate probabilities for the perturbed system using CellRank.
3. Use rgv.mt.cellfate_perturbation to compare the perturbed results against the baseline.
4. Visualize the results with rgv.pl.cellfate_perturbation.

import regvelo as rgv
import cellrank as cr

TF_candidates = ["Gabpa"]
adata_perturb_dict = {}

# 1. Simulate perturbation
for tf in TF_candidates:
    adata_perturb, reg_vae_perturb = rgv.tl.in_silico_block_simulation(
        model="regvelo_model", adata=adata, TF=tf, cutoff=0,
    )
    adata_perturb_dict[tf] = adata_perturb

# 2. Recompute fate for perturbed system
# (Assuming ct_indices is prepared from terminal states)
for tf, adata_target_perturb in adata_perturb_dict.items():
    vkp = cr.kernels.VelocityKernel(adata_target_perturb).compute_transition_matrix()
    estimator_p = cr.estimators.GPCCA(vkp)
    estimator_p.set_terminal_states(ct_indices)
    estimator_p.compute_fate_probabilities(solver="direct")

# 3. Compare and Plot
df = rgv.mt.cellfate_perturbation(perturbed=adata_perturb_dict, baseline=adata, terminal_state=TERMINAL_STATES)
rgv.pl.cellfate_perturbation(
    adata=adata, df=df, fontsize=14, figsize=(8, 4),
    legend_loc="center left", legend_bbox=(1.02, 0.5),
    color_label="cell_type",
)

To use RegVelo, you must provide an anndata.AnnData object and a prior Gene Regulatory Network (GRN).

Requirements:

• AnnData object: Must contain spliced and unspliced counts in adata.layers["spliced"] and adata.layers["unspliced"].
• Prior GRN: A binary pandas.DataFrame or numpy.ndarray where rows are regulators (TFs) and columns are target genes. A value of 1 indicates a regulatory connection, and 0 indicates no connection.

Preprocessing Workflow:

1. Perform standard scRNA-seq preprocessing (QC, normalization, HVG selection, PCA, UMAP) using scvelo and scanpy.
2. Use rgv.pp.preprocess_data to apply velocity gene filtering and min-max scaling.
3. Use rgv.pp.set_prior_grn to integrate the GRN. The resulting skeleton is stored in adata.uns["skeleton"] (regulators as rows, targets as columns).

import scvelo as scv
import scanpy as sc
import pandas as pd
import numpy as np
import regvelo as rgv
import anndata as ad

# Load data
adata = ad.read_h5ad("adata.h5ad")
GRN = pd.read_parquet("prior_GRN.parquet")

# Standard preprocessing
scv.pp.filter_genes(adata, min_shared_counts=20)
scv.pp.normalize_per_cell(adata)
sc.pp.filter_genes_dispersion(adata, n_top_genes=3000)
sc.pp.log1p(adata)
sc.tl.pca(adata, n_comps=30)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=50)
sc.tl.umap(adata)
scv.pp.moments(adata, n_pcs=None, n_neighbors=None)

# RegVelo specific preprocessing
adata = rgv.pp.preprocess_data(adata)
adata = rgv.pp.set_prior_grn(adata, GRN.T)

# Identify TFs present in the skeleton
TF = adata.var_names[adata.uns["skeleton"].sum(1) != 0]
adata.var["TF"] = adata.var_names.isin(TF)

The ModelComparison class is the primary interface for identifying the optimal training strategy and regularization parameters using multi-view biological information. It automates the workflow of training multiple configurations, evaluating them against biological metrics, and visualizing the results.

Workflow

1. Initialize: Create the object with your adata and optional biological metadata (terminal states, transitions).
2. Train: Run train() with a list of strategies and regularization values.
3. Evaluate: Run evaluate() using different biological perspectives (e.g., Pseudotime, Stemness, TSI, CBC).
4. Plot: Run plot_results() to visualize which configuration performs best.

from regvelo import ModelComparison

# Initialize
comp = ModelComparison(adata=adata, terminal_states=TERMINAL_STATES, state_transition=STATE_TRANSITION, n_states=N_STATES)

# Train
comp.train(model_list=['soft', 'hard', 'soft_regularized'], lam2=[0.3, 0.5, 0.8], n_repeat=3)

# Evaluate and Plot
comp.evaluate(side_information='Pseudo_Time')
comp.plot_results(side_information='Pseudo_Time')