Working with cytome datasets

Cytome is a SQLite-based single-file format for single-cell omics data. It stores expression matrices, cell/gene metadata, embeddings, graphs, and ATAC fragment data in a single portable .cytome file.

Key advantages over H5AD:

Streaming access: Read data chunk-by-chunk without loading the full matrix into RAM
SQL metadata queries: Filter cells with SQL WHERE clauses
Single portable file: All modalities (RNA, ATAC, peaks, embeddings) in one file
Compressed storage: Zstd/LZ4 compression on CSR chunks
Provenance tracking: Automatic logging of operations

This tutorial covers:

Creating cytome files from AnnData
Reading and querying data
Streaming access for large datasets
Converting back to AnnData
Using cytome in PIASO pipelines
Developer guide: building methods on cytome

import cytome
import numpy as np
import scipy.sparse as sp
import pandas as pd
import anndata
import piaso
import os
import tempfile
piaso.settings.set_figure_params(style="cell")   # one house style across every figure

1. Creating Cytome Files

1.1 From AnnData

The most common way to create a cytome file is to convert from an existing AnnData object.

# Create a toy AnnData for demonstration
n_cells, n_genes = 5000, 2000
np.random.seed(42)

X = sp.random(n_cells, n_genes, density=0.05, format="csr", dtype=np.float32)
obs = pd.DataFrame({
    "barcode": [f"CELL_{i:04d}" for i in range(n_cells)],
    "cell_type": np.random.choice(["Exc", "Inh", "Astro", "Oligo", "Micro"], n_cells),
    "sample_id": np.random.choice(["S1", "S2", "S3"], n_cells),
    "n_genes": np.array((X > 0).sum(axis=1)).flatten(),
})
var = pd.DataFrame({
    "gene_id": [f"Gene_{i}" for i in range(n_genes)],
    "highly_variable": np.random.choice([True, False], n_genes, p=[0.1, 0.9]),
})

adata = anndata.AnnData(X=X, obs=obs, var=var)
adata.layers["counts"] = X.copy()
adata.obsm["X_pca"] = np.random.randn(n_cells, 50).astype(np.float32)
adata.obsm["X_umap"] = np.random.randn(n_cells, 2).astype(np.float32)

print(adata)

# Convert to cytome
tmpdir = tempfile.mkdtemp(prefix="cytome_tutorial_")
output_path = os.path.join(tmpdir, "tutorial_demo.cytome")
ds = cytome.from_anndata(adata, modality="RNA", output=output_path)
ds.close()  # close so we can reopen cleanly later

print(f"Cytome file created: {output_path}")
print(f"File size: {os.path.getsize(output_path) / 1e6:.1f} MB")

1.2 From H5AD file (backed mode for large files)

For large datasets that don’t fit in RAM, use from_h5ad with backed=True. This streams the expression matrix chunk-by-chunk during conversion.

from cytome.io.convert_anndata import from_h5ad

# Save the toy data as h5ad first
h5ad_path = os.path.join(tmpdir, "tutorial_demo.h5ad")
adata.write_h5ad(h5ad_path)

# Convert from h5ad with backed mode (no full matrix load)
output_backed = os.path.join(tmpdir, "tutorial_demo_backed.cytome")
ds_backed = from_h5ad(
    h5ad_path,
    output=output_backed,
    modality="RNA",
    backed=True,
    chunk_size=1024,      # rows per read chunk
    storage_chunk_size=128  # rows per on-disk blob
)
print(f"Backed conversion complete: {ds_backed.n_cells} cells, {ds_backed.n_genes} genes")
ds_backed.close()

1.3 From scratch

You can also build a cytome file programmatically.

scratch_path = os.path.join(tmpdir, "scratch.cytome")
ds2 = cytome.create(scratch_path)

# Add cell metadata
ds2.set_entity("cells", pd.DataFrame({
    "barcode": ["A", "B", "C", "D", "E"],
    "cell_type": ["Exc", "Inh", "Exc", "Astro", "Oligo"],
}))

# Add gene metadata
ds2.set_entity("genes", pd.DataFrame({
    "gene_id": ["TP53", "EGFR", "MYC"],
}))

# Add an expression matrix
mat = sp.csr_matrix(np.array([
    [10, 0, 5],
    [0, 3, 0],
    [7, 1, 2],
    [0, 0, 8],
    [4, 6, 0],
], dtype=np.float32))
ds2.add_matrix("RNA_counts", mat)

