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

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_obsm_X_pca"]
print(f"PCA shape: {pca.shape}")

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

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())

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(...)`