DOTpy

July 13, 2026 · View on GitHub

Deconvolution by Optimal Transport for Spatial Transcriptomics

A Python implementation of the DOT algorithm for transferring cell type annotations from single-cell RNA-seq reference data to spatial transcriptomics data using multi-objective optimization.

Features

🚀 GPU acceleration using PyTorch for fast computation
🧬 AnnData integration - seamlessly works with scanpy workflows
🎯 Multi-objective optimization using Frank-Wolfe algorithm
📊 High & low resolution support - works with both subcellular (Xenium, MERFISH, CosMx) and spot-based (Visium, ST) technologies
🎨 Built-in visualization tools for spatial cell type mapping
💾 Checkpointing for long-running optimizations
⚡ Mixed precision support for memory-efficient GPU training

Installation

From source

git clone https://github.com/earmingol/DOTpy.git
cd DOTpy
pip install -e .

Requirements

Python >= 3.8
PyTorch >= 1.10.0 (with CUDA support for GPU acceleration)
scanpy >= 1.9.0
anndata >= 0.8.0
numpy >= 1.20.0
matplotlib >= 3.5.0
scikit-learn >= 1.0.0
scipy >= 1.7.0

Quick Start

Basic Usage

import scanpy as sc
from dotpy import DOT, setup_reference, setup_spatial, plot_spatial_weights

# Load data
ref_adata = sc.read_h5ad('reference.h5ad')
spatial_adata = sc.read_h5ad('spatial.h5ad')

# Process reference and spatial data
ref_processed = setup_reference(
    ref_adata,
    cell_type_key='cell_type',
    subcluster_size=10,
    max_genes=5000,
    verbose=True
)

spatial_processed = setup_spatial(
    spatial_adata,
    spatial_key='spatial',
    th_spatial=0.84,
    verbose=True
)

# Run DOT with batching
dot = DOT(
    spatial_processed, 
    ref_processed,
    batch_size=500  # Adjust for your GPU memory
)

dot.fit(
    mode='highres',
    iterations=100,
    checkpoint_dir='./checkpoints',  # Save checkpoints
    checkpoint_freq=10,
    verbose=True
)

# Get results
weights = dot.get_weights(normalize=True)
cell_types = dot.get_cell_types()

# Visualize results
plot_spatial_weights(
    spatial_adata.obsm['spatial'],
    weights,
    cell_types=cell_types,
    ncols=4,
    save_path='cell_type_maps.png'
)

Resume from Checkpoint

dot.fit(
    mode='highres',
    iterations=100,
    resume_from='./checkpoints/checkpoint_iter_50.pkl',
    verbose=True
)

Command-Line Interface

For production workflows or batch processing, DOTpy includes a CLI script (run_dot_cli.py) that runs the full pipeline from the terminal without writing any Python code.

Basic usage

python run_dot_cli.py --ref reference.h5ad --spatial spatial.h5ad

When a single AnnData object contains multiple tissue sections or slides, the CLI can iterate through each sample automatically using --sample-key. Each sample is preprocessed, deconvolved, and saved independently, with GPU memory freed between runs:

python run_dot_cli.py \
    --ref reference.h5ad \
    --spatial spatial_multi_slide.h5ad \
    --sample-key slide_id \
    --save-combined \
    -v

This produces per-sample results (weights.csv, annotations.csv, plots) and optionally a combined output with --save-combined.

Resolution modes

# High-resolution (Xenium, MERFISH, CosMx) — default
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad --mode highres

# Low-resolution (Visium, ST)
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad --mode lowres --ratios-weight 0.3

GPU acceleration and memory options

# Automatic GPU detection (default)
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad --device auto

# Force CPU
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad --device cpu

# Mixed precision for large datasets on GPU
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad --device cuda --mixed-precision

Checkpointing for long runs

python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad \
    --checkpoint-dir ./checkpoints --checkpoint-freq 10

# Resume from a checkpoint
python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad \
    --resume-from ./checkpoints/sample_1/checkpoint_iter_50.pkl

Lineage-level annotation

When the reference contains a higher-level grouping (e.g., lineage or broad class), the CLI can map cell types to that level automatically:

python run_dot_cli.py --ref ref.h5ad --spatial spatial.h5ad \
    --cell-type-key cell_type --lineage-key lineage

Full example

python run_dot_cli.py \
    --ref reference.h5ad \
    --spatial visium_slides.h5ad \
    --sample-key slide_id \
    --cell-type-key cell_subclass \
    --lineage-key cell_class \
    --mode lowres \
    --ratios-weight 0.3 \
    --max-genes 5000 \
    --subcluster-size 10 \
    --batch-size 5000 \
    --iterations 100 \
    --device auto \
    --output dot_results \
    --output-dir ./results \
    --save-combined \
    -v

All CLI options

Flag	Default	Description
`--ref`	(required)	Path to reference scRNA-seq h5ad file
`--spatial`	(required)	Path to spatial transcriptomics h5ad file
`--sample-key`	`None`	Column in `obs` to split spatial data by sample/slide
`--cell-type-key`	`cell_type`	Column in reference `obs` with cell type labels
`--lineage-key`	`None`	Column in reference `obs` for higher-level grouping
`--counts-layer`	`counts`	Layer in spatial data with raw counts (`"X"` to use `.X` directly)
`--ref-counts-layer`	`X`	Layer in reference data with raw counts (`"X"` to use `.X` directly)
`--mode`	`highres`	`highres` or `lowres`
`--ratios-weight`	`0.0`	Weight for matching reference cell-type abundances
`--max-genes`	`5000`	Maximum genes for reference preprocessing
`--subcluster-size`	`10`	Maximum subclusters per cell type
`--th-spatial`	`0.84`	Cosine similarity threshold for spatial pairs
`--batch-size`	`5000`	Batch size for GPU processing
`--iterations`	`100`	Maximum Frank-Wolfe iterations
`--device`	`auto`	`auto`, `cuda`, or `cpu`
`--mixed-precision`	off	Use float16 intermediates on GPU
`--checkpoint-dir`	`None`	Directory for checkpoints
`--checkpoint-freq`	`10`	Checkpoint every N iterations
`--resume-from`	`None`	Resume from a checkpoint file
`--output`	`dot_results`	Output file prefix
`--output-dir`	`.`	Output directory
`--save-combined`	off	Merge per-sample results into one file
`--no-plots`	off	Skip plot generation
`--no-h5ad`	off	Skip saving per-sample h5ad files
`-v`, `--verbose`	off	Print detailed progress

Output files

For each sample, the CLI produces:

results/
  dot_results_slide1_weights.csv         # Cell-type weights per spot (S x K)
  dot_results_slide1_annotations.csv     # Dominant cell type per spot
  dot_results_slide1.h5ad                # Full AnnData with results
  figures/
    dot_results_slide1_cell_types.png    # Spatial cell type map
    dot_results_slide1_weights.png       # Per-type weight heatmaps
    dot_results_slide1_convergence.png   # Optimization convergence plot
  dot_results_combined_weights.csv       # (with --save-combined)
  dot_results_combined.h5ad              # (with --save-combined)

High-Resolution Data (Xenium, MERFISH, CosMx)

For subcellular resolution data where each spot typically contains 1 cell:

dot.fit(
    mode='highres',
    ratios_weight=0.0,
    iterations=100,
    verbose=True
)

Low-Resolution Data (Visium, ST)

For spot-based technologies where spots contain multiple cells:

dot.fit(
    mode='lowres',
    max_spot_size=20,  # Maximum cells per spot
    ratios_weight=0.3,  # Weight for matching cell type proportions
    iterations=100,
    verbose=True
)

Algorithm Overview

DOT uses multi-objective optimization to find cell type assignments that:

Match gene expression - Predicted expression should match observed spatial data
Preserve spatial coherence - Neighboring spots should have similar composition
Respect cell type abundances - Overall proportions should match reference (optional)
Enforce sparsity - Limit mixing of cell types per spot

The optimization is performed using the Frank-Wolfe algorithm, which efficiently handles the constrained optimization problem on GPUs.

Advanced Usage

Custom Parameters

# Setup reference with custom parameters
ref_processed = setup_reference(
    ref_adata,
    cell_type_key='cell_type',
    subcluster_size=15,      # More subclusters per cell type
    max_genes=10000,         # Use more genes
    remove_mt=True,          # Remove mitochondrial genes
    th_inner_logfold=0.75,   # Log-fold threshold for gene selection
    random_state=42,         # For reproducibility
    verbose=True
)

# Setup spatial with custom thresholds
spatial_processed = setup_spatial(
    spatial_adata,
    spatial_key='spatial',
    th_spatial=0.80,         # Adjust spatial similarity threshold
    th_gene_low=0.01,        # Minimum gene expression frequency
    th_gene_high=0.99,       # Maximum gene expression frequency
    radius='auto',           # Or specify numeric value
    remove_mt=True,          # Remove mitochondrial genes
    verbose=True
)

# DOT with custom device and optimization settings
import torch
device = 'cuda' if torch.cuda.is_available() else 'cpu'

dot = DOT(
    spatial_processed,
    ref_processed,
    batch_size=500,          # Adjust for GPU memory
    device=device            # Explicitly set device
)

# Fine-tune optimization
dot.fit(
    mode='highres',
    ratios_weight=0.2,       # Weight for abundance matching
    iterations=200,          # More iterations
    gap_threshold=0.001,     # Tighter convergence
    use_mixed_precision=True,  # Use float16 on GPU
    verbose=True
)

GPU/CPU Selection

# Check if CUDA is available
import torch
device = 'cuda' if torch.cuda.is_available() else 'cpu'
print(f"Using device: {device}")

# Pass device to DOT
dot = DOT(
    spatial_processed,
    ref_processed,
    device=device
)

Saving Results

# Add results to spatial AnnData
spatial_adata.obsm['dot_weights'] = weights

# Add individual cell type columns
for i, ct in enumerate(cell_types):
    spatial_adata.obs[f'dot_{ct}'] = weights[:, i]

# Save
spatial_adata.write('spatial_with_deconvolution.h5ad')

Visualization

Spatial Cell Type Maps

from dotpy.visualization import plot_spatial_weights

fig = plot_spatial_weights(
    coords=spatial_adata.obsm['spatial'],
    weights=weights,
    cell_types=cell_types,
    ncols=4,
    point_size=10,
    cmap='magma',
    flip_y=True,
    save_path='cell_type_maps.png',
    dpi=300
)

Optimization History

from dotpy.visualization import plot_optimization_history

fig = plot_optimization_history(
    dot.history,
    save_path='optimization_history.png'
)

Cell Type Proportions

from dotpy.visualization import plot_cell_type_proportions

fig = plot_cell_type_proportions(
    weights,
    cell_types=cell_types,
    save_path='proportions.png'
)

Performance Tips

GPU Acceleration

DOTpy automatically uses CUDA if available. For best performance:

# Check GPU memory
import torch
if torch.cuda.is_available():
    print(f"GPU: {torch.cuda.get_device_name(0)}")
    print(f"Memory: {torch.cuda.get_device_properties(0).total_memory / 1e9:.2f} GB")

Memory Management

For very large datasets:

# Reduce number of genes
ref_processed = setup_reference(
    ref_adata,
    max_genes=2000,  # Use fewer genes
    ...
)

# Reduce subclustering
ref_processed = setup_reference(
    ref_adata,
    subcluster_size=5,  # Fewer subclusters
    ...
)

# Use smaller batch size
dot = DOT(spatial, ref, batch_size=100)

# Enable mixed precision on GPU
dot.fit(
    mode='highres',
    use_mixed_precision=True,
    iterations=100
)

Speed vs Accuracy

# Faster (fewer iterations)
dot.fit(mode='highres', iterations=50)

# More accurate (more iterations, tighter convergence)
dot.fit(
    mode='highres',
    iterations=200,
    gap_threshold=0.001
)

Comparison with R Implementation

This PyTorch implementation provides:

✅ Faster computation through GPU acceleration
✅ Same algorithm and mathematical formulation
✅ AnnData integration for Python/scanpy workflows
✅ Memory efficiency through PyTorch's optimized operations

Key differences:

Uses PyTorch tensors instead of R matrices
Integrates with scanpy/AnnData instead of Seurat
Supports GPU acceleration out of the box

Citation

If you use DOT in your research, please cite:

Rahimi, A., Vale-Silva, L.A., Fälth Savitski, M. et al. 
DOT: a flexible multi-objective optimization framework for transferring features across single-cell and spatial omics. 
Nat Commun 15, 4994 (2024). https://doi.org/10.1038/s41467-024-48868-z

Troubleshooting

"CUDA out of memory"

# Solution 1: Reduce batch size
dot = DOT(spatial, ref, batch_size=100)

# Solution 2: Enable mixed precision
dot.fit(mode='highres', use_mixed_precision=True)

# Solution 3: Use CPU
dot = DOT(spatial, ref, device='cpu')

"Too slow on CPU"

# Solution: Reduce data size
ref = setup_reference(adata, max_genes=2000, subcluster_size=5)

No common genes found

# Check gene names
print(f"Ref genes: {ref_adata.var_names[:10]}")
print(f"Spatial genes: {spatial_adata.var_names[:10]}")

# Ensure gene names match (e.g., both use same gene ID system)

"Convergence issues"

# Solution: More iterations or looser threshold
dot.fit(iterations=200, gap_threshold=0.05)