Vector Field Analysis and Visualization

Introduction

One of the most powerful features of CellODE is its ability to infer and visualize cellular velocity fields. The vector field represents the instantaneous direction and magnitude of cellular state changes, providing insights into:

Differentiation trajectories
Cell fate decisions
Developmental dynamics

This vignette demonstrates how to use CellODE’s built-in functions for vector field analysis and visualization.

Quick Start: Using CellODE Functions

library(CellODE)
library(ggplot2)
library(Seurat)

Create Example Data

First, let’s create a Seurat object with simulated trajectory data:

set.seed(42)
n_cells <- 500
n_genes <- 100

# Simulate trajectory
t_true <- seq(0, 1, length.out = n_cells)

# Create expression data (genes change along trajectory)
counts <- matrix(0, nrow = n_genes, ncol = n_cells)
for (i in 1:n_genes) {
    phase <- runif(1, 0, 2 * pi)
    amplitude <- runif(1, 3, 10)
    counts[i, ] <- rpois(n_cells, lambda = pmax(1, amplitude * (1 + sin(t_true * 2 * pi + phase))))
}
rownames(counts) <- paste0("Gene", seq_len(n_genes))
colnames(counts) <- paste0("Cell", seq_len(n_cells))

# Create Seurat object
srt <- CreateSeuratObject(counts = counts)

# Add UMAP embedding (curved trajectory)
umap1 <- t_true * 10 + rnorm(n_cells, sd = 0.5)
umap2 <- 4 * sin(t_true * 2 * pi) + rnorm(n_cells, sd = 0.3)

umap_embeddings <- matrix(c(umap1, umap2), ncol = 2)
rownames(umap_embeddings) <- colnames(srt)
colnames(umap_embeddings) <- c("UMAP_1", "UMAP_2")

srt[["umap"]] <- CreateDimReducObject(
    embeddings = umap_embeddings,
    key = "UMAP_",
    assay = "RNA"
)

# Add pseudotime
srt$pseudotime <- t_true + rnorm(n_cells, sd = 0.02)
srt$pseudotime <- pmax(0, pmin(1, srt$pseudotime))

Simulate CellODE Model Output

In a real workflow, you would train a CellODE model. Here we simulate the output:

# Simulate latent space (normally from model output)
n_latent <- 5
zs <- matrix(rnorm(n_cells * n_latent), ncol = n_latent)

# Make latent space correlate with trajectory
zs[, 1] <- t_true * 2 + rnorm(n_cells, sd = 0.2)
zs[, 2] <- sin(t_true * pi) + rnorm(n_cells, sd = 0.1)
rownames(zs) <- colnames(srt)

# Simulate vector field in latent space
vf <- matrix(0, nrow = n_cells, ncol = n_latent)
vf[, 1] <- 1 + rnorm(n_cells, sd = 0.1)  # Constant velocity in z1
vf[, 2] <- cos(t_true * pi) * pi / 2 + rnorm(n_cells, sd = 0.1)  # Derivative of sin
rownames(vf) <- colnames(srt)

# Store in Seurat object (as CellODE does)
srt@misc$X_zs <- zs
srt@misc$X_VF <- vf

Vector Field Computation with CellODE

Step 1: Compute Cosine Similarity

Use cosine_similarity() to compute transition probabilities in latent space:

# Compute cosine similarity matrix
T_mat <- cosine_similarity(
    zs = zs,
    vf = vf,
    n_neigh = 20
)

cat("Transition matrix dimensions:", dim(T_mat), "\n")
#> Transition matrix dimensions: 500 500
cat("Non-zero entries:", sum(T_mat != 0), "\n")
#> Non-zero entries: 52714

Step 2: Compute Embedding Velocity

Use vector_field_embedding() to project velocities to UMAP space:

# Get embedding coordinates
E <- Embeddings(srt, "umap")

# Compute velocity in embedding space
V_embed <- vector_field_embedding(
    T_mat = T_mat,
    E = E,
    scale = 10
)

cat("Embedding velocity dimensions:", dim(V_embed), "\n")
#> Embedding velocity dimensions: 500 2
cat("Mean velocity magnitude:", mean(sqrt(rowSums(V_embed^2))), "\n")
#> Mean velocity magnitude: 0.240518

Step 3: Grid-Based Vector Field

Use vector_field_embedding_grid() for cleaner visualization:

# Compute grid-based vector field
grid_result <- vector_field_embedding_grid(
    E = E,
    V = V_embed,
    smooth = 0.5,
    stream = FALSE,
    density = 1.0
)

cat("Grid points:", nrow(grid_result$E_grid), "\n")
#> Grid points: 830