# Add an embedding
ds2.add_embedding("RNA_pca", np.random.randn(5, 3).astype(np.float32))

# Flush to persist all buffered writes atomically
ds2.flush()
print(f"Created: {ds2.n_cells} cells, {ds2.n_genes} genes")
ds2.close()

1.4 From Cell Ranger output

# From standard Cell Ranger
# ds = cytome.from_cellranger("/path/to/cellranger/outs", output="sample.cytome")

# From Cell Ranger ARC (multiome: RNA + ATAC)
# ds = cytome.from_cellranger_arc("/path/to/arc/outs", output="multiome.cytome",
#                                  import_fragments=True, build_index=True)

2. Reading and Querying Data

2.1 Opening a cytome file

ds = cytome.open(output_path)
print(f"Cells:      {ds.n_cells}")
print(f"Genes:      {ds.n_genes}")
print(f"Modalities: {ds.modalities}")

2.2 Cell and gene metadata (EntityTable)

Metadata is stored as SQL tables and accessed through EntityTable objects.

# List columns
print("Cell columns:", ds.cells.columns)
print("Gene columns:", ds.genes.columns)

# Read a column as numpy array
cell_types = ds.cells["cell_type"]
print(f"\nUnique cell types: {np.unique(cell_types)}")
print(f"Cell type counts: {pd.Series(cell_types).value_counts().to_dict()}")

# SQL queries on metadata
exc_cells = ds.cells.query("cell_type = 'Exc'")
print(f"Excitatory cells: {len(exc_cells)}")
print(exc_cells.head())

# Get a boolean mask for subsetting (useful for streaming)
exc_mask = ds.cells.query_mask("cell_type = 'Exc'")
print(f"Mask shape: {exc_mask.shape}, True cells: {exc_mask.sum()}")

# Update metadata
full_counts = ds.RNA.counts.to_memory()
ds.cells["n_genes"] = np.array((full_counts > 0).sum(axis=1)).flatten()
print(f"n_genes range: {ds.cells['n_genes'].min()} - {ds.cells['n_genes'].max()}")

2.3 Expression matrix access (MeasurementLayer)

Expression data is accessed through modality accessors. The matrix is stored as compressed CSR chunks on disk.

# Access the RNA counts layer
counts = ds.RNA.counts
print(f"Shape: {counts.shape}")
print(f"Dtype: {counts.dtype}")

# Slice rows and columns (only reads needed chunks from disk)
subset = counts[0:100, 0:50]  # first 100 cells, first 50 genes
print(f"Subset shape: {subset.shape}, nnz: {subset.nnz}")

# Integer index subsetting (convert boolean mask to indices)
exc_indices = np.where(exc_mask)[0]
exc_counts = counts.rows(exc_indices)
print(f"Exc cells shape: {exc_counts.shape}")

# Read specific rows
rows = counts.rows([0, 10, 100, 1000])
print(f"Selected rows: {rows.shape}")

# Load the full matrix into memory (only for small datasets!)
full_mat = counts.to_memory()
print(f"Full matrix: {full_mat.shape}, {full_mat.nnz} non-zeros")

2.4 Embeddings and graphs

# List available embeddings
print("Cell embeddings:", list(ds.embeddings.keys()))

# Access an embedding
pca = ds.embeddings["RNA_pca"]
print(f"PCA shape: {pca.shape}")

umap = ds.embeddings["RNA_umap"]
print(f"UMAP shape: {umap.shape}")

An obsm array is stored as {modality}_{key} with scanpy’s X_ prefix dropped, so obsm['X_pca'] in the RNA modality becomes RNA_pca. Files written before cytome 0.2.6 use the longer RNA_obsm_X_pca form; both are readable, and to_anndata restores the original obsm key either way. If you are unsure which a file uses, list(ds.embeddings.keys()) above shows it.

2.5 Metadata and provenance

# Metadata store (key-value, JSON-serialized)
print("Metadata keys:", list(ds.metadata.keys())[:10])

# Provenance log
print("\nProvenance log:")
print(ds.provenance.show())

2.x Category order and colours, persisted in the file

set_categories stores a categorical column’s display order and palette inside the cytome, so every plot in every later session honours them without re-declaring anything — one advantage of settings that travel with the data:

ds.set_categories(
    "cell_type",
    order=["Astro", "Oligo", "OPC", "Micro", "Endo"],
    colors={"Astro": "#4c72b0", "Oligo": "#dd8452", "OPC": "#55a868",
            "Micro": "#c44e52", "Endo": "#8172b3"},
)
# piaso.pl.plotEmbedding(ds, color="cell_type") now uses this order + palette,
# including after ds.close() / cytome.open() in a fresh session.

(AnnData gets the same effect per session via the categorical dtype’s order plus adata.uns["cell_type_colors"] — see the plotting tutorial’s “Category order and colours” section.)

3. Streaming Access (Large Datasets)

The key advantage of cytome over h5ad is streaming: you can process millions of cells with bounded RAM by iterating over chunks.

3.1 Basic chunk iteration

# Iterate over chunks (each chunk is a (csr_matrix, row_indices) tuple)
total_nnz = 0
n_chunks = 0

for chunk, row_idx in ds.iter_chunks(modality="RNA", layer="counts"):
    total_nnz += chunk.nnz
    n_chunks += 1

print(f"Chunks: {n_chunks}")
print(f"Total non-zeros: {total_nnz}")
print(f"Matches full matrix: {total_nnz == full_mat.nnz}")

3.2 Streaming with cell filtering

Pass a cell_mask to skip chunks that don’t contain cells of interest. Chunks with no matching cells are skipped entirely (no disk I/O).

# Stream only excitatory cells
exc_mask = ds.cells.query_mask("cell_type = 'Exc'")
exc_nnz = 0
for chunk, row_idx in ds.iter_chunks(modality="RNA", layer="counts", cell_mask=exc_mask):
    exc_nnz += chunk.nnz

print(f"Exc cells non-zeros: {exc_nnz}")

3.3 Streaming with batch size control

Use batch_size to control the trade-off between RAM usage and compute efficiency. Larger batches reduce Python overhead but use more memory.

# Larger batches for compute efficiency
n_chunks_large = 0
for chunk, row_idx in ds.iter_chunks(modality="RNA", layer="counts", batch_size=1024):
    n_chunks_large += 1

print(f"With batch_size=1024: {n_chunks_large} batches (vs {n_chunks} raw chunks)")

3.4 Example: streaming mean expression per gene

Here’s a practical example of computing gene means without loading the full matrix.

# Streaming mean computation
gene_sums = np.zeros(ds.n_genes, dtype=np.float64)
total_cells = 0

for chunk, row_idx in ds.iter_chunks(modality="RNA", layer="counts", batch_size=1024):
    gene_sums += np.asarray(chunk.sum(axis=0)).flatten()
    total_cells += chunk.shape[0]

gene_means = gene_sums / total_cells

# Verify against full matrix
expected = np.asarray(full_mat.mean(axis=0)).flatten()
print(f"Max difference: {np.max(np.abs(gene_means - expected)):.2e}")
print(f"Top 5 expressed genes: {np.argsort(gene_means)[-5:][::-1]}")

3.5 Row-level and column-level iteration

The MeasurementLayer also supports low-level chunk iteration.

# Row-level iteration (CSR chunks)
for row_start, row_end, csr_chunk in ds.RNA.counts.iter_rows():
    # Process each on-disk chunk
    pass

print(f"Last chunk: rows {row_start}-{row_end}, shape {csr_chunk.shape}")

4. Converting Back to AnnData

Cytome supports full round-trip conversion to AnnData. All layers, embeddings, metadata, and graphs are preserved.

from cytome.io.convert_anndata import to_anndata

# Full round-trip
adata_rt = to_anndata(ds, modality="RNA")
print(adata_rt)
print(f"\nobs columns: {list(adata_rt.obs.columns)}")
print(f"var columns: {list(adata_rt.var.columns)}")
print(f"obsm keys:   {list(adata_rt.obsm.keys())}")
print(f"layers:      {list(adata_rt.layers.keys())}")

# Partial export with cell mask (e.g., only excitatory cells)
exc_mask = ds.cells.query_mask("cell_type = 'Exc'")
adata_exc = to_anndata(ds, modality="RNA", cell_mask=exc_mask)
print(f"Excitatory subset: {adata_exc.shape}")

5. Using Cytome in PIASO Pipelines

PIASO’s streaming functions — score, runGDR, COSG — work directly on cytome files without loading the full expression matrix.

5.1 Streaming score normalization

# score() accepts a cytome path directly
# piaso.tl.calculateScoreParallel(
#     adata,
#     cytome_path="/path/to/dataset.cytome",
#     modality="RNA",
#     batch_size=1024,
# )
#
# This streams through the expression matrix in chunks:
#   Pass 1: compute per-gene mean/variance (for KNN parameters)
#   Pass 2: apply fused normalization and accumulate scores
#
# Peak RAM stays at O(batch_size x n_genes) regardless of dataset size.
print("See PIASO_tutorial.ipynb for full score/GDR pipeline examples.")

5.2 Streaming COSG marker genes

# COSG on cytome: single-pass streaming, bounded RAM
# piaso.tl.runGDRParallel(
#     adata,
#     cytome_path="/path/to/dataset.cytome",
#     modality="RNA",
#     batch_size=2048,
# )
#
# Internally calls run_cosg_cytome_cpu() which accumulates
# dot-product scores in O(n_genes x n_groups) float32 buffers.
print("See cosg_cytome_benchmark.ipynb for detailed COSG benchmarks.")

6. Developer Guide: Building Methods on Cytome

This section is for developers who want to build new streaming algorithms on top of cytome.

6.1 Pattern: streaming accumulation

The core pattern for cytome-aware algorithms is streaming accumulation: iterate over chunks, update fixed-size accumulators, and finalize after the last chunk.

def streaming_variance(ds, modality="RNA", layer="counts", batch_size=1024):
    """Compute per-gene variance in a single streaming pass (Welford's method).

    Peak RAM: O(n_genes) regardless of n_cells.
    """
    n_genes = ds.n_genes
    col_sum = np.zeros(n_genes, dtype=np.float64)
    col_sq_sum = np.zeros(n_genes, dtype=np.float64)
    total_cells = 0

    for chunk, row_idx in ds.iter_chunks(modality=modality, layer=layer, batch_size=batch_size):
        col_sum += np.asarray(chunk.sum(axis=0)).flatten()
        # For squared sum, operate on .data to avoid densifying
        chunk_sq = chunk.copy()
        chunk_sq.data **= 2
        col_sq_sum += np.asarray(chunk_sq.sum(axis=0)).flatten()
        total_cells += chunk.shape[0]

    mean = col_sum / total_cells
    variance = col_sq_sum / total_cells - mean ** 2
    return mean, variance


# Test it
mean_stream, var_stream = streaming_variance(ds)

# Verify against full matrix
dense = full_mat.toarray()
print(f"Mean max error:     {np.max(np.abs(mean_stream - dense.mean(axis=0))):.2e}")
print(f"Variance max error: {np.max(np.abs(var_stream - dense.var(axis=0))):.2e}")

6.2 Pattern: streaming with group labels

Many algorithms need per-group statistics (e.g., marker genes per cluster). Read labels from the SQL metadata and accumulate per group.

def streaming_group_means(ds, groupby, modality="RNA", layer="counts", batch_size=1024):
    """Compute per-group mean expression in a single streaming pass.

    Peak RAM: O(n_groups x n_genes).
    """
    # Read labels from SQL metadata
    labels = ds.cells[groupby]
    unique_labels = np.unique(labels)
    label_to_idx = {l: i for i, l in enumerate(unique_labels)}
    n_groups = len(unique_labels)
    n_genes = ds.n_genes

    group_sums = np.zeros((n_groups, n_genes), dtype=np.float64)
    group_counts = np.zeros(n_groups, dtype=np.int64)

    for chunk, row_idx in ds.iter_chunks(modality=modality, layer=layer, batch_size=batch_size):
        chunk_labels = labels[row_idx]

        for label in unique_labels:
            mask = chunk_labels == label
            if mask.any():
                idx = label_to_idx[label]
                group_sums[idx] += np.asarray(chunk[mask].sum(axis=0)).flatten()
                group_counts[idx] += mask.sum()

    group_means = group_sums / group_counts[:, None]
    return pd.DataFrame(group_means, index=unique_labels,
                        columns=ds.genes["gene_id"])


means_df = streaming_group_means(ds, groupby="cell_type")
print(f"Group means shape: {means_df.shape}")
print(means_df.iloc[:, :5])

6.3 Pattern: streaming matrix write

Use create_layer_writer() to write large matrices chunk-by-chunk without materializing the full result.

# Example: streaming log1p normalization
writer = ds.create_layer_writer(
    layer_name="RNA_log1p",
    n_rows=ds.n_cells,
    n_cols=ds.n_genes,
    dtype=np.float32,
    compression="zstd",
    storage_chunk_size=128,
)

for chunk, row_idx in ds.iter_chunks(modality="RNA", layer="counts", batch_size=1024):
    # Normalize: log1p of library-size-normalized counts
    lib_size = np.asarray(chunk.sum(axis=1)).flatten()
    lib_size[lib_size == 0] = 1  # avoid division by zero
    median_lib = np.median(lib_size[lib_size > 0])

    # Sparse-friendly: scale in-place on .data, then log1p
    normalized = chunk.astype(np.float64).multiply(median_lib / lib_size[:, None])
    normalized = normalized.tocsr()
    np.log1p(normalized.data, out=normalized.data)
    normalized = normalized.astype(np.float32)

    writer.write_chunk(normalized, row_idx[0])

writer.finalize()
print(f"Wrote RNA_log1p layer: {ds.n_cells} x {ds.n_genes}")

6.4 Key design principles for cytome-aware methods

Never call to_memory() on large datasets. Use iter_chunks() or iter_rows().
Pre-allocate accumulators as float64 arrays sized O(n_genes) or O(n_groups x n_genes).
Read labels once from ds.cells[groupby] before the streaming loop.
Use cell_mask to skip irrelevant chunks when only a subset of cells is needed.
Flush after writes: call ds.flush() or writer.finalize() to persist changes.
Log provenance with ds.provenance.log() after completing an operation.
Sparse operations: use .data, .indices, .indptr for zero-copy CSR manipulation.

6.5 SQLite internals

Cytome stores everything in a single SQLite database. Developers can access it directly for advanced use cases.

# Use the dataset's own connection — no need to open a separate sqlite3 handle.
# `ds._conn` is the raw SQLite connection if you want low-level access.
ds = cytome.open(output_path)

# Print row counts via the public API
print("Entity tables:")
print(f"  cells: {ds.n_cells:,} rows")
print(f"  genes: {ds.n_genes:,} rows")
print(f"  peaks: {ds.n_peaks:,} rows")

# Embeddings & matrices
print(f"\nMatrices: {ds.list_matrices()}")
print(f"Embeddings: {ds.list_embeddings()}")

# For advanced introspection, the raw SQLite connection is still available.
# This lists every table in the underlying database.
tables = ds._conn.execute(
    "SELECT name FROM sqlite_master WHERE type='table' ORDER BY name"
).fetchall()
print(f"\nAll SQLite tables: {len(tables)}")
for (t,) in tables[:10]:
    count = ds._conn.execute(f'SELECT COUNT(*) FROM "{t}"').fetchone()[0]
    print(f"  {t}: {count} rows")
if len(tables) > 10:
    print(f"  ... and {len(tables) - 10} more (mostly per-chromosome fragment tables)")

ds.close()

7. Cleanup

ds.close()

# Clean up temp directory
import shutil
if os.path.exists(tmpdir):
    shutil.rmtree(tmpdir)
    print(f"Removed {tmpdir}")

Summary

Operation	Function / Method
Create from AnnData	`cytome.from_anndata(adata, modality, output)`
Create from H5AD (backed)	`from_h5ad(path, output, backed=True)` from `cytome.io.convert_anndata`
Create from Cell Ranger	`cytome.from_cellranger(path, output)`
Create empty	`cytome.create(path)`
Open existing	`cytome.open(path)`
Merge datasets	`cytome.merge(inputs, output)`
Cell metadata	`ds.cells[col]`, `ds.cells.query(sql)`, `ds.cells.query_mask(sql)`
Gene metadata	`ds.genes[col]`
Expression matrix	`ds.RNA.counts[rows, cols]`, `.to_memory()`, `.rows(idx)`
Stream chunks	`ds.iter_chunks(modality, layer, cell_mask, batch_size)`
Row iteration	`ds.RNA.counts.iter_rows()`
Embeddings	`ds.embeddings[name]`
Write matrix	`ds.add_matrix(name, sparse)` then `ds.flush()`
Stream write	`ds.create_layer_writer(...)` then `.write_chunk()` then `.finalize()`
Convert to AnnData	`to_anndata(ds, modality)` from `cytome.io.convert_anndata`
Provenance	`ds.provenance.log(...)`