3. Preprocessing#

Pipeline phase: Phase 3

What you will learn:

  • Back up raw counts before normalization

  • Normalize and log-transform

  • Select highly variable genes

  • Run PCA, neighbors, UMAP, and Leiden clustering

  • Use the modality-aware preprocess recipe

Note on integration: The default integration for Visium is scVI. This tutorial uses Harmony (lighter dependency) to stay dependency-free. See the API docs for all sc_tools.pp options.

import numpy as np
import anndata as ad
import sc_tools.pp as pp

Create synthetic data#

np.random.seed(7)
n_spots, n_genes = 400, 300

X = np.random.negative_binomial(4, 0.45, size=(n_spots, n_genes)).astype(float)
gene_names = [f"GENE_{i}" for i in range(n_genes)]

adata = ad.AnnData(
    X=X,
    obs={
        "sample": [f"S{i % 2}" for i in range(n_spots)],
        "library_id": [f"lib_{i % 2}" for i in range(n_spots)],
    },
)
adata.var_names = gene_names
print(adata)

Option A: Step-by-step preprocessing#

This gives you full control over each step.

# 1. Backup raw counts before any transformation
pp.backup_raw(adata)
print("adata.raw backed up:", adata.raw is not None)

# 2. Optionally filter mitochondrial / ribosomal genes before normalization
pp.filter_genes_by_pattern(adata, patterns=["^MT-", "^RPL", "^RPS"])
print(f"After MT/RP filter: {adata.shape}")

# 3. Normalize to 10,000 counts per spot, then log1p transform
pp.normalize_total(adata, target_sum=1e4)
pp.log_transform(adata)

print("X max after log1p:", adata.X.max().round(3))

# 4. PCA on top-N genes (seurat_v3 HVG by default)
pp.pca(adata, n_comps=30, n_top_genes=150)
print("PCA done. obsm keys:", list(adata.obsm.keys()))

# 5. Build neighbor graph and compute UMAP
pp.neighbors(adata, n_neighbors=15)
pp.umap(adata)
print("UMAP done. Shape:", adata.obsm["X_umap"].shape)

# 6. Leiden clustering
pp.leiden(adata, resolution=0.5)
print("Clusters:", adata.obs["leiden"].value_counts().to_dict())

Option B: Recipe (one-call preprocessing)#

The preprocess recipe handles all steps above for a given modality. For IMC data, it uses arcsinh_transform instead of log_transform.

# Re-create fresh adata
adata2 = ad.AnnData(X=X.copy())
adata2.var_names = gene_names
adata2.obs["library_id"] = [f"lib_{i % 2}" for i in range(n_spots)]

pp.preprocess(
    adata2,
    modality="visium",
    integration="harmony",   # scvi is default; harmony used here for lighter deps
    batch_key="library_id",
    n_top_genes=150,
    n_pcs=30,
    resolution=0.5,
)

print("Done. obsm:", list(adata2.obsm.keys()))
print("Clusters:", adata2.obs["leiden"].value_counts().to_dict())

Save Phase 3 checkpoint#

adata.write_h5ad("results/adata.normalized.p3.h5ad")

The checkpoint must have:

  • adata.raw backed up (raw counts)

  • A batch-corrected embedding (e.g. obsm['X_scvi'] or obsm['X_pca_harmony'])

  • obs['leiden'] cluster labels

See Architecture.md §2.2 in the repository for the full requirement table.