Visualization with CellODE

Plot Pseudotime with Direction Arrows

Use plot_pseudotime() to visualize pseudotime with trajectory arrows:

# Basic pseudotime plot
p1 <- plot_pseudotime(
    srt, 
    time_key = "pseudotime",
    arrow_show = TRUE,
    arrow_density = 1.5,
    point_palette = "Spectral",
    title = "Pseudotime with Direction Arrows"
)
print(p1)

Plot Vector Field - Grid Mode

Use plot_vector_field() with grid mode for clean visualization:

# Vector field - grid mode
p2 <- plot_vector_field(
    srt,
    plot_type = "grid",
    color_by = "pseudotime",
    arrow_density = 0.8,
    arrow_color = "grey20",
    point_palette = "Spectral",
    title = "Vector Field (Grid Mode)"
)
print(p2)

Plot Vector Field - Stream Mode

Stream mode provides colored arrows based on velocity magnitude:

# Vector field - stream mode with colored arrows
p3 <- plot_vector_field(
    srt,
    plot_type = "stream",
    color_by = "pseudotime",
    arrow_density = 1.2,
    point_palette = "viridis",
    title = "Vector Field (Stream Mode)"
)
print(p3)

Plot Vector Field - Raw Mode

Raw mode shows arrows at each cell position:

# Vector field - raw mode (per-cell arrows)
p4 <- plot_vector_field(
    srt,
    plot_type = "raw",
    color_by = "pseudotime",
    arrow_density = 0.3,  # Subsample for clarity
    arrow_color = "navy",
    point_palette = "magma",
    title = "Vector Field (Raw Mode)"
)
print(p4)

Customization Options

Arrow Styling

# Custom arrow styling
p_custom <- plot_vector_field(
    srt,
    plot_type = "grid",
    color_by = "pseudotime",
    arrow_type = "closed",  # Filled arrowheads
    arrow_color = "#E63946",  # Red arrows
    arrow_density = 0.6,
    arrow_length = 0.2,
    arrow_width = 0.8,
    arrow_angle = 30,
    point_palette = "RdYlBu",
    point_size = 0.8,
    title = "Custom Arrow Styling"
)
print(p_custom)

Clean Publication Figure

# Clean figure for publication
p_pub <- plot_pseudotime(
    srt,
    time_key = "pseudotime",
    arrow_show = TRUE,
    arrow_density = 1.0,
    arrow_type = "open",
    arrow_color = "black",
    point_palette = "viridis",
    point_size = 1.2,
    show_axes = FALSE,
    title = NULL
)
print(p_pub)

Mathematical Background

From Latent Dynamics to Velocity

The Neural ODE in CellODE defines the velocity in latent space:

$\mathbf{v}_\text{latent}(t) = \frac{d\mathbf{z}}{dt} = f_\theta(\mathbf{z}(t), t)$

Cosine Similarity

The cosine_similarity() function computes:

$\text{cos\_sim}_{ij} = \frac{\mathbf{v}_i \cdot (\mathbf{z}_j - \mathbf{z}_i)}{||\mathbf{v}_i|| \cdot ||\mathbf{z}_j - \mathbf{z}_i||}$

Transition Probability

The transition probabilities are computed as:

$T_{ij} = \text{sign}(\text{cos\_sim}_{ij}) \cdot (\exp(|\text{cos\_sim}_{ij}| \cdot s) - 1)$

Embedding Velocity

The vector_field_embedding() function computes:

$\mathbf{v}_i^\text{embed} = \sum_j T_{ij} \cdot \frac{\mathbf{e}_j - \mathbf{e}_i}{||\mathbf{e}_j - \mathbf{e}_i||}$

Comparison with RNA Velocity

Aspect	CellODE	RNA Velocity (scVelo)
Input	Gene expression only	Spliced/unspliced counts
Model	Neural ODE	Steady-state/dynamic model
Time inference	Joint with dynamics	Separate
R package	Native R	Python (reticulate needed)

Best Practices

1. Neighbor Selection

# More neighbors = smoother but less detailed
T_smooth <- cosine_similarity(zs, vf, n_neigh = 30)

# Fewer neighbors = more detail but noisier
T_detail <- cosine_similarity(zs, vf, n_neigh = 10)

2. Visualization Density

# Publication figure (fewer arrows)
plot_vector_field(srt, arrow_density = 0.3)

# Exploration (more arrows)
plot_vector_field(srt, arrow_density = 1.0)

3. Self-Transition

# Include self-transition for smoother results
V_with_self <- vector_field_embedding(T_mat, E, self_transition = TRUE)

Summary

