Multi-resolution deconvolution of spatial transcriptomics with DestVI and Harreman#
In this tutorial we apply the new version of DestVI (Deconvolution of Spatial Transcriptomics using Variational Inference) to a 10x Visium lymph node dataset, and then feed the deconvolution results into Harreman to infer cell-type-aware metabolic crosstalk.
Why DestVI? Visium spots (55 µm) typically contain several cells. Most deconvolution methods return discrete cell-type proportions, but cells of the same type can exist in different activation states that differ functionally. DestVI models both cell-type proportions and within-cell-type continuous state variation (the “gamma” latent space), enabling downstream analyses that go beyond what cell-type labels alone can reveal.
Tutorial plan:
Load and preprocess the scRNA-seq reference and the Visium spatial dataset.
Train the single-cell Latent Variable Model (scLVM / CondSCVI) on scRNA-seq data.
Train the spatial Latent Variable Model (stLVM / DestVI) to deconvolve each Visium spot.
Visualize cell-type proportions in tissue space.
Explore intra-cell-type gamma variation with spatially weighted PCA.
Perform cell-type-specific differential expression (B cells example).
Integrate DestVI proportions with Harreman to infer cell-type-aware metabolic cell-cell communication (CCC).
#!pip install --quiet --force-reinstall --no-deps git+https://github.com/yoseflab/destvi_utils.git@main
import os
import tempfile
import destvi_utils
import matplotlib.pyplot as plt
import numpy as np
import scanpy as sc
import scvi
import scviva
import seaborn as sns
import torch
from scvi.model import CondSCVI
from scviva.model._destvi import DestVI
import gc
gc.collect()
torch.cuda.empty_cache()
try:
import rapids_singlecell as rsc
print("RAPIDS SingleCell is installed and can be imported")
HAS_RSC = True
except ImportError:
HAS_RSC = False
RAPIDS SingleCell is installed and can be imported
/home/access/miniconda3/lib/python3.13/site-packages/rapids_singlecell/__init__.py:38: UserWarning:
Multiple rapids_singlecell packages are installed: rapids-singlecell, rapids-singlecell-cu13
Please uninstall all versions and reinstall only one:
pip uninstall rapids-singlecell rapids-singlecell-cu13
_detect_duplicate_installation()
scviva.settings.seed = 0
print("Last run with scVIVA-Tools version:", scviva.__version__)
Seed set to 0
Last run with scVIVA-Tools version: 0.1.6
scvi.settings.seed = 0
print("Last run with scvi-tools version:", scvi.__version__)
Seed set to 0
Last run with scvi-tools version: 1.5.0.post1
Note
Modify save_dir below to change where downloaded data files are cached between runs.
sc.set_figure_params(figsize=(6, 6), frameon=False)
sns.set_theme()
torch.set_float32_matmul_precision("high") # use TF32 on Ampere+ GPUs for speed
save_dir = tempfile.TemporaryDirectory()
%config InlineBackend.print_figure_kwargs={"facecolor": "w"}
%config InlineBackend.figure_format="retina"
We work with data from a comparative study of murine lymph nodes: wild-type mice vs. mice injected with mycobacteria (Mycobacterium smegmatis). The dataset consists of:
A 10x Visium spatial transcriptomics section of a lymph node (
ST-LN-compressed.h5ad).A matching 10x Chromium scRNA-seq dataset from the same tissue (
scRNA-LN-compressed.h5ad).
Both files are hosted in the DestVI reproducibility repository and will be downloaded automatically on first run.
out1_path = os.path.join(save_dir.name, "ST-LN-compressed.h5ad")
st_adata = sc.read(
out1_path, backup_url="https://exampledata.scverse.org/scvi-tools/ST-LN-compressed.h5ad"
)
st_adata
AnnData object with n_obs × n_vars = 1092 × 13948
obs: 'in_tissue', 'array_row', 'array_col', 'batch', 'LN', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'n_counts', 'leiden', 'lymph_node'
var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'LN_colors'
obsm: 'X_pca', 'X_umap', 'location', 'modules', 'spatial'
layers: None (.X)
out2_path = os.path.join(save_dir.name, "scRNA-LN-compressed.h5ad")
sc_adata = sc.read(
out2_path, backup_url="https://exampledata.scverse.org/scvi-tools/scRNA-LN-compressed.h5ad"
)
sc_adata
AnnData object with n_obs × n_vars = 14989 × 12854
obs: 'n_genes', 'cell_types', 'batch', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'pred_cell_types', 'doublet_scores', 'doublet_predictions', 'MS', 'louvain_r0.5', 'louvain_r0.7', 'louvain_r1.0', 'leiden_r0.5', 'leiden_r0.7', 'leiden_r1.0', 'DC_A', 'DC_B', 'mono_1', 'mono_2', 'louvain_sub_0.1', 'louvain_sub_0.2', 'louvain_sub_0.3', 'louvain_sub', 'louvain_sub_1', 'louvain_sub_2', 'louvain_sub_3', 'SCANVI_pred_cell_types', 'SCVI_pred_cell_types', 'broad_cell_types'
var: 'gene_ids-0', 'genome-0', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable'
uns: 'batch_colors', 'leiden', 'leiden_r1.0_colors', 'louvain', 'louvain_r0.5_colors', 'louvain_r0.7_colors', 'louvain_r1.0_colors', 'louvain_sub_0.2_colors', 'louvain_sub_0.3_colors', 'louvain_sub_1_colors', 'neighbors', 'pred_cell_types_colors', 'umap'
obsm: 'X_scVI', 'X_umap'
obsp: 'connectivities', 'distances'
layers: None (.X)
Data loading and preprocessing#
Single-cell reference#
We load the scRNA-seq data, which contains immune cells from murine lymph nodes profiled with 10x Chromium. The data come pre-annotated with both broad (broad_cell_types) and fine-grained (cell_types) cluster labels; DestVI will use the broad labels to define cell types and the fine labels to model within-cell-type heterogeneity.
Key rule of thumb: DestVI assumes that at most one cell state per cell type occupies a given spot. If cells of a type span two biologically distinct states (e.g., resting vs. inflamed monocytes) that could co-exist in the same spot, consider splitting them into separate cell types before proceeding.
The UMAP below shows the broad cell-type annotation used as input to DestVI. Twelve major immune populations are represented. DestVI will learn a cell-type-conditioned latent space for each population.
sc.pl.umap(sc_adata, color="broad_cell_types")
# Retain genes detected in at least 10 cells across the dataset
G = 2000
sc.pp.filter_genes(sc_adata, min_counts=10)
sc_adata.layers["counts"] = sc_adata.X.copy() # preserve raw counts before normalization
# Select highly variable genes for model training (seurat_v3 uses raw counts)
sc.pp.highly_variable_genes(
sc_adata, n_top_genes=G, subset=True, layer="counts", flavor="seurat_v3"
)
sc.pp.normalize_total(sc_adata, target_sum=10e4)
sc.pp.log1p(sc_adata)
sc_adata.raw = sc_adata # store normalized data for plotting
Spatial data#
We load the Visium spatial dataset. DestVI requires both datasets to share the same gene set, so we will intersect the variable genes after loading. Raw counts are preserved in a dedicated counts layer (DestVI always operates on raw counts internally).
st_adata.layers["counts"] = st_adata.X.copy() # raw counts for DestVI
st_adata.obsm["spatial"] = st_adata.obsm["location"] # rename to standard key
# Normalize and log-transform for visualization (raw counts are preserved in the layer)
sc.pp.normalize_total(st_adata, target_sum=10e4)
sc.pp.log1p(st_adata)
st_adata.raw = st_adata
# Restrict both datasets to their shared gene set — a prerequisite for DestVI
intersect = np.intersect1d(sc_adata.var_names, st_adata.var_names)
st_adata = st_adata[:, intersect].copy()
sc_adata = sc_adata[:, intersect].copy()
G = len(intersect)
print(f"Shared genes: {G}")
Shared genes: 1888
sc.pl.embedding(st_adata, basis="spatial", color="lymph_node", s=80)
Step 1 — Fit the scLVM (CondSCVI)#
CondSCVI is a cell-type-conditional Variational Autoencoder that learns a low-dimensional latent space for each cell type separately. The latent coordinates (gamma) capture continuous within-cell-type variation (e.g., an interferon-stimulated sub-state of B cells) that DestVI will later map back onto the spatial data.
CondSCVI.setup_anndata registers the raw counts layer, the broad cell-type labels (used as the conditioning variable during training), and the fine-grained labels (used to define vamp_prior_p clusters that shape the prior distribution in the spatial model).
sc_adata.obs.head()
| n_genes | cell_types | batch | n_genes_by_counts | total_counts | total_counts_mt | pct_counts_mt | pred_cell_types | doublet_scores | doublet_predictions | ... | louvain_sub_0.1 | louvain_sub_0.2 | louvain_sub_0.3 | louvain_sub | louvain_sub_1 | louvain_sub_2 | louvain_sub_3 | SCANVI_pred_cell_types | SCVI_pred_cell_types | broad_cell_types | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AAACCCAAGCAGAAAG-1-0 | 1505 | Tregs | 0 | 1505 | 5862.0 | 307.0 | 5.237121 | Tregs | 0.090909 | 0.0 | ... | 6 | 6 | 6 | Tregs | Tregs | 6,1 | 6,1 | Tregs | Tregs | Tregs |
| AAACCCAAGGTCACAG-1-0 | 1883 | Mature B cells | 0 | 1883 | 6673.0 | 346.0 | 5.185074 | Mature B | 0.116751 | 0.0 | ... | 0 | 0 | 0 | B cells,0 | B cells,0 | 0 | 0 | Cycling B/T cells | Mature B | B cells |
| AAACCCACATCGTTCC-1-0 | 2080 | CD8 T cells | 0 | 2080 | 7764.0 | 438.0 | 5.641422 | CD122+ CD8 T | 0.041626 | 0.0 | ... | 1 | 1-3,0 | 1 | CD8 T cells | CD8 T cells | 1-3,0 | 1-3,0 | CD8 T | CD8 T | CD8 T cells |
| AAACCCAGTCGTCTCT-1-0 | 1647 | Cycling B/T cells | 0 | 1647 | 4619.0 | 265.0 | 5.737173 | Mature B | 0.098634 | 0.0 | ... | 5 | 5 | 5 | B cells,7 | B cells,7 | 5 | 5 | Ifit3-high B | Mature B | B cells |
| AAACCCAGTGTAAACA-1-0 | 1857 | Mature B cells | 0 | 1857 | 5965.0 | 339.0 | 5.683152 | Mature B | 0.153374 | 0.0 | ... | 0 | 0 | 0 | B cells,3 | B cells,3 | 0 | 0 | B-macrophage doublets | Mature B | B cells |
5 rows × 31 columns
CondSCVI.setup_anndata(
sc_adata,
layer="counts", # raw counts
labels_key="broad_cell_types", # broad cell-type labels (conditioning variable)
fine_labels_key="cell_types", # fine labels (used for vamp prior clustering)
batch_key="batch",
)
We train CondSCVI for 100 epochs without cell-type reweighting (weight_obs=False), which is appropriate when cell types are roughly balanced. We use a mixture-of-Gaussians prior (prior='mog', num_classes_mog=10) to better capture multi-modal distributions within each cell type.
Training takes ~5 minutes on a GPU.
sc_model = CondSCVI(
sc_adata,
weight_obs=False, # no cell-type reweighting (balanced dataset)
prior="mog", # mixture-of-Gaussians prior for multi-modal cell states
num_classes_mog=10, # 10 mixture components per cell type
)
sc_model.view_anndata_setup()
Anndata setup with scvi-tools version 1.5.0.post1.
Setup via `CondSCVI.setup_anndata` with arguments:
{ │ 'batch_key': 'batch', │ 'labels_key': 'broad_cell_types', │ 'fine_labels_key': 'cell_types', │ 'layer': 'counts', │ 'unlabeled_category': 'unlabeled', │ 'size_factor_key': None }
Summary Statistics ┏━━━━━━━━━━━━━━━━━━┳━━━━━━━┓ ┃ Summary Stat Key ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━┩ │ n_batch │ 4 │ │ n_cells │ 14989 │ │ n_fine_labels │ 16 │ │ n_labels │ 12 │ │ n_vars │ 1888 │ └──────────────────┴───────┘
Data Registry ┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Registry Key ┃ scvi-tools Location ┃ ┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ X │ adata.layers['counts'] │ │ batch │ adata.obs['_scvi_batch'] │ │ fine_labels │ adata.obs['_scvi_fine_labels'] │ │ labels │ adata.obs['_scvi_labels'] │ └──────────────┴────────────────────────────────┘
labels State Registry ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓ ┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩ │ adata.obs['broad_cell_types'] │ B cells │ 0 │ │ │ CD4 T cells │ 1 │ │ │ CD8 T cells │ 2 │ │ │ GD T cells │ 3 │ │ │ Macrophages │ 4 │ │ │ Migratory DCs │ 5 │ │ │ Monocytes │ 6 │ │ │ NK cells │ 7 │ │ │ Tregs │ 8 │ │ │ cDC1s │ 9 │ │ │ cDC2s │ 10 │ │ │ pDCs │ 11 │ └───────────────────────────────┴───────────────┴─────────────────────┘
batch State Registry ┏━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓ ┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃ ┡━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩ │ adata.obs['batch'] │ 0 │ 0 │ │ │ 1 │ 1 │ │ │ 2 │ 2 │ │ │ 3 │ 3 │ └────────────────────┴────────────┴─────────────────────┘
fine_labels State Registry ┏━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓ ┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩ │ adata.obs['cell_types'] │ CD4 T cells │ 0 │ │ │ CD8 T cells │ 1 │ │ │ Cxcl9-high monocytes │ 2 │ │ │ Cycling B/T cells │ 3 │ │ │ GD T cells │ 4 │ │ │ Ifit3-high B cells │ 5 │ │ │ Ly6-high monocytes │ 6 │ │ │ Macrophages │ 7 │ │ │ Mature B cells │ 8 │ │ │ Migratory DCs │ 9 │ │ │ NK cells │ 10 │ │ │ Tregs │ 11 │ │ │ cDC1s │ 12 │ │ │ cDC2s │ 13 │ │ │ pDCs │ 14 │ │ │ unlabeled │ 15 │ └─────────────────────────┴──────────────────────┴─────────────────────┘
sc_model.train(max_epochs=100)
Trainer will use only 1 of 2 GPUs because it is running inside an interactive / notebook environment. You may try to set `Trainer(devices=2)` but please note that multi-GPU inside interactive / notebook environments is considered experimental and unstable. Your mileage may vary.
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
💡 Tip: For seamless cloud logging and experiment tracking, try installing [litlogger](https://pypi.org/project/litlogger/) to enable LitLogger, which logs metrics and artifacts automatically to the Lightning Experiments platform.
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0,1]
/home/access/miniconda3/lib/python3.13/site-packages/lightning/pytorch/utilities/_pytree.py:21: `isinstance(treespec, LeafSpec)` is deprecated, use `isinstance(treespec, TreeSpec) and treespec.is_leaf()` instead.
/home/access/miniconda3/lib/python3.13/site-packages/lightning/pytorch/trainer/connectors/data_connector.py:434: The 'train_dataloader' does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` to `num_workers=63` in the `DataLoader` to improve performance.
`Trainer.fit` stopped: `max_epochs=100` reached.
sc_model.history["elbo_train"].iloc[5:].plot()
plt.show()
The training loss converges quickly. Reducing max_epochs below 200 degrades the quality of the learned gamma representations and is not recommended.
Step 2 — Train the stLVM (DestVI)#
DestVI transfers the decoder network learned by CondSCVI to the spatial domain. For each spot, it infers:
Cell-type proportions — what fraction of the spot’s RNA comes from each cell type.
Cell-type-specific gamma values — where in the CondSCVI latent space the cells of each type in that spot fall.
DestVI.from_rna_model constructs the spatial model directly from the trained CondSCVI model. The smoothed layer (a spatially smoothed version of the raw counts, constructed via a 5-NN graph) is used internally to improve proportion estimation in low-count spots.
# Remove spots with fewer than 10 total counts (empty or near-empty capture areas)
st_adata = st_adata[st_adata.layers["counts"].sum(1) > 10].copy()
st_adata.obs["batch"] = "spatial" # distinguish from scRNA-seq batches in the registry
def spatial_nn_gex_smth(stadata, n_neighs):
# Spatially smooth raw counts by averaging over k nearest neighbors.
if not HAS_RSC:
sc.pp.neighbors(stadata, n_neighs, use_rep="spatial", key_added="Xspatial")
else:
rsc.pp.neighbors(stadata, n_neighs, use_rep="spatial", key_added="Xspatial")
stadata.obsp["Xspatial_connectivities"] = stadata.obsp["Xspatial_connectivities"].ceil()
stadata.obsp["Xspatial_connectivities"].setdiag(1)
return stadata.obsp["Xspatial_connectivities"].dot(stadata.layers["counts"])
# Smoothed counts improve proportion estimation in low-count spots
st_adata.layers["smoothed"] = spatial_nn_gex_smth(st_adata, n_neighs=5)
Key hyperparameters and their effects:
vamp_prior_p(fromsc_model): number of k-means clusters used to shape the variational prior per cell type. More clusters → more gradual gamma transitions across tissue.l1_sparsity: encourages sparser cell-type proportion estimates. Higher values → fewer non-zero cell types per spot.beta_weighting_prior: controls anchoring strength to the scRNA-seq prior. Adjust if library preparation differed substantially between assays.add_celltypes=2: adds catch-all cell-type slots to absorb expression patterns not in the reference.
Training takes ~5 minutes on a GPU.
st_model = DestVI.from_rna_model(
st_adata,
sc_model,
add_celltypes=2,
n_latent_amortization=None,
anndata_setup_kwargs={"smoothed_layer": "smoothed"},
)
st_model.view_anndata_setup()
Anndata setup with scvi-tools version 1.5.0.post1.
Setup via `DestVI.setup_anndata` with arguments:
{'layer': 'counts', 'smoothed_layer': 'smoothed', 'batch_key': 'batch'}
Summary Statistics ┏━━━━━━━━━━━━━━━━━━┳━━━━━━━┓ ┃ Summary Stat Key ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━┩ │ n_batch │ 5 │ │ n_cells │ 1092 │ │ n_vars │ 1888 │ │ n_x_smoothed │ 1888 │ └──────────────────┴───────┘
Data Registry ┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Registry Key ┃ scvi-tools Location ┃ ┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ X │ adata.layers['counts'] │ │ batch │ adata.obs['_scvi_batch'] │ │ ind_x │ adata.obs['_indices'] │ │ x_smoothed │ adata.layers['smoothed'] │ └──────────────┴──────────────────────────┘
batch State Registry ┏━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓ ┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃ ┡━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩ │ adata.obs['batch'] │ 0 │ 0 │ │ │ 1 │ 1 │ │ │ 2 │ 2 │ │ │ 3 │ 3 │ │ │ spatial │ 4 │ └────────────────────┴────────────┴─────────────────────┘
st_model.train(max_epochs=250)
Trainer will use only 1 of 2 GPUs because it is running inside an interactive / notebook environment. You may try to set `Trainer(devices=2)` but please note that multi-GPU inside interactive / notebook environments is considered experimental and unstable. Your mileage may vary.
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
💡 Tip: For seamless cloud logging and experiment tracking, try installing [litlogger](https://pypi.org/project/litlogger/) to enable LitLogger, which logs metrics and artifacts automatically to the Lightning Experiments platform.
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0,1]
/home/access/miniconda3/lib/python3.13/site-packages/lightning/pytorch/utilities/_pytree.py:21: `isinstance(treespec, LeafSpec)` is deprecated, use `isinstance(treespec, TreeSpec) and treespec.is_leaf()` instead.
/home/access/miniconda3/lib/python3.13/site-packages/lightning/pytorch/trainer/connectors/data_connector.py:434: The 'train_dataloader' does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` to `num_workers=63` in the `DataLoader` to improve performance.
/home/access/miniconda3/lib/python3.13/site-packages/lightning/pytorch/loops/fit_loop.py:321: The number of training batches (9) is smaller than the logging interval Trainer(log_every_n_steps=10). Set a lower value for log_every_n_steps if you want to see logs for the training epoch.
`Trainer.fit` stopped: `max_epochs=250` reached.
The training loss (ELBO) should decrease monotonically and plateau before max_epochs. We recommend at least 1,000 epochs for production runs; 250 is sufficient for this tutorial.
st_model.history["elbo_train"].iloc[10:].plot()
plt.show()
Step 3 — Extract and visualize cell-type proportions#
DestVI returns two types of output:
Cell-type proportions (broad resolution): a spot × cell-type matrix of non-negative values summing to 1.
Gamma values (fine resolution): a 5-dimensional latent vector per cell type per spot, capturing within-cell-type state variation.
Cell-type proportions#
get_proportions() extracts the estimated cell-type proportions for every spot. We display a subset of cell types (B cells, CD8 T cells, Monocytes) to illustrate the expected spatial compartmentalization of the lymph node.
# Extract estimated cell-type proportions (spots x cell types, rows sum to 1)
st_adata.obsm["proportions"] = st_model.get_proportions()
st_adata.obsm["proportions"].head(5)
| B cells | CD4 T cells | CD8 T cells | GD T cells | Macrophages | Migratory DCs | Monocytes | NK cells | Tregs | cDC1s | cDC2s | pDCs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AAACCGGGTAGGTACC-1-0 | 0.718202 | 0.013814 | 0.077247 | 0.002520 | 0.061948 | 0.024428 | 0.040254 | 0.001708 | 0.021843 | 0.032114 | 0.001446 | 0.004477 |
| AAACCTCATGAAGTTG-1-0 | 0.685641 | 0.034963 | 0.066043 | 0.005822 | 0.035378 | 0.042503 | 0.032396 | 0.001153 | 0.058732 | 0.034320 | 0.000948 | 0.002100 |
| AAAGACTGGGCGCTTT-1-0 | 0.641712 | 0.034497 | 0.115063 | 0.003290 | 0.048898 | 0.053918 | 0.017097 | 0.003457 | 0.032648 | 0.045555 | 0.001791 | 0.002073 |
| AAAGGGCAGCTTGAAT-1-0 | 0.118764 | 0.102144 | 0.413737 | 0.010971 | 0.037121 | 0.150463 | 0.030685 | 0.008396 | 0.073395 | 0.051870 | 0.000849 | 0.001604 |
| AAAGTCGACCCTCAGT-1-0 | 0.859868 | 0.017001 | 0.011791 | 0.001749 | 0.056792 | 0.008404 | 0.013093 | 0.000551 | 0.005423 | 0.023398 | 0.000614 | 0.001317 |
# Clip at the 99th percentile to reduce the visual impact of extreme outlier spots
ct_list = ["B cells", "CD8 T cells", "Monocytes"]
for ct in ct_list:
data = st_adata.obsm["proportions"][ct].values
st_adata.obs[ct] = np.clip(data, 0, np.quantile(data, 0.99))
sc.pl.embedding(st_adata, basis="spatial", color=ct_list, cmap="Reds", s=80)
... storing 'batch' as categorical
Because stLVM proportion estimates are never exactly zero, follow-up analyses that use cell-type-specific gene expression must first threshold the proportions to distinguish spots where a cell type is genuinely present from spots where its estimated proportion is residual noise.
destvi_utils.automatic_proportion_threshold determines a data-driven cutoff for each cell type. This utility is part of the destvi_utils companion package (installable from GitHub; see the top of this notebook).
ct_thresholds = destvi_utils.automatic_proportion_threshold(
st_adata, ct_list=ct_list, kind_threshold="secondary"
)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:00<00:00, 434.68it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:00<00:00, 44253.05it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:00<00:00, 51250.05it/s]
The results confirm the expected spatial compartmentalization: B cells in the follicle zone, CD8 T cells in T-cell zones, and condition-dependent monocyte distribution (refer to the DestVI paper for full details).
Step 4 — Explore intra-cell-type variation (gamma values)#
Beyond proportions, DestVI provides gamma — a low-dimensional (5D) latent vector per cell type per spot that captures continuous within-cell-type variation. A spot enriched in IFN-stimulated B cells will have different gamma values for the “B cells” cell type than a spot enriched in resting B cells, even if the overall B-cell proportion is similar.
# Store 5D gamma latent vectors for every cell type as obsm entries
# gamma[i, j] = latent coordinate j for cell type i in each spot
for ct, g in st_model.get_gamma().items():
st_adata.obsm[f"{ct}_gamma"] = g
st_adata.obsm["B cells_gamma"].head(5)
| 0 | 1 | 2 | 3 | 4 | |
|---|---|---|---|---|---|
| AAACCGGGTAGGTACC-1-0 | -0.521573 | 0.429240 | 0.059924 | -0.800617 | 0.964216 |
| AAACCTCATGAAGTTG-1-0 | -0.514134 | 0.241981 | 0.335989 | -1.057969 | 1.453098 |
| AAAGACTGGGCGCTTT-1-0 | 0.307449 | -0.105390 | 0.506047 | -0.594057 | 0.739734 |
| AAAGGGCAGCTTGAAT-1-0 | -0.590200 | 0.273401 | 0.387136 | -0.574227 | 0.486441 |
| AAAGTCGACCCTCAGT-1-0 | -0.514111 | 0.129744 | 0.117499 | -0.733029 | 1.413773 |
destvi_utils.explore_gamma_space provides an automated pipeline to interpret the gamma space for each cell type:
Select spots above the proportion threshold for each cell type.
Compute spatially-weighted PCA on the gamma values of selected spots to find the axes of spatial variation.
Color each spot and each reference single cell by their coordinate in the first two sPCs.
Plot (A) colors in tissue coordinates, (B) colors in sPC space for spots, (C) colors in sPC space for reference single cells.
Annotate enriched genes along each sPC axis using EnrichR.
For B cells, the first sPC corresponds to an interferon response (Ifit1/3, Isg15, Oas3) concentrated in the interfollicular area of mycobacteria-exposed lymph nodes.
destvi_utils.explore_gamma_space(st_model, sc_model, ct_list=ct_list, ct_thresholds=ct_thresholds)
Genes associated with SpatialPC1
Positively
---------------------------------------------------------------------------------------
Ifit3, Stat1, Rtp4, Slfn5, Zbp1, Irf7, Trim30a, Usp18, Ifi47, Ifit3b, Parp14, Ifit1, Ifit2, Oasl2, Isg15, Ifi27l2a, Rnf213, Rsad2, Igtp, Gbp7, Serpina3g, Bst2, Phf11b, Oas3, Eif2ak2, Gbp4, Pkib, Isg20, Herc6, Ly6a, Lgals3bp, Cmpk2, Malat1, Ifitm3, Trafd1, Irf1, Cybb, Socs1, Ms4a4c, Sdc3, Ms4a4b, Plac8, Psmb9, Serpina3f, Gbp2, Epsti1, Mpeg1, Tcf4, Phf11a, B2m
---------------------------------------------------------------------------------------
Interferon signaling*, Interferon alpha/beta signaling*, Immune system signaling by interferons, interleukins, prolactin, and growth hormones*, Immune system*, Type II interferon signaling (interferon-gamma)*, Interferon-gamma signaling pathway*, Antiviral mechanism by interferon-stimulated genes*, Interferon alpha signaling regulation*, Toll-like receptor signaling pathway regulation*, Interferon gamma signaling regulation*
Negatively
---------------------------------------------------------------------------------------
Rpsa, Rps2, Rplp0, Rpl41, Ppia, Fth1, Rplp1, Tmsb4x, Ftl1, Ptma, Chchd10, Nme2, Actg1, Npm1, Rps12, Gm26532, Tubb4b, Calr, Hmgb1, Cd74, Atp5g1, Ybx1, Slc25a5, Calm1, Impdh2, Marcksl1, Ppp1r14b, Tuba1a, Pxdc1, Cd83, Mif, Nr4a1, Fdps, Pdia6, H2-Aa, Bri3, Tubb5, Junb, Hs3st1, Gcsh, Dctpp1, Plaur, Rexo2, H2-Ab1, Vpreb3, Prdx1, Tuba1b, Eif4a1, Ran, Fkbp2
---------------------------------------------------------------------------------------
Protein metabolism*, Disease*, Influenza infection*, T cell receptor regulation of apoptosis*, Translation*, Cytoplasmic ribosomal proteins*, Influenza viral RNA transcription and replication*, Post-chaperonin tubulin folding pathway*, Pathogenic Escherichia coli infection*, Activation of mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S*
Genes associated with SpatialPC2
Positively
---------------------------------------------------------------------------------------
Cd79a, Cd74, H2-Eb1, H2-Aa, H2-Ab1, Ighm, Ly6d, Fth1, Malat1, Mef2c, Ms4a1, Iglc3, Iglc2, Igkc, Fcer2a, Fcmr, Fcrla, Rnase6, Vpreb3, Ctsh, Mzb1, Itm2b, Napsa, Gm42418, Serpinb1a, S100a10, Lsp1, Crip1, Cst3, Srgn, Fcgrt, Hexb, Cd38, Bcl11a, Tspo, Ly86, Tmsb4x, Hist1h1c, Cd81, Rplp1, Jun, Ciita, Hspa1b, Wfdc17, Id3, Neat1, Vim, Pglyrp1, Capg, Unc93b1
---------------------------------------------------------------------------------------
Antigen-activated B-cell receptor generation of second messengers*, Immune system*, Innate immune system*, Signaling by the B cell receptor (BCR)*, Toll receptor cascades*, Adaptive immune system*, Immunoregulatory interactions between a lymphoid and a non-lymphoid cell*, Complement cascade, Antigen processing and presentation, Nuclear events mediated by MAP kinases
Negatively
---------------------------------------------------------------------------------------
Birc5, Ccna2, Tpx2, Ube2c, Rrm2, Mki67, Kif11, Hmmr, Cdk1, Ccnb2, Cdca3, Ncapg, Bub1, Kif4, Plk1, Ncaph, Cdca8, Kif2c, Pbk, Hist1h3c, Uhrf1, Esco2, Kif22, Neil3, Spc24, Cks1b, Ska1, Cdkn3, Cdca5, Cenpf, Ccnb1, Nusap1, Prc1, E2f8, Top2a, Ppp1r14b, Cdca2, Nuf2, Melk, Tyms, Kif15, Mxd3, Kif20a, Spag5, Cep55, Ccnd2, Ckap2l, Kif18b, Mcm10, Ddx21
---------------------------------------------------------------------------------------
Cell cycle*, M phase pathway*, FOXM1 transcription factor network*, DNA replication*, Polo-like kinase 1 (PLK1) pathway*, Aurora B signaling*, Mitotic G1-G1/S phases*, Cyclin A/B1-associated events during G2/M transition*, Kinesins*, Mitotic prometaphase*
Genes associated with SpatialPC1
Positively
---------------------------------------------------------------------------------------
Rgcc, Ccr9, Cd8b1, Ramp1, Cd3g, Npc2, Cd8a, Itgae, Igfbp4, Ppia, Lef1, Actn2, Cd3e, Trbc2, Actb, Lat, Cd3d, Inpp4b, Trac, Actn1, Cd226, Plekho1, Timp2, Sox4, Ifi27l2a, Tmsb4x, Ddit4, Dtx1, Gsn, Bcl11b, Gpr68, Tesc, Lztfl1, Ptma, Tubb5, Rplp0, Izumo1r, Malat1, Sh3bp1, Thy1, Itm2b, Lck, Rpl41, Trat1, Dapl1, Cyb5a, Cd247, Auts2, Tdrp, Gria3
---------------------------------------------------------------------------------------
T helper cell surface molecules*, Generation of second messenger molecules*, T cell receptor signaling in naive CD8+ T cells*, Lck and Fyn tyrosine kinases in initiation of T cell receptor activation*, PD-1 signaling*, Interleukin-17 signaling pathway*, MEF2D role in T cell apoptosis*, Inhibition of T cell receptor signaling by activated Csk*, T cell activation co-stimulatory signal*, HIV-induced T cell apoptosis*
Negatively
---------------------------------------------------------------------------------------
Ccl5, Ly6c2, Ctla2a, Xcl1, Cxcr3, Samd3, S100a6, 1700025G04Rik, Il18rap, Il2rb, Klrc2, Ccr5, Hopx, Ms4a4c, Slamf7, Klrc1, Cd44, Klre1, Ifngr1, Eomes, Plek, Lrrk1, Il18r1, Gzmm, Dennd4a, Ctla2b, Klrk1, Ahnak, Sidt1, Klra7, Serpina3g, Fcer1g, Ifitm10, Fasl, Anxa2, Ccr2, Pglyrp1, Fosb, Runx2, Itgb1, Gem, Il7r, Gpr183, Ccl4, Arl4d, AW112010, Tnfaip3, Spry2, Ifng, Pim1
---------------------------------------------------------------------------------------
Interleukin-2 signaling pathway*, Cytokine-cytokine receptor interaction*, Selective expression of chemokine receptors during T-cell polarization*, T cell receptor regulation of apoptosis*, Binding of chemokines to chemokine receptors*, Leptin influence on immune response*, Interleukin-12-mediated signaling events*, Interleukin-12/STAT4 pathway*, TNF-alpha effects on cytokine activity, cell motility, and apoptosis*, Natural killer cell-mediated cytotoxicity*
Genes associated with SpatialPC2
Positively
---------------------------------------------------------------------------------------
Isg15, Stat1, Ifit1, Ifit3, Zbp1, Igtp, Iigp1, Slfn5, Gbp2, Isg20, Gbp7, Rtp4, Rsad2, Oasl2, Ifit3b, Gbp4, Ifi27l2a, Ccr9, Usp18, Irf7, Ifi47, Ly6a, Parp14, Bst2, Rnf213, Gbp8, Slfn1, B2m, Cmpk2, Ddx60, Malat1, Trim30a, Herc6, Oas3, Cd274, Ifih1, Rgcc, Ramp1, Eif2ak2, Lgals3bp, Phf11b, Itgae, Trafd1, Gbp5, Pmepa1, AW011738, Ppa1, Psmb9, Art2b, Igfbp4
---------------------------------------------------------------------------------------
Interferon signaling*, Immune system signaling by interferons, interleukins, prolactin, and growth hormones*, Interferon alpha/beta signaling*, Immune system*, Interferon-gamma signaling pathway*, Antiviral mechanism by interferon-stimulated genes*, Type II interferon signaling (interferon-gamma)*, TRAF3-dependent IRF activation pathway, RIG-I-like receptor signaling pathway, RIG-I/MDA5-mediated induction of interferon-alpha/beta pathways
Negatively
---------------------------------------------------------------------------------------
Ccl5, Rpsa, Ctla2a, Ly6c2, Cxcr3, Rps12, Samd3, Hopx, Xcl1, Rplp0, Rps2, Il2rb, 1700025G04Rik, Cd44, Ahnak, Gzmm, Ifitm10, Klrd1, Pglyrp1, Plek, Ctla2b, Rpl41, Rplp1, Ifngr1, S100a6, H2afz, Il18rap, Zyx, Ccr5, Myo1f, Sidt1, Gpr183, Klra7, Bcl2, Lrrk1, Itgb1, Bbc3, Itm2a, Klre1, Nkg7, Il18r1, Klrc2, Crip1, Runx2, Anxa2, Eomes, Ccr2, Klrc1, Actg1, Npm1
---------------------------------------------------------------------------------------
Interleukin-2 signaling pathway*, Cytokine-cytokine receptor interaction*, T cell receptor regulation of apoptosis*, Binding of chemokines to chemokine receptors*, Cytoplasmic ribosomal proteins*, Ras-independent pathway in NK cell-mediated cytotoxicity*, Selective expression of chemokine receptors during T-cell polarization*, Interleukin-12-mediated signaling events*, Influenza viral RNA transcription and replication*, Influenza infection*
Genes associated with SpatialPC1
Positively
---------------------------------------------------------------------------------------
Ifit3, Cybb, Ifi204, Oas3, Plac8, Igkc, Ms4a6b, Rnf213, Rtp4, Ifit3b, Itgal, Ifit2, Oasl2, Ms4a1, Igtp, Rnase6, Cr2, Ms4a4c, Mzb1, Gbp2, Sp140, Iglc2, Pou2af1, Iglc3, Ifi47, Vpreb3, Cxcl10, Slfn5, Fcmr, Ifit1, Ly6d, Gnb4, Arhgap31, Fcgr1, Themis2, Icosl, Trim30a, Zbp1, Lsp1, Blk, Arhgef10l, Flt1, Kctd14, Cst3, Mef2c, Ighm, Fcrla, Gbp7, Smpdl3b, Hp
---------------------------------------------------------------------------------------
Immune system*, Interferon alpha/beta signaling*, Antigen-activated B-cell receptor generation of second messengers*, Interferon signaling*, Signaling by the B cell receptor (BCR)*, Immune system signaling by interferons, interleukins, prolactin, and growth hormones*, Immunoregulatory interactions between a lymphoid and a non-lymphoid cell*, Type II interferon signaling (interferon-gamma)*, B lymphocyte cell surface molecules*, Adaptive immune system
Negatively
---------------------------------------------------------------------------------------
Ppia, Fabp5, Cd63, Tuba1a, Cdk4, Sec61b, Mif, Tuba1c, Tmsb4x, Hebp1, Atpif1, Itgax, Spint1, Clec4n, Tspan3, Plxnc1, Rgs10, Dnase1l3, Hspe1, Eif4ebp1, Hnrnpll, Etfb, Pdia3, S100a1, Cfp, Sema7a, Itga8, Cyb5a, Clec4a2, Ftl1, Rps2, Basp1, Tgfbi, Atp5g1, Tuba1b, Sbf2, Uqcc2, Tfec, Cnih4, Ctsd, Tubb2a, Dbf4, Nhp2, Eif1ax, Lsm2, Kif2c, Anxa4, Hnrnpab, Hdgf, Tarm1
---------------------------------------------------------------------------------------
Post-chaperonin tubulin folding pathway*, Protein metabolism*, Cooperation of prefoldin and TriC/CCT in actin and tubulin folding*, Protein folding*, Pathogenic Escherichia coli infection*, Phagosome*, Gap junction pathway*, Activation of mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S*, Translation, Other semaphorin interactions
Genes associated with SpatialPC2
Positively
---------------------------------------------------------------------------------------
Tmem26, Vcam1, Chst2, Cd38, Gfra2, Cd81, Prkar1b, Etv5, Ly75, Prg3, Gpr157, Phyhd1, Rai14, Actn1, Tspan4, Sirpa, Asb2, Nuak1, Fam20c, Inpp4b, Gpm6b, Timp2, C1qc, Cacna1b, Cd63, Ehf, C1qa, Art2b, Mertk, C1qb, Lrp8, Pdgfa, Hebp1, Pltp, Zc3h12c, Nav2, Sucnr1, Itga9, Jup, Tspan3, Nr1h3, Kcnj10, Sema6d, Nrgn, Chek1, Il1r1, Ifng, 3830403N18Rik, Hbegf, Tgm1
---------------------------------------------------------------------------------------
Complement activation, classical pathway*, Hypertrophy pathway*, Systemic lupus erythematosus*, TGF-beta regulation
of extracellular matrix*, Prion diseases*, Chagas disease*, Hemostasis pathway*, Malaria, Mechanism of gene
regulation by peroxisome proliferators via PPAR-alpha, Interleukin-4 regulation of apoptosis
Negatively
---------------------------------------------------------------------------------------
Lsp1, Ms4a6c, Sap30, Ms4a6b, Fcer1g, Ms4a4c, Gpr141, S100a6, Ifitm3, Ifi205, Plac8, Napsa, Hp, Pnp, Oas3, Actb, Isg15, S100a11, Klrk1, Fgl2, Nme2, Ifi47, Tyrobp, Psme2, Sp140, Slfn1, Lmnb1, Klra2, Glrx, BC028528, Itgal, Pgk1, Anxa1, Ly6c2, Rassf4, Rtp4, Pid1, Wfdc17, Nupr1, Zyx, Anxa2, Emb, S100a4, Herc6, Cst3, Igtp, Zbp1, Srgn, Prdx5, Ms4a6d
---------------------------------------------------------------------------------------
T cell receptor regulation of apoptosis*, FSH regulation of apoptosis*, Prostaglandin biosynthesis and regulation*,
Interleukin-2 signaling pathway*, Natural killer cell-mediated cytotoxicity, Interferon alpha/beta signaling,
Nucleotide metabolism, Nucleotide di- and triphosphate biosynthesis and interconversion, Immune system,
Immunoregulatory interactions between a lymphoid and a non-lymphoid cell
The spatially-weighted PCA and its functional annotation provide a systematic way to formulate hypotheses about cell-state variation in spatial context — e.g., which sub-states are concentrated in which tissue compartments and how they change between conditions.
Step 4b — Pathway/signature scoring with VISION#
DestVI’s cell-type proportions and gamma values tell us which cell types are
present in each spot and how their internal state varies, but not what
pathways are active there. We use scviva.tools.vision.VisionAnalysis — the
same standalone VISION pipeline demonstrated in
Vision_tutorial — to score every spot against the
same KEGG metabolic pathway gene sets used later in Step 6’s Harreman
metabolic-module analysis, then ask which pathways track which cell-type
proportions. This is a complementary, curated-pathway-first counterpart to
Step 6’s data-driven Hotspot module discovery.
import requests
from scviva.tools.vision import VisionAnalysis
# st_adata.X is already log-normalized (see Data loading and preprocessing above),
# so norm_data_key=None tells VisionAnalysis to score directly from adata.X.
va = VisionAnalysis(st_adata, norm_data_key=None)
va.setup(compute_neighbors_on_key="spatial", num_neighbors=5)
resp = requests.get("https://exampledata.scverse.org/scvi-tools/KEGG_metab_m.json")
resp.raise_for_status()
kegg_metab = resp.json()
print(f"KEGG metabolic pathways available: {len(kegg_metab)}")
va.load_signatures(dicts=[kegg_metab], min_signature_genes=5, sig_gene_threshold=0.001)
va.compute_signatures()
autocorr_sig = st_adata.uns["vision_signature_scores"].sort_values("c_prime", ascending=False)
print(f"vision_signatures: {st_adata.obsm['vision_signatures'].shape}")
print("\nTop 10 most spatially autocorrelated KEGG pathways:")
print(autocorr_sig.head(10).round(4).to_string())
INFO Building KNN graph: n_cells=1092, K=5.
2026-09-07 14:43:54 | [INFO] Building KNN graph: n_cells=1092, K=5.
INFO KNN weight matrix built — nnz=5460.
2026-09-07 14:43:54 | [INFO] KNN weight matrix built — nnz=5460.
INFO KNN weights stored in adata.obsp['weights'] — shape (1092, 1092).
2026-09-07 14:43:54 | [INFO] KNN weights stored in adata.obsp['weights'] — shape (1092, 1092).
INFO VisionAnalysis: setup complete.
2026-09-07 14:43:54 | [INFO] VisionAnalysis: setup complete.
KEGG metabolic pathways available: 69
INFO Computing signatures for each cell...
2026-09-07 14:43:55 | [INFO] Computing signatures for each cell...
INFO Signatures computed in 0.016 s.
/home/access/miniconda3/lib/python3.13/site-packages/scviva/tools/vision/tools/signature.py:244: UserWarning: Sparse invariant checks are implicitly disabled. Memory errors (e.g. SEGFAULT) will occur when operating on a sparse tensor which violates the invariants, but checks incur performance overhead. To silence this warning, explicitly opt in or out. See `torch.sparse.check_sparse_tensor_invariants.__doc__` for guidance. (Triggered internally at /__w/pytorch/pytorch/aten/src/ATen/Context.cpp:848.)
X_pt = torch.sparse_csr_tensor(
/home/access/miniconda3/lib/python3.13/site-packages/scviva/tools/vision/tools/signature.py:244: UserWarning: Sparse CSR tensor support is in beta state. If you miss a functionality in the sparse tensor support, please submit a feature request to https://github.com/pytorch/pytorch/issues. (Triggered internally at /__w/pytorch/pytorch/aten/src/ATen/SparseCsrTensorImpl.cpp:49.)
X_pt = torch.sparse_csr_tensor(
2026-09-07 14:43:55 | [INFO] Signatures computed in 0.016 s.
INFO Computing Geary's C for all the signatures...
2026-09-07 14:43:55 | [INFO] Computing Geary's C for all the signatures...
INFO Generating the null distribution...
2026-09-07 14:43:55 | [INFO] Generating the null distribution...
INFO Computing signatures for each cell...
2026-09-07 14:43:55 | [INFO] Computing signatures for each cell...
INFO Signatures computed in 0.215 s.
2026-09-07 14:43:55 | [INFO] Signatures computed in 0.215 s.
INFO Finished computing Geary's C for all the signatures in 0.532 seconds.
2026-09-07 14:43:55 | [INFO] Finished computing Geary's C for all the signatures in 0.532 seconds.
vision_signatures: (1092, 19)
Top 10 most spatially autocorrelated KEGG pathways:
c_prime pvals fdr
GLYCOSAMINOGLYCAN DEGRADATION 0.2197 0.0696 0.8231
OXIDATIVE PHOSPHORYLATION 0.1882 0.0866 0.8231
GLYCOLYSIS GLUCONEOGENESIS 0.1533 0.1420 0.8990
ETHER LIPID METABOLISM 0.1088 0.2632 0.9900
PYRIMIDINE METABOLISM 0.0987 0.3302 0.9900
GLYCEROLIPID METABOLISM 0.0955 0.3229 0.9900
PURINE METABOLISM 0.0825 0.4212 0.9900
GLYCEROPHOSPHOLIPID METABOLISM 0.0788 0.4095 0.9900
AMINO SUGAR AND NUCLEOTIDE SUGAR METABOLISM 0.0685 0.4832 0.9900
ONE CARBON POOL BY FOLATE 0.0545 0.5675 0.9900
Correlating signature scores with cell-type proportions#
Cell-type proportions are compositional — they sum to 1 per spot and aren’t normally distributed — so we use Spearman rank correlation rather than Pearson between each KEGG pathway’s per-spot VISION score and each cell type’s per-spot proportion.
import pandas as pd
from scipy.stats import spearmanr
sig_scores = st_adata.obsm["vision_signatures"]
proportions = st_adata.obsm["proportions"]
corr = pd.DataFrame(index=sig_scores.columns, columns=proportions.columns, dtype=float)
for ct in proportions.columns:
for sig in sig_scores.columns:
corr.loc[sig, ct] = spearmanr(sig_scores[sig], proportions[ct]).statistic
# Keep only the pathways most strongly associated with at least one cell type
top_sigs = corr.abs().max(axis=1).sort_values(ascending=False).head(20).index
plt.figure(figsize=(8, 8))
sns.heatmap(corr.loc[top_sigs].astype(float), cmap="RdBu_r", center=0, vmin=-0.5, vmax=0.5)
plt.title("Spearman correlation: KEGG pathway score vs. cell-type proportion")
plt.tight_layout()
plt.show()
Per-cell-type pathway enrichment#
destvi_utils.automatic_proportion_threshold already gave us a data-driven cutoff (ct_thresholds) separating spots where a cell type is genuinely present from spots where its estimated proportion is residual noise (Step 3). We use it to test (per KEGG pathway, per cell type in ct_list) whether the pathway’s VISION score is significantly higher in spots where the cell type is present (one-sided Mann-Whitney U, BH-corrected).
from scipy.stats import mannwhitneyu
from statsmodels.stats.multitest import multipletests
enrichment_records = []
for ct in ct_list:
present = proportions[ct].values > ct_thresholds[ct]
for sig in sig_scores.columns:
stat, pval = mannwhitneyu(
sig_scores.loc[present, sig], sig_scores.loc[~present, sig], alternative="greater"
)
enrichment_records.append({"cell_type": ct, "signature": sig, "pval": pval})
enrichment = pd.DataFrame(enrichment_records)
enrichment["fdr"] = multipletests(enrichment["pval"], method="fdr_bh")[1]
print("Top KEGG pathways enriched in cell-type-present spots (FDR < 0.05):")
print(
enrichment[enrichment["fdr"] < 0.05]
.sort_values(["cell_type", "fdr"])
.groupby("cell_type")
.head(3)
.to_string(index=False)
)
Top KEGG pathways enriched in cell-type-present spots (FDR < 0.05):
cell_type signature pval fdr
B cells GLYCOSAMINOGLYCAN DEGRADATION 2.484476e-14 1.416152e-12
B cells ETHER LIPID METABOLISM 1.655367e-10 4.717796e-09
B cells AMINO SUGAR AND NUCLEOTIDE SUGAR METABOLISM 2.497193e-07 3.558501e-06
CD8 T cells GLYCOLYSIS GLUCONEOGENESIS 2.087917e-07 3.558501e-06
CD8 T cells ONE CARBON POOL BY FOLATE 3.093918e-06 3.527067e-05
CD8 T cells FRUCTOSE AND MANNOSE METABOLISM 1.040005e-04 8.908317e-04
Step 5 — Cell-type-specific differential expression (B cells)#
We focus on B cells, which the gamma-space analysis flagged as spatially variable for interferon response genes. DestVI can impute cell-type-specific gene expression for any spot using get_scale_for_ct. This integrates proportion weights and gamma values to estimate what the B-cell-specific transcriptome looks like in each spot — even in spots that are not exclusively B cells.
Below, we visualize the spatial distribution of IFN-response genes (Ifit3, Ifit1, Isg15, etc.) restricted to spots with ≥20% B cells.
plt.figure(figsize=(8, 8))
ct_name = "B cells"
gene_name = ["Ifit3", "Ifit3b", "Ifit1", "Isg15", "Oas3", "Usp18", "Isg20"]
# Restrict to spots with at least 20% B cells to reduce imputation noise
indices = np.where(st_adata.obsm["proportions"][ct_name].values > 0.2)[0]
# Impute B-cell-specific expression and sum across IFN genes
specific_expression = np.sum(st_model.get_scale_for_ct(ct_name, indices=indices)[gene_name], 1)
specific_expression = np.log(1 + 1e4 * specific_expression)
# Plot all spots as faint background; color foreground spots by IFN expression
plt.scatter(st_adata.obsm["location"][:, 0], st_adata.obsm["location"][:, 1], alpha=0.05)
plt.scatter(
st_adata.obsm["location"][indices][:, 0],
st_adata.obsm["location"][indices][:, 1],
c=specific_expression,
s=10,
cmap="Reds",
)
plt.colorbar()
plt.title(f"Imputation of {gene_name} in {ct_name}")
plt.show()
We apply a Kolmogorov-Smirnov test on imputed B-cell expression to compare two spatial regions in the treated lymph nodes:
IFN-rich zone: treated sections (TC or BD) where
log(1 + 1e5 × Ifit3)exceeds a threshold of 4.Comparison zone: same sections but below the threshold.
The volcano plot highlights IFN-response genes (Ifit3, Isg15, Usp18, etc.) as the top upregulated genes in the IFN-rich interfollicular area, consistent with the spatially-weighted PCA result.
ct = "B cells"
imputation = st_model.get_scale_for_ct(ct)
color = np.log(1 + 1e5 * imputation["Ifit3"].values)
threshold = 4 # separates IFN-high from IFN-low B-cell spots
# IFN-rich zone: treated sections (TC or BD) with high Ifit3 expression
mask = np.logical_and(
np.logical_or(st_adata.obs["LN"] == "TC", st_adata.obs["LN"] == "BD"),
color > threshold,
).values
# Comparison zone: same sections but low Ifit3 expression
mask2 = np.logical_and(
np.logical_or(st_adata.obs["LN"] == "TC", st_adata.obs["LN"] == "BD"),
color < threshold,
).values
# Run KS test on imputed B-cell expression; results stored in st_adata.uns["IFN_rich"]
#_ = destvi_utils.de_genes(
# st_model, mask=mask, mask2=mask2, threshold=ct_thresholds[ct], ct=ct, key="IFN_rich"
#)
_ = destvi_utils.de_genes(
st_model, mask=mask, mask2=mask2, threshold=ct_thresholds[ct], ct=ct, key="IFN_rich",
st_adata=st_adata,
)
display(st_adata.uns["IFN_rich"]["de_results"].head(10))
| log2FC | score | pval | |
|---|---|---|---|
| Ifit3b | 4.266329 | 0.740779 | 0.0 |
| Ifit3 | 4.186249 | 0.735299 | 0.0 |
| Ifit1 | 3.779947 | 0.723096 | 0.0 |
| Rsad2 | 3.556458 | 0.632184 | 0.0 |
| Ifit2 | 3.078209 | 0.629725 | 0.0 |
| Usp18 | 3.230167 | 0.621738 | 0.0 |
| Irf7 | 2.517198 | 0.570902 | 0.0 |
| Isg15 | 2.749901 | 0.570579 | 0.0 |
| Ifitm3 | 2.729520 | 0.564957 | 0.0 |
| Oas3 | 3.011421 | 0.563694 | 0.0 |
destvi_utils.plot_de_genes(
st_adata,
interesting_genes=["Ifit3", "Ifit3b", "Ifit1", "Isg15", "Oas3", "Usp18", "Isg20"],
key="IFN_rich",
)
Step 6 — Cell-type-aware metabolic crosstalk with Harreman#
Having deconvolved the spatial data with DestVI, we now use Harreman to infer cell-type-aware metabolic cell-cell communication (CCC). When a DestVI model is passed to HarremanAnalysis, Harreman automatically:
Retrieves cell-type proportions via
model.get_proportions().Constructs proportion-weighted, cell-type-specific expression layers for every cell type.
Enables
ct_specific=Truemode, scoring interactions separately for each sending/receiving cell-type combination.
This contrasts with the cell-type-agnostic mode used in the Visium colon tutorial, where all expressed genes contribute equally regardless of cell type.
import numpy as np
from scipy.stats import zscore
from scviva.tools import HarremanAnalysis
# Assign the dominant cell type to each spot using Z-normalized proportions.
# Z-normalization is important: without it, abundant types (e.g., B cells) would
# always dominate regardless of local enrichment in a given spot.
proportions = st_model.get_proportions()
z_proportions = proportions.apply(zscore, axis=0)
st_adata.obs["dominant_cell_type"] = z_proportions.idxmax(axis=1)
st_adata.obs["dominant_cell_type"].value_counts()
dominant_cell_type
B cells 330
Migratory DCs 101
CD8 T cells 92
cDC1s 91
GD T cells 79
pDCs 76
cDC2s 71
CD4 T cells 65
Monocytes 63
Tregs 63
Macrophages 36
NK cells 25
Name: count, dtype: int64
# Passing model=st_model triggers the DestVI integration path:
# - get_proportions() is called automatically
# - proportion-weighted, cell-type-specific expression layers are attached to adata
# is_deconvolved=True confirms the DestVI layers were set up successfully
ha = HarremanAnalysis(st_adata, model=st_model)
ha.setup(
compute_neighbors_on_key="spatial",
species="mouse",
database="both", # HarremanDB (transporters) + CellChatDB (LR pairs)
n_neighbors=5,
cell_type_key="dominant_cell_type",
)
print(f"is_deconvolved={ha.is_deconvolved}")
print(ha)
INFO HarremanAnalysis: attached 12 DestVI cell-type layers.
2026-09-07 14:43:59 | [INFO] HarremanAnalysis: attached 12 DestVI cell-type layers.
INFO HarremanAnalysis: setup complete (species=mouse, database=both).
2026-09-07 14:45:17 | [INFO] HarremanAnalysis: setup complete (species=mouse, database=both).
is_deconvolved=True
HarremanAnalysis(is_set_up=True, is_deconvolved=True, completed_steps=['setup'])
# ct_specific=True: enumerate gene pairs for each sender/receiver cell-type combination
ha.compute_gene_pairs(ct_specific=True)
# Infer cell-type-aware metabolic CCC:
# mode="cell_type" uses proportion-weighted expression layers
# Both parametric (DANB) and non-parametric (1000-permutation) tests are run
ha.compute_cell_communication(mode="cell_type", n_permutations=1000, test="both")
/home/access/miniconda3/lib/python3.13/site-packages/scviva/tools/harreman/tools/cell_communication/_ct_ccc_scoring.py:657: UserWarning: Creating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with numpy.array() before converting to a tensor. (Triggered internally at /__w/pytorch/pytorch/torch/csrc/utils/tensor_new.cpp:252.)
indices = torch.tensor([weights.row, weights.col], dtype=torch.long, device=device)
Permutation test: 100%|██████████████████████████████████████████████████████████████████████████████████████| 1000/1000 [00:21<00:00, 46.86it/s]
assert ha.results.ct_cell_communication is not None # ha.results reflects the real compute_cell_communication output
# Retain interactions with FDR < 0.05 (non-parametric test by default)
ha.select_significant_interactions(fdr_threshold=0.05)
Spatial autocorrelation and metabolic module discovery#
We run the Hotspot pipeline to identify spatially co-varying metabolic gene modules. The resulting module scores are used later to interpret which metabolic zones are linked to specific cell-type-aware interactions.
# Compute spatial autocorrelation restricted to mouse metabolic enzymes (DANB model)
ha.hs.compute_local_autocorrelation(
layer_key="counts", model="danb", species="mouse", use_metabolic_genes=True
)
/home/access/miniconda3/lib/python3.13/site-packages/scipy/sparse/_index.py:216: SparseEfficiencyWarning: Changing the sparsity structure of a csr_array is expensive. lil and dok are more efficient.
self._set_arrayXarray(i, j, x)
# Pairwise local correlation on all spatially autocorrelated metabolic genes
ha.hs.compute_local_correlation()
/home/access/miniconda3/lib/python3.13/site-packages/scipy/sparse/_index.py:216: SparseEfficiencyWarning: Changing the sparsity structure of a csr_array is expensive. lil and dok are more efficient.
self._set_arrayXarray(i, j, x)
# Cluster genes into metabolic modules using agglomerative clustering
ha.hs.create_modules(min_gene_threshold=10)
# Compute per-spot module scores; device="cpu" avoids GPU memory conflicts
ha.hs.calculate_module_scores(device="cpu")
0%| | 0/2 [00:00<?, ?it/s]/home/access/miniconda3/lib/python3.13/site-packages/scipy/sparse/_index.py:216: SparseEfficiencyWarning: Changing the sparsity structure of a csr_array is expensive. lil and dok are more efficient.
self._set_arrayXarray(i, j, x)
/home/access/miniconda3/lib/python3.13/site-packages/scipy/sparse/_index.py:216: SparseEfficiencyWarning: Changing the sparsity structure of a csr_array is expensive. lil and dok are more efficient.
self._set_arrayXarray(i, j, x)
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 70.17it/s]
ha.pl.local_correlation_plot(mod_cmap="tab10")
Correlate interaction scores with metabolic module activity#
We compute per-spot interaction scores in cell-type-aware mode (mode='cell_type') and Pearson-correlate them with the metabolic module scores discovered above. A high correlation indicates that a specific cell-type-aware metabolic exchange preferentially occurs in spots dominated by a particular metabolic gene program.
# Compute per-spot interaction scores in cell-type-aware mode
# Both parametric and non-parametric scores are computed
ha.compute_interacting_cell_scores(
mode="cell_type", test="both", device="cpu", n_permutations=1000
)
# Pearson-correlate metabolite interaction scores with module scores
# ct_aware=True uses the cell-type-specific interaction scores
ha.tl.compute_interaction_module_correlation(
cor_method="pearson",
interaction_type="metabolite",
test="non-parametric",
ct_aware=True,
)
ha.pl.plot_interaction_module_correlation(threshold=0.1, figsize = (10, 80))
import gc
gc.collect()
import torch
torch.cuda.empty_cache()