01. Cell embedding with MULTIVI¶

Bolero conditions every prediction on an atlas-scale cell-state embedding, and it is trained on metacells rather than noisy single cells. Two tutorial pages build those inputs for a dataset, using ChromiumPBMC as the worked example:

1. This page — turn a single-cell multiome dataset into a low-dimensional cell embedding with a joint RNA + ATAC MULTIVI model.
2. Meta cells — collapse cells into SEACells metacells on top of that embedding.

Where the data comes from. Datasets are registered in the companion package bolerodata. DATASETS["ChromiumPBMC"] resolves the on-disk single-cell RNA / ATAC AnnData files and the standardized cell metadata from the lab data lake, so you never manage file paths by hand.

Prerequisites. A CUDA GPU, the bolerodata data lake mounted, and the bolero runtime environment (pixi shell). Steps 3–5 are heavy (a full MULTIVI fit over the whole dataset). All intermediate artifacts are written next to this notebook; re-running a step reuses the cached file if it already exists, so a second pass only makes the figures.

Setup¶

In [1]:

from pathlib import Path

import anndata
import joblib
import matplotlib.pyplot as plt
import mudata
import numpy as np
import pandas as pd
import scanpy as sc
import scvi
import seaborn as sns
import torch
from bolerodata import DATASETS
from scipy.sparse import csr_matrix
from scvi.model import MULTIVI
from tqdm.auto import tqdm

scvi.settings.seed = 0
torch.set_float32_matmul_precision("high")
sc.set_figure_params(figsize=(4, 4), frameon=False)

from pathlib import Path import anndata import joblib import matplotlib.pyplot as plt import mudata import numpy as np import pandas as pd import scanpy as sc import scvi import seaborn as sns import torch from bolerodata import DATASETS from scipy.sparse import csr_matrix from scvi.model import MULTIVI from tqdm.auto import tqdm scvi.settings.seed = 0 torch.set_float32_matmul_precision("high") sc.set_figure_params(figsize=(4, 4), frameon=False)

[rank: 0] Seed set to 0

In [2]:

# --- Configuration ---------------------------------------------------------------
# Intermediate artifacts (adata.rna.h5ad, adata.atac.h5ad, model/, *.feather, ...)
# are written to the current directory. Run this notebook from its own folder.
DATASET_NAME = "ChromiumPBMC"

# Cap cells per (subclass x tissue) group so a few abundant states do not dominate
# training of the embedding model.
EXP_BATCH = "donor"  # experimental-batch covariate to regress out
DOWNSAMPLE_GROUPBY = ["subclass", "tissue"]
MAX_CELLS_PER_GROUP = 1500

N_TOP_GENES = 5000  # highly-variable genes kept for the RNA view
N_LATENT = 30  # cell-embedding dimensionality
MAX_EPOCHS = 200  # lower (e.g. 10) for a quick smoke test

dataset = DATASETS[DATASET_NAME]
dataset.name

# --- Configuration --------------------------------------------------------------- # Intermediate artifacts (adata.rna.h5ad, adata.atac.h5ad, model/, *.feather, ...) # are written to the current directory. Run this notebook from its own folder. DATASET_NAME = "ChromiumPBMC" # Cap cells per (subclass x tissue) group so a few abundant states do not dominate # training of the embedding model. EXP_BATCH = "donor" # experimental-batch covariate to regress out DOWNSAMPLE_GROUPBY = ["subclass", "tissue"] MAX_CELLS_PER_GROUP = 1500 N_TOP_GENES = 5000 # highly-variable genes kept for the RNA view N_LATENT = 30 # cell-embedding dimensionality MAX_EPOCHS = 200 # lower (e.g. 10) for a quick smoke test dataset = DATASETS[DATASET_NAME] dataset.name

Out[2]:

'ChromiumPBMC'

Step 1 — Build the training cell list¶

Start from the standardized per-cell metadata and:

• drop extreme-coverage outliers (likely multiplets / debris), then
• downsample each subclass x tissue group to at most MAX_CELLS_PER_GROUP cells.

The downsampled list is used only to train the embedding model efficiently; the trained model is later applied to all cells (Step 5).

In [3]:

cell_metadata = dataset.make_cell_metadata()

