Getting Started: Single-Slice Analysis

This tutorial learns the probabilistic relationship between continuous cell states and local tissue contexts in a deterministic synthetic dataset. It then uses that relationship to identify niches, compare cell-state-dependent niche preferences, and estimate the probability of a target population in each context. The self-contained script run_tutorial.py demonstrates the same core workflow.

1. Create spatial data

Mievformer requires a continuous cell-state representation in adata.obsm['X_pca'] and spatial coordinates in adata.obsm['spatial']. The tutorial generator creates 120 cells with three spatially organized cell types. The labels are used for interpretation and population selection after training; the model learns from the continuous molecular states.

import numpy as np
import pandas as pd
import scanpy as sc

rng = np.random.default_rng(7)
side_x, side_y = 12, 10
spatial = np.stack(
    np.meshgrid(np.arange(side_x), np.arange(side_y)), axis=-1
).reshape(-1, 2)
spatial = spatial + rng.normal(scale=0.08, size=spatial.shape)
cell_type_index = np.digitize(spatial[:, 0], bins=[4.0, 8.0])
cell_types = np.asarray(['TypeA', 'TypeB', 'TypeC'])[cell_type_index]

means = np.full((spatial.shape[0], 30), 1.2)
for group in range(3):
    means[cell_type_index == group, group * 5:(group + 1) * 5] += 3.0
counts = rng.poisson(means).astype(np.float32)
u, s, _ = np.linalg.svd(
    counts - counts.mean(axis=0, keepdims=True), full_matrices=False
)

adata = sc.AnnData(
    counts,
    obs=pd.DataFrame(
        {'cell_type': cell_types},
        index=[f'cell_{i:03d}' for i in range(spatial.shape[0])],
    ),
)
adata.var_names = [f'Gene_{i:02d}' for i in range(adata.n_vars)]
adata.obsm['spatial'] = spatial.astype(np.float32)
adata.obsm['X_pca'] = (u[:, :10] * s[:10]).astype(np.float32)

For real data, replace this block with your AnnData object while keeping the same two obsm keys.

2. Learn microenvironments and identify niches

optimize_nicheformer masks each central cell’s molecular state and learns its relationship to the surrounding cells and their relative positions. It returns a niche representation, UMAP coordinates, niche labels, and the score-function weights needed for downstream analysis:

import torch
import mievformer as mf

adata = mf.optimize_nicheformer(
    adata,
    model_path='tutorial_model.pth',
    ngpu=1 if torch.cuda.is_available() else 0,
    max_epochs=1,          # Increase for real analysis
    batch_size=32,
    latent_dim=6,
    neighbor_num=5,       # Use a larger context for real datasets
    ca_reference_num=30,  # Omit for the adaptive default up to 1000
    ca_n_components=5,    # Omit to use automatic dimension selection
    niche_n_neighbors=10,
    random_state=7,
)

The main outputs for exploring tissue organization are:

AnnData key	Meaning
`obsm['e']`	Niche representation derived from reference cell-state profiles
`obsm['X_umap']`	UMAP of the niche representation
`obs['leiden_e']`	Inferred niche labels

sc.pl.umap(adata, color=['leiden_e', 'cell_type'])

The API reference describes the remaining stored outputs and representation parameters.

3. Compare cell-state-dependent niche preferences

The learned density ratio also relates each cell’s molecular state to candidate microenvironments. These soft membership profiles describe which niches are preferred by that state, complementing the niche label of its observed location. The score-function weights are already available after training:

adata = mf.calculate_niche_density_ratio(
    adata, ref_num=30, stratify_key='leiden_e'
)
adata = mf.calculate_niche_cluster_membership(adata)

The relative membership profiles are stored in adata.obsm['dist_e_agg']. They can reveal variation within a population, such as cells with a strong preference for one niche and cells with broader preferences. For example, group the TypeB cells by their profiles:

type_b = adata[adata.obs['cell_type'] == 'TypeB'].copy()
type_b = mf.analyze_niche_membership(
    type_b,
    n_clusters=3,
    file_path='niche_composition_clustermap.png',
)

The current API averages reference scores within each niche cluster, so these profiles are relative values rather than rows of normalized cluster probabilities. See the membership API for the sampling and aggregation details.

4. Map conditional population probabilities

In the other direction, the same learned relationship estimates the probability of a target population given a local tissue context. Here, TypeA is the target:

adata = mf.estimate_population_density(
    adata,
    group='TypeA',
    cluster_key='cell_type',
    max_cell_num=30,
)
sc.pl.embedding(adata, basis='spatial', color='TypeA_density')

adata.obs['TypeA_density'] contains the estimated conditional probability at each cell location. The historical density name refers to a population probability, not a count per unit area. Estimate it on the full tissue before selecting a query cell population, so the reference distribution includes all cell states.

5. Find expression signatures associated with colocalization

To ask which genes vary within TypeB cells as their surroundings become more enriched for TypeA, correlate their expression with the TypeA probability:

query = adata[adata.obs['cell_type'] == 'TypeB'].copy()
sc.pp.normalize_total(query, target_sum=1e4)
sc.pp.log1p(query)
corrs = mf.analyze_density_correlation(
    query,
    density_col='TypeA_density',
    file_path='density_correlation.png',
)
print(corrs.nlargest(5))

Normalization and log transformation are applied here because the synthetic dataset contains counts. For real data, use the expression values appropriate for your analysis without repeating preprocessing already performed. The correlations describe expression associated with TypeA-enriched contexts within TypeB, rather than expression differences between cell types.

For multiple spatial samples, continue with the multi-batch tutorial.