CellODE provides a complete workflow for vector field analysis:

Function	Purpose
`cosine_similarity()`	Compute transition probabilities
`vector_field_embedding()`	Project to embedding space
`vector_field_embedding_grid()`	Grid-based smoothing
`plot_vector_field()`	Visualize vector field
`plot_pseudotime()`	Visualize pseudotime with arrows

Session Info

sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.6.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] C
#> 
#> time zone: Asia/Shanghai
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] SeuratObject_4.1.4 Seurat_4.4.0       ggplot2_4.0.1      CellODE_1.0.0     
#> 
#> loaded via a namespace (and not attached):
#>   [1] deldir_2.0-4           pbapply_1.7-4          gridExtra_2.3         
#>   [4] rlang_1.1.7            magrittr_2.0.4         RcppAnnoy_0.0.23      
#>   [7] otel_0.2.0             spatstat.geom_3.7-0    matrixStats_1.5.0     
#>  [10] ggridges_0.5.7         compiler_4.4.0         reshape2_1.4.5        
#>  [13] png_0.1-8              systemfonts_1.3.1      callr_3.7.6           
#>  [16] vctrs_0.7.1            stringr_1.6.0          pkgconfig_2.0.3       
#>  [19] fastmap_1.2.0          labeling_0.4.3         promises_1.5.0        
#>  [22] rmarkdown_2.30         ps_1.9.1               ragg_1.5.0            
#>  [25] purrr_1.2.1            torch_0.16.3           bit_4.6.0             
#>  [28] xfun_0.56              cachem_1.1.0           jsonlite_2.0.0        
#>  [31] goftest_1.2-3          later_1.4.5            spatstat.utils_3.2-1  
#>  [34] irlba_2.3.5.1          parallel_4.4.0         cluster_2.1.8.1       
#>  [37] R6_2.6.1               ica_1.0-3              spatstat.data_3.1-9   
#>  [40] stringi_1.8.7          bslib_0.9.0            RColorBrewer_1.1-3    
#>  [43] reticulate_1.44.1      spatstat.univar_3.1-6  parallelly_1.46.1     
#>  [46] scattermore_1.2        lmtest_0.9-40          jquerylib_0.1.4       
#>  [49] Rcpp_1.1.1             knitr_1.51             tensor_1.5.1          
#>  [52] future.apply_1.20.1    zoo_1.8-15             sctransform_0.4.3     
#>  [55] httpuv_1.6.16          Matrix_1.7-4           splines_4.4.0         
#>  [58] igraph_2.2.1           tidyselect_1.2.1       abind_1.4-8           
#>  [61] dichromat_2.0-0.1      yaml_2.3.12            spatstat.random_3.4-4 
#>  [64] spatstat.explore_3.7-0 codetools_0.2-20       miniUI_0.1.2          
#>  [67] processx_3.8.6         listenv_0.10.0         plyr_1.8.9            
#>  [70] lattice_0.22-7         tibble_3.3.1           shiny_1.12.1          
#>  [73] withr_3.0.2            S7_0.2.1               ROCR_1.0-12           
#>  [76] evaluate_1.0.5         Rtsne_0.17             future_1.69.0         
#>  [79] desc_1.4.3             survival_3.8-3         polyclip_1.10-7       
#>  [82] fitdistrplus_1.2-6     pillar_1.11.1          KernSmooth_2.23-26    
#>  [85] plotly_4.12.0          generics_0.1.4         sp_2.2-0              
#>  [88] scales_1.4.0           coro_1.1.0             globals_0.18.0        
#>  [91] xtable_1.8-4           glue_1.8.0             lazyeval_0.2.2        
#>  [94] tools_4.4.0            data.table_1.18.0      RANN_2.6.2            
#>  [97] fs_1.6.6               leiden_0.4.3.1         cowplot_1.2.0         
#> [100] grid_4.4.0             tidyr_1.3.2            nlme_3.1-168          
#> [103] patchwork_1.3.2        cli_3.6.5              spatstat.sparse_3.1-0 
#> [106] textshaping_1.0.4      viridisLite_0.4.2      dplyr_1.1.4           
#> [109] uwot_0.2.4             gtable_0.3.6           sass_0.4.10           
#> [112] digest_0.6.39          progressr_0.18.0       ggrepel_0.9.6         
#> [115] htmlwidgets_1.6.4      farver_2.1.2           htmltools_0.5.9       
#> [118] pkgdown_2.1.3          lifecycle_1.0.5        httr_1.4.7            
#> [121] mime_0.13              bit64_4.6.0-1          MASS_7.3-65

Zaoqu Liu

2026-01-26