# Drop extreme-coverage outliers (2x the 99th percentile of UMIs / fragments).
max_umi = cell_metadata["n_umis"].quantile(0.99) * 2
max_frag = cell_metadata["n_fragments"].quantile(0.99) * 2
keep = (cell_metadata["n_umis"] < max_umi) & (cell_metadata["n_fragments"] < max_frag)
cell_metadata = cell_metadata[keep].copy()
print(f"{keep.sum()} / {keep.size} cells pass the coverage filter")

fig, axes = plt.subplots(figsize=(8, 2), nrows=2, constrained_layout=True)
sns.histplot(cell_metadata["n_umis"], bins=100, ax=axes[0])
sns.histplot(cell_metadata["n_fragments"], bins=100, ax=axes[1])

cell_metadata = dataset.make_cell_metadata() # Drop extreme-coverage outliers (2x the 99th percentile of UMIs / fragments). max_umi = cell_metadata["n_umis"].quantile(0.99) * 2 max_frag = cell_metadata["n_fragments"].quantile(0.99) * 2 keep = (cell_metadata["n_umis"] < max_umi) & (cell_metadata["n_fragments"] < max_frag) cell_metadata = cell_metadata[keep].copy() print(f"{keep.sum()} / {keep.size} cells pass the coverage filter") fig, axes = plt.subplots(figsize=(8, 2), nrows=2, constrained_layout=True) sns.histplot(cell_metadata["n_umis"], bins=100, ax=axes[0]) sns.histplot(cell_metadata["n_fragments"], bins=100, ax=axes[1])

45178 / 45270 cells pass the coverage filter

Out[3]:

<Axes: xlabel='n_fragments', ylabel='Count'>

In [4]:

# One categorical group per (subclass x tissue) combination.
cell_metadata = cell_metadata.dropna(subset=DOWNSAMPLE_GROUPBY)
emb_group = cell_metadata[DOWNSAMPLE_GROUPBY[0]].astype(str)
for col in DOWNSAMPLE_GROUPBY[1:]:
    emb_group += "+" + cell_metadata[col].astype(str)
cell_metadata["emb_group"] = emb_group.astype("category")

use_cells = []
for _, gdf in cell_metadata.groupby("emb_group", observed=True):
    if gdf.shape[0] <= MAX_CELLS_PER_GROUP:
        use_cells.extend(gdf.index)
    else:
        use_cells.extend(gdf.sample(MAX_CELLS_PER_GROUP, random_state=0).index)

train_meta = cell_metadata.loc[use_cells].copy()
train_meta["exp_batch"] = train_meta[EXP_BATCH] if EXP_BATCH else DATASET_NAME
train_meta = train_meta[
    [
        "sample",
        "cluster",
        "subclass",
        "tissue",
        "DissectionRegion",
        "exp_batch",
        "emb_group",
    ]
]
train_meta.to_feather("cell_metadata.feather")
print(f"{cell_metadata.shape[0]} cells -> {train_meta.shape[0]} sampled for training")

# One categorical group per (subclass x tissue) combination. cell_metadata = cell_metadata.dropna(subset=DOWNSAMPLE_GROUPBY) emb_group = cell_metadata[DOWNSAMPLE_GROUPBY[0]].astype(str) for col in DOWNSAMPLE_GROUPBY[1:]: emb_group += "+" + cell_metadata[col].astype(str) cell_metadata["emb_group"] = emb_group.astype("category") use_cells = [] for _, gdf in cell_metadata.groupby("emb_group", observed=True): if gdf.shape[0] <= MAX_CELLS_PER_GROUP: use_cells.extend(gdf.index) else: use_cells.extend(gdf.sample(MAX_CELLS_PER_GROUP, random_state=0).index) train_meta = cell_metadata.loc[use_cells].copy() train_meta["exp_batch"] = train_meta[EXP_BATCH] if EXP_BATCH else DATASET_NAME train_meta = train_meta[ [ "sample", "cluster", "subclass", "tissue", "DissectionRegion", "exp_batch", "emb_group", ] ] train_meta.to_feather("cell_metadata.feather") print(f"{cell_metadata.shape[0]} cells -> {train_meta.shape[0]} sampled for training")

45178 cells -> 15968 sampled for training

Step 2 — RNA view: highly-variable genes¶

Load the gene-count AnnData, restrict it to the training cells, and keep the top N_TOP_GENES highly-variable genes. Raw counts are preserved in X because MULTIVI uses a count likelihood for the RNA modality.

In [5]:

if Path("adata.rna.h5ad").exists():
    rna = anndata.read_h5ad("adata.rna.h5ad")
else:
    rna = scvi.data.read_h5ad(dataset.gene_adata_path)
    rna = rna[rna.obs_names.isin(train_meta.index)].copy()
    sc.pp.filter_genes(rna, min_counts=10)
    sc.pp.calculate_qc_metrics(rna, inplace=True, log1p=True)

    rna.layers["counts"] = rna.X.copy()  # preserve raw counts
    sc.pp.normalize_total(rna, target_sum=1e4)  # HVG selection wants normalized data...
    sc.pp.log1p(rna)
    sc.pp.highly_variable_genes(
        rna, n_top_genes=N_TOP_GENES, subset=True, layer="counts", flavor="seurat_v3"
    )
    rna.X = rna.layers.pop("counts")  # ...but the model wants raw counts back
    rna.write_h5ad("adata.rna.h5ad")

sc.pl.violin(rna, ["n_genes_by_counts", "total_counts"], jitter=0.4, multi_panel=True)
print(rna.shape)

if Path("adata.rna.h5ad").exists(): rna = anndata.read_h5ad("adata.rna.h5ad") else: rna = scvi.data.read_h5ad(dataset.gene_adata_path) rna = rna[rna.obs_names.isin(train_meta.index)].copy() sc.pp.filter_genes(rna, min_counts=10) sc.pp.calculate_qc_metrics(rna, inplace=True, log1p=True) rna.layers["counts"] = rna.X.copy() # preserve raw counts sc.pp.normalize_total(rna, target_sum=1e4) # HVG selection wants normalized data... sc.pp.log1p(rna) sc.pp.highly_variable_genes( rna, n_top_genes=N_TOP_GENES, subset=True, layer="counts", flavor="seurat_v3" ) rna.X = rna.layers.pop("counts") # ...but the model wants raw counts back rna.write_h5ad("adata.rna.h5ad") sc.pl.violin(rna, ["n_genes_by_counts", "total_counts"], jitter=0.4, multi_panel=True) print(rna.shape)

(15968, 5000)

Step 3 — ATAC view: peak counts¶

The peak matrix is stored as per-chunk AnnData files. Keep peaks covered in a reasonable fraction of cells; MULTIVI models the peak counts directly for the accessibility modality.

This writes a large adata.atac.h5ad (~1.6 GB for ChromiumPBMC) next to the notebook.

In [6]:

if Path("adata.atac.h5ad").exists():
    atac = anndata.read_h5ad("adata.atac.h5ad")
else:
    peak_paths = sorted(Path(dataset.peak_adata_path).glob("*.h5ad"))

    # Peak selection: covered above ~1% of cells (and at least the top-100k by signal).
    peak_sum, n_cells = None, 0
    for p in tqdm(peak_paths, desc="scan peaks"):
        a = anndata.read_h5ad(p)
        a = a[a.obs_names.isin(train_meta.index)]
        s = pd.Series(np.asarray(a.X.sum(0)).ravel(), index=a.var_names)
        peak_sum = s if peak_sum is None else peak_sum + s
        n_cells += a.shape[0]
    min_cov = max(
        n_cells * 0.01, peak_sum.sort_values(ascending=False)[:100_000].values[-1]
    )
    use_peaks = peak_sum > min_cov
    print(f"keeping {int(use_peaks.sum())} peaks")

    parts = []
    for p in tqdm(peak_paths, desc="load peaks"):
        a = anndata.read_h5ad(p)[:, use_peaks]
        a = a[a.obs_names.isin(train_meta.index)].copy()
        parts.append(a)
    atac = anndata.concat(parts)
    for col, data in train_meta.items():
        atac.obs[col] = data
    atac.write_h5ad("adata.atac.h5ad")

print(atac.shape)

if Path("adata.atac.h5ad").exists(): atac = anndata.read_h5ad("adata.atac.h5ad") else: peak_paths = sorted(Path(dataset.peak_adata_path).glob("*.h5ad")) # Peak selection: covered above ~1% of cells (and at least the top-100k by signal). peak_sum, n_cells = None, 0 for p in tqdm(peak_paths, desc="scan peaks"): a = anndata.read_h5ad(p) a = a[a.obs_names.isin(train_meta.index)] s = pd.Series(np.asarray(a.X.sum(0)).ravel(), index=a.var_names) peak_sum = s if peak_sum is None else peak_sum + s n_cells += a.shape[0] min_cov = max( n_cells * 0.01, peak_sum.sort_values(ascending=False)[:100_000].values[-1] ) use_peaks = peak_sum > min_cov print(f"keeping {int(use_peaks.sum())} peaks") parts = [] for p in tqdm(peak_paths, desc="load peaks"): a = anndata.read_h5ad(p)[:, use_peaks] a = a[a.obs_names.isin(train_meta.index)].copy() parts.append(a) atac = anndata.concat(parts) for col, data in train_meta.items(): atac.obs[col] = data atac.write_h5ad("adata.atac.h5ad") print(atac.shape)

(15968, 99919)

Step 4 — Train the MULTIVI embedding¶

Recent scvi-tools versions require a MuData object for multiome models: the RNA and ATAC views are kept as two modalities (no var_name collisions, no manual concatenation). setup_mudata's modalities argument maps each model input to the modality that holds it. The exp_batch (donor) covariate is modeled so the latent space reflects cell state rather than experimental batch.

In [7]:

# Align the two views on the same cells and assemble the MuData.
rna = rna[atac.obs_names].copy()
rna.obs["exp_batch"] = atac.obs["exp_batch"].values

rna_var_names = rna.var_names.copy()
atac_var_names = atac.var_names.copy()

mdata = mudata.MuData({"rna": rna, "atac": atac})

MULTIVI.setup_mudata(
    mdata,
    batch_key="exp_batch",
    modalities={"rna_layer": "rna", "atac_layer": "atac", "batch_key": "rna"},
)

# Align the two views on the same cells and assemble the MuData. rna = rna[atac.obs_names].copy() rna.obs["exp_batch"] = atac.obs["exp_batch"].values rna_var_names = rna.var_names.copy() atac_var_names = atac.var_names.copy() mdata = mudata.MuData({"rna": rna, "atac": atac}) MULTIVI.setup_mudata( mdata, batch_key="exp_batch", modalities={"rna_layer": "rna", "atac_layer": "atac", "batch_key": "rna"}, )

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/mudata/_core/mudata.py:1467: FutureWarning: From 0.4 .update() will not pull obs/var columns from individual modalities by default anymore. Set mudata.set_options(pull_on_update=False) to adopt the new behaviour, which will become the default. Use new pull_obs/pull_var and push_obs/push_var methods for more flexibility.
  self._update_attr("var", axis=0, join_common=join_common)
/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/mudata/_core/mudata.py:1318: FutureWarning: From 0.4 .update() will not pull obs/var columns from individual modalities by default anymore. Set mudata.set_options(pull_on_update=False) to adopt the new behaviour, which will become the default. Use new pull_obs/pull_var and push_obs/push_var methods for more flexibility.
  self._update_attr("obs", axis=1, join_common=join_common)

In [8]:

model_dir = "model"
if Path(model_dir).exists():
    model = MULTIVI.load(model_dir, adata=mdata)
else:
    model = MULTIVI(mdata, n_latent=N_LATENT)
    model.train(max_epochs=MAX_EPOCHS, batch_size=512, early_stopping=True)
    model.save(model_dir, overwrite=True)
    joblib.dump(
        {"rna_var_names": rna_var_names, "atac_var_names": atac_var_names},
        f"{model_dir}/features.joblib",
    )
    pd.DataFrame({k: v.squeeze() for k, v in model.history.items()}).to_csv(
        "train_history.csv"
    )

model_dir = "model" if Path(model_dir).exists(): model = MULTIVI.load(model_dir, adata=mdata) else: model = MULTIVI(mdata, n_latent=N_LATENT) model.train(max_epochs=MAX_EPOCHS, batch_size=512, early_stopping=True) model.save(model_dir, overwrite=True) joblib.dump( {"rna_var_names": rna_var_names, "atac_var_names": atac_var_names}, f"{model_dir}/features.joblib", ) pd.DataFrame({k: v.squeeze() for k, v in model.history.items()}).to_csv( "train_history.csv" )

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/torch/nn/init.py:453: UserWarning: Initializing zero-element tensors is a no-op
  warnings.warn("Initializing zero-element tensors is a no-op")
