---

SNAP: Spatial Niche Architecture Profiler

SNAP is an R toolkit for analysis after spatial Regions have been annotated. It connects spot-level local texture, Region descriptors, Region-level graphs, expression programs and graph-aware comparisons while retaining ordinary Seurat, matrix, data-frame and igraph inputs.

The examples below show the compact pipelines intended for routine use. The lower-level functions remain available when an analysis requires custom signals, alternative graph weights or additional interpretation.

Installation

SNAP has been functionally tested under R 4.4.1. During development, install the package directly from GitHub.

install.packages("remotes")
remotes::install_github(
  "coy6/SNAP",
  dependencies = TRUE,
  upgrade = "never"
)

Alternatively, with pak:

install.packages("pak")
pak::pak("coy6/SNAP")

Region descriptors

coords must contain x, y, barcode and annot, in the same spot order as label and signal. A boundary spot is a spot participating in at least one neighbor pair whose Region annotations differ.

library(SNAP)

desc_res <- region_descriptor_pipeline(
  seurat_obj = seurat_obj,
  coords = coords,
  label = seurat_obj[[]][coords$barcode, "cell_state", drop = TRUE],
  signal = as.numeric(expr["ExampleGene", coords$barcode]),
  expr_matrix = expr,
  assay = "SCT",
  k = 6,
  verbose = TRUE
)

region_descriptors <- desc_res$descriptors

Seurat annotation to NMF meta-programs

Signal selection is particularly important in regional analysis which reduces computational complexity. The SNAP package adapts and accelerates the meta-program algorithm from previous work (Gavish A, et al., Nature (2023), 10.1038/s41586-023-06130-4), computing NMF programs on a regional basis and integrating them into meta-programs.

The meta-program pipeline extracts a non-negative Seurat assay layer, runs repeated NMF within Regions, filters unstable and redundant programs, and merges cross-Region-supported programs into meta-programs. For large data sets, pass a pre-selected gene set.

meta_res <- region_metaprogram_pipeline(
  seurat_obj = seurat_obj,
  region_col = "region",
  assay = "SCT",
  layer = "data",
  gene_subset = Seurat::VariableFeatures(seurat_obj),
  ranks = 3:8,
  nrun = 10,
  top_n = 50,
  intra_min = 35,
  intra_max = 10,
  inter_filter = "shared",
  inter_min = 10,
  min_intersect_initial = 10,
  min_intersect_cluster = 10,
  min_meta_group_size = 5,
  seed = 123,
  verbose = TRUE
)

# Named top-gene signatures
meta_res$meta_programs

# Constituent NMF programs and their source Regions
meta_res$program_membership

# Programs that did not enter a meta-program
meta_res$unclustered_programs

The default inter_filter = "shared" is appropriate for meta-program discovery: a retained program must overlap a program from another Region by at least inter_min genes. Use "specific" when the aim is Region-restricted programs, or "none" when cross-Region selection is performed externally. NMF input must be non-negative; do not use scale.data. If a Seurat v5 assay remains split across multiple sample-specific layers, join the relevant layers before running the pipeline so that the selected layer contains every analyzed spot.

Compact FGW alignment

Graphs must be undirected and simple. Feature matrices have Regions on rows, the same descriptors on columns, and row names matching graph vertex names. When descriptors have different units, standardize them jointly across the two samples before alignment.

g1 <- igraph::simplify(igraph::as_undirected(region_graphs[[1]]))
g2 <- igraph::simplify(igraph::as_undirected(region_graphs[[2]]))

feature_cols <- c(
  "main_fraction", "boundary_frac", "anisotropy_mean",
  "boundary_label_entropy_mean", "boundary_signal_variance_mean"
)

x1 <- descriptors[[1]][igraph::V(g1)$name, feature_cols, drop = FALSE]
x2 <- descriptors[[2]][igraph::V(g2)$name, feature_cols, drop = FALSE]
stopifnot(!anyNA(x1), !anyNA(x2))
x_all <- scale(rbind(x1, x2))
x1 <- x_all[seq_len(nrow(x1)), , drop = FALSE]
x2 <- x_all[nrow(x1) + seq_len(nrow(x2)), , drop = FALSE]

fgw_res <- fgw_pipeline(
  g1 = g1,
  g2 = g2,
  features1 = x1,
  features2 = x2,
  alpha = 0.5,
  p = region_mass_1,
  q = region_mass_2
)

fgw_res$loss
fgw_res$best_match
fgw_res$transport
fgw_res$edge_conservation$global
fgw_res$neighborhood_concordance$bidirectional

best_probability is the conditional transport probability given a source Region, not a hard one-to-one assignment. Edge conservation and neighborhood concordance are returned in both directions because the source and target graphs can have different sizes and degrees. Compare loss only across runs using the same feature set, graph construction and parameterization.

Compact graph OT across ordered slices

Each expression matrix must have features on rows and graph Region names on columns. Graph OT normalizes every analyzed feature within each slice; it tests redistribution among Regions, not total-abundance change.

ot_hits <- graph_ot_pipeline(
  g_list = region_graphs,
  expr_list = region_expression,
  genes = genes_of_interest,
  mass_threshold = 0.20,
  slice_names = c("E9.5", "E10.5", "E11.5", "E12.5")
)

head(ot_hits)

# Features with supported transport events in each transition
split(ot_hits$feature, ot_hits$step) |>
  lapply(unique)

# Inspect events for one feature
subset(ot_hits, feature == "ExampleGene")

The returned data frame contains only plan rows with mass > mass_threshold, together with total_cost and bigM_ratio. A large bigM_ratio indicates that substantial mass used an artificial high-cost connection caused by graph disconnection and should be inspected before biological interpretation.

For compatibility with plotting functions that require the original nested result, keep the lower-level representation and filter it without truncating individual plans:

ot_raw <- compute_ot_across_slices(region_graphs, region_expression,
                                   genes = genes_of_interest)
ot_nested <- filter_ot_list(ot_raw$step_results, mass_threshold = 0.20)