/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/lightning/fabric/plugins/environments/slurm.py:204: The `srun` command is available on your system but is not used. HINT: If your intention is to run Lightning on SLURM, prepend your python command with `srun` like so: srun python /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pi ...
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.
/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/lightning/fabric/plugins/environments/slurm.py:204: The `srun` command is available on your system but is not used. HINT: If your intention is to run Lightning on SLURM, prepend your python command with `srun` like so: srun python /large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pi ...
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/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=47` in the `DataLoader` to improve performance.
/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/.pixi/envs/default/lib/python3.11/site-packages/lightning/pytorch/trainer/connectors/data_connector.py:434: The 'val_dataloader' does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` to `num_workers=47` in the `DataLoader` to improve performance.

Training:   0%|          | 0/200 [00:00<?, ?it/s]

`Trainer.fit` stopped: `max_epochs=200` reached.

In [9]:

# Training curve (available whenever the model was trained in this session or cached to CSV).
if model.history:
    hist = pd.DataFrame({k: v.squeeze() for k, v in model.history.items()})
elif Path("train_history.csv").exists():
    hist = pd.read_csv("train_history.csv", index_col=0)
else:
    hist = None

val_col = None
if hist is not None:
    for c in ["validation_loss", "elbo_validation", "reconstruction_loss_validation"]:
        if c in hist.columns:
            val_col = c
            break

if val_col:
    plt.figure(figsize=(4, 3))
    plt.plot(hist[val_col])
    plt.xlabel("epoch")
    plt.ylabel(val_col)
else:
    print(
        "No training history (model was loaded from cache); train fresh to plot the curve."
    )

# Training curve (available whenever the model was trained in this session or cached to CSV). if model.history: hist = pd.DataFrame({k: v.squeeze() for k, v in model.history.items()}) elif Path("train_history.csv").exists(): hist = pd.read_csv("train_history.csv", index_col=0) else: hist = None val_col = None if hist is not None: for c in ["validation_loss", "elbo_validation", "reconstruction_loss_validation"]: if c in hist.columns: val_col = c break if val_col: plt.figure(figsize=(4, 3)) plt.plot(hist[val_col]) plt.xlabel("epoch") plt.ylabel(val_col) else: print( "No training history (model was loaded from cache); train fresh to plot the curve." )

Step 5 — Apply the model to all cells¶

Training used a downsampled cell list. Now compute the latent embedding for every cell in the dataset, streaming over the peak chunks (one MuData per chunk) to keep memory bounded, and save it as a feather.

In [10]:

emb_path = "multivi.latent_embedding.feather"
if Path(emb_path).exists():
    latent_emb = pd.read_feather(emb_path)
else:
    all_meta = dataset.make_cell_metadata()
    all_meta["exp_batch"] = all_meta[EXP_BATCH] if EXP_BATCH else DATASET_NAME

    rna_all = anndata.read_h5ad(dataset.gene_adata_path, backed="r")
    rna_all = rna_all[:, rna_var_names].to_memory()

    parts = []
    for p in tqdm(sorted(Path(dataset.peak_adata_path).glob("*.h5ad")), desc="infer"):
        a = anndata.read_h5ad(p)
        a = a[a.obs_names.isin(all_meta.index), atac_var_names].copy()
        this_rna = rna_all[a.obs_names].copy()
        this_rna.obs["exp_batch"] = all_meta.loc[a.obs_names, "exp_batch"].values

        mdata_chunk = mudata.MuData({"rna": this_rna, "atac": a})
        MULTIVI.setup_mudata(
            mdata_chunk,
            batch_key="exp_batch",
            modalities={"rna_layer": "rna", "atac_layer": "atac", "batch_key": "rna"},
        )
        z = model.get_latent_representation(adata=mdata_chunk)
        parts.append(pd.DataFrame(z, index=a.obs_names))
    latent_emb = pd.concat(parts)
    latent_emb.to_feather(emb_path)

print(latent_emb.shape)

emb_path = "multivi.latent_embedding.feather" if Path(emb_path).exists(): latent_emb = pd.read_feather(emb_path) else: all_meta = dataset.make_cell_metadata() all_meta["exp_batch"] = all_meta[EXP_BATCH] if EXP_BATCH else DATASET_NAME rna_all = anndata.read_h5ad(dataset.gene_adata_path, backed="r") rna_all = rna_all[:, rna_var_names].to_memory() parts = [] for p in tqdm(sorted(Path(dataset.peak_adata_path).glob("*.h5ad")), desc="infer"): a = anndata.read_h5ad(p) a = a[a.obs_names.isin(all_meta.index), atac_var_names].copy() this_rna = rna_all[a.obs_names].copy() this_rna.obs["exp_batch"] = all_meta.loc[a.obs_names, "exp_batch"].values mdata_chunk = mudata.MuData({"rna": this_rna, "atac": a}) MULTIVI.setup_mudata( mdata_chunk, batch_key="exp_batch", modalities={"rna_layer": "rna", "atac_layer": "atac", "batch_key": "rna"}, ) z = model.get_latent_representation(adata=mdata_chunk) parts.append(pd.DataFrame(z, index=a.obs_names)) latent_emb = pd.concat(parts) latent_emb.to_feather(emb_path) print(latent_emb.shape)

(45270, 30)

Step 6 — Neighbors, UMAP and a first look¶

Assemble an AnnData whose obsm["X_multivi"] holds the embedding, build the kNN graph, and run UMAP + Leiden. This with_coords file is the hand-off to the metacell page.

In [11]:

coords_path = "adata.multivi.with_coords.h5ad"
if Path(coords_path).exists():
    adata = anndata.read_h5ad(coords_path)
else:
    all_meta = dataset.make_cell_metadata()
    cells = all_meta.index.intersection(latent_emb.index)
    adata = anndata.AnnData(
        X=csr_matrix((cells.size, 0)),
        obs=all_meta.reindex(cells),
        obsm={"X_multivi": latent_emb.reindex(cells).values},
    )
    sc.pp.neighbors(adata, use_rep="X_multivi")
    sc.tl.umap(adata, min_dist=0.2)
    sc.tl.leiden(
        adata,
        key_added="multivi_cluster",
        resolution=0.2,
        flavor="igraph",
        n_iterations=2,
    )
    adata.write_h5ad(coords_path)

adata

coords_path = "adata.multivi.with_coords.h5ad" if Path(coords_path).exists(): adata = anndata.read_h5ad(coords_path) else: all_meta = dataset.make_cell_metadata() cells = all_meta.index.intersection(latent_emb.index) adata = anndata.AnnData( X=csr_matrix((cells.size, 0)), obs=all_meta.reindex(cells), obsm={"X_multivi": latent_emb.reindex(cells).values}, ) sc.pp.neighbors(adata, use_rep="X_multivi") sc.tl.umap(adata, min_dist=0.2) sc.tl.leiden( adata, key_added="multivi_cluster", resolution=0.2, flavor="igraph", n_iterations=2, ) adata.write_h5ad(coords_path) adata

Out[11]:

AnnData object with n_obs × n_vars = 45270 × 0
    obs: 'sample', 'cluster', 'n_fragments', 'tsse', 'n_umis', 'n_genes', 'age', 'age_int', 'donor', 'sex', 'tissue', 'DissectionRegion', 'subclass', 'class', 'pmvi_cluster', 'group', 'meta_cell'
    uns: 'neighbors', 'pmvi_cluster', 'umap'
    obsm: 'X_multivi', 'X_umap'
    obsp: 'connectivities', 'distances'

In [15]:

adata.obs["age_int"] = adata.obs["age_int"].astype(float)
sc.pl.umap(
    adata,
    color=["class", "donor"],
    ncols=1,
    wspace=0.3,
    size=8,
)

adata.obs["age_int"] = adata.obs["age_int"].astype(float) sc.pl.umap( adata, color=["class", "donor"], ncols=1, wspace=0.3, size=8, )

---

The embedding is now saved (adata.multivi.with_coords.h5ad and multivi.latent_embedding.feather). Continue to Meta cells to collapse cells into SEACells metacells — the units Bolero actually trains on.