04. Parquet dataset¶

The metacell AnnDatas from the previous page summarize accessibility per peak. To train on DNA sequence, Bolero instead needs accessibility at base resolution across the whole genome. This page builds that store — a parquet dataset (a "parquet_db") — for ChromiumPBMC, then shows what it contains and how to read it.

The build has two stages: sum each metacell's SnapATAC2 fragment insertions into a per-chromosome sparse matrix, then tile the genome into fixed windows and write each window as one compressed row of a parquet file. We build the 32 bp resolution used in the paper; the same builder produces other resolutions by changing three numbers.

Where the data comes from. DATASETS["ChromiumPBMC"] resolves the per-sample SnapATAC2 fragment files (.snap_file_table) and the genome (hg38). The cell → metacell assignment is read from adata.multivi.with_coords.h5ad produced two pages back.

This is a heavy step. It reads genome-wide fragments for every cell and writes a multi-GB dataset, so it wants many CPUs. Run it in place from its own folder.

Setup¶

In [1]:

Copied!





from pathlib import Path

import anndata
import joblib
import numpy as np
import pandas as pd
import pyarrow.parquet as pq
from scipy.sparse import csr_matrix

from bolero import init
from bolero.pp.ray_chunk_dataset import GenomeChunkDatasetGenerator
from bolero.pp.snap_adata import AdataPseudobulkMerger
from bolero.tl.dataset.parquet_db import GenomeParquetDB
from bolero.tl.dataset.sc_transforms import compressed_bytes_to_array
from bolerodata import DATASETS

init(num_cpus=32)
from pathlib import Path

import anndata
import joblib
import numpy as np
import pandas as pd
import pyarrow.parquet as pq
from scipy.sparse import csr_matrix

from bolero import init
from bolero.pp.ray_chunk_dataset import GenomeChunkDatasetGenerator
from bolero.pp.snap_adata import AdataPseudobulkMerger
from bolero.tl.dataset.parquet_db import GenomeParquetDB
from bolero.tl.dataset.sc_transforms import compressed_bytes_to_array
from bolerodata import DATASETS

init(num_cpus=32)

2026-07-12 20:08:06,232	INFO worker.py:1781 -- Started a local Ray instance.

In [2]:

Copied!





# --- Configuration ---------------------------------------------------------------
DATASET_NAME = "ChromiumPBMC"
INPUT_ADATA = "adata.multivi.with_coords.h5ad"  # per-cell, carries obs["meta_cell"]

# 32 bp coverage -- the resolution used for the paper.
WINDOW_SIZE = 600_000  # bp spanned by each stored meta-region (one parquet row)
STEP_SIZE = 200_000  # stride between meta-regions (windows overlap)
RESOLUTION = 32  # bp per bin -> WINDOW_SIZE // RESOLUTION bins per row

CSC_DIR = f"{DATASET_NAME}_csc_dataset"  # intermediate per-metacell CSC store
PARQUET_DIR = f"{DATASET_NAME}-MetaCell-{RESOLUTION}bp"  # the final parquet_db
PREFIX = f"{DATASET_NAME}.MetaCell"  # parquet column-name prefix

dataset = DATASETS[DATASET_NAME]
# --- Configuration ---------------------------------------------------------------
DATASET_NAME = "ChromiumPBMC"
INPUT_ADATA = "adata.multivi.with_coords.h5ad"  # per-cell, carries obs["meta_cell"]

# 32 bp coverage -- the resolution used for the paper.
WINDOW_SIZE = 600_000  # bp spanned by each stored meta-region (one parquet row)
STEP_SIZE = 200_000  # stride between meta-regions (windows overlap)
RESOLUTION = 32  # bp per bin -> WINDOW_SIZE // RESOLUTION bins per row

CSC_DIR = f"{DATASET_NAME}_csc_dataset"  # intermediate per-metacell CSC store
PARQUET_DIR = f"{DATASET_NAME}-MetaCell-{RESOLUTION}bp"  # the final parquet_db
PREFIX = f"{DATASET_NAME}.MetaCell"  # parquet column-name prefix

dataset = DATASETS[DATASET_NAME]

Step 1 — Cell → metacell map¶

AdataPseudobulkMerger needs a table mapping each cell (by sample and barcode) to the metacell it belongs to. We read the per-cell metadata, keep sample + meta_cell, and strip the {sample}+ prefix off the barcode so it matches the cell names in the SnapATAC2 files.

In [3]:

Copied!





adata = anndata.read_h5ad(INPUT_ADATA)
cell_metadata = adata.obs
n_meta_cells = cell_metadata["meta_cell"].nunique()

pseudobulk_meta = (
    cell_metadata[["sample", "meta_cell"]].reset_index().iloc[:, [1, 2, 0]]
)
pseudobulk_meta.columns = ["sample", "group", "bc"]
pseudobulk_meta["bc"] = pseudobulk_meta["bc"].map(lambda name: name.split("+")[-1])

snap_file_dict = (
    dataset.snap_file_table.set_index("sample")["snap_path"].dropna().to_dict()
)
snap_file_dict = {
    s: p for s, p in snap_file_dict.items() if s in set(pseudobulk_meta["sample"])
}
print(f"{n_meta_cells} metacells across {len(snap_file_dict)} samples")
pseudobulk_meta.head()
adata = anndata.read_h5ad(INPUT_ADATA)
cell_metadata = adata.obs
n_meta_cells = cell_metadata["meta_cell"].nunique()

pseudobulk_meta = (
    cell_metadata[["sample", "meta_cell"]].reset_index().iloc[:, [1, 2, 0]]
)
pseudobulk_meta.columns = ["sample", "group", "bc"]
pseudobulk_meta["bc"] = pseudobulk_meta["bc"].map(lambda name: name.split("+")[-1])

snap_file_dict = (
    dataset.snap_file_table.set_index("sample")["snap_path"].dropna().to_dict()
)
snap_file_dict = {
    s: p for s, p in snap_file_dict.items() if s in set(pseudobulk_meta["sample"])
}
print(f"{n_meta_cells} metacells across {len(snap_file_dict)} samples")
pseudobulk_meta.head()

619 metacells across 4 samples

Out[3]:

	sample	group	bc
0	pbmc_sorted_10k	group5.1+SEACell-34	AAACAGCCAAGGAATC-1
1	pbmc_sorted_10k	group6+SEACell-1	AAACAGCCAATCCCTT-1
2	pbmc_sorted_10k	group5.1+SEACell-34	AAACAGCCAATGCGCT-1
3	pbmc_sorted_10k	group5.1+SEACell-27	AAACAGCCACACTAAT-1
4	pbmc_sorted_10k	group3+SEACell-28	AAACAGCCACCAACCG-1

Step 2 — Sum fragments into a metacell CSC store¶

AdataPseudobulkMerger streams the SnapATAC2 fragment matrices and sums the cells of each metacell, writing one sparse metacell × base matrix per chromosome (insertion_<chrom>.csc.joblib) under an intermediate directory. This is the expensive, genome-wide aggregation.

In [4]:

Copied!





merger = AdataPseudobulkMerger(
    pseudobulk_meta=pseudobulk_meta,
    adata_path_dict=snap_file_dict,
    output_dir=CSC_DIR,
    sparse_format="csc",
)
merger.merge()
del merger
merger = AdataPseudobulkMerger(
    pseudobulk_meta=pseudobulk_meta,
    adata_path_dict=snap_file_dict,
    output_dir=CSC_DIR,
    sparse_format="csc",
)
merger.merge()
del merger

pseudobulk meta table: 4 samples, 619 groups, 45270 cells

dump insertion_chr3 groups merge

dump insertion_chr14 groups merge

dump insertion_chr15 groups merge

dump insertion_chr2 groups merge

dump insertion_chr18 groups merge

dump insertion_chr21 groups merge

dump insertion_chr22 groups merge

dump insertion_chr4 groups merge

dump insertion_chr16 groups merge

dump insertion_chr9 groups merge

dump insertion_chr10 groups merge

dump insertion_chrY groups merge

dump insertion_chr6 groups merge

dump insertion_chr5 groups merge

dump insertion_chr8 groups merge

dump insertion_chr20 groups merge

dump insertion_chr11 groups merge

dump insertion_chr13 groups merge

dump insertion_chrX groups merge

dump insertion_chr7 groups merge

dump insertion_chr1 groups merge

dump insertion_chr19 groups merge

dump insertion_chr17 groups merge

dump insertion_chr12 groups merge

Step 3 — Build the 32 bp parquet dataset¶

GenomeChunkDatasetGenerator tiles the genome into WINDOW_SIZE-bp windows (stepping every STEP_SIZE bp) and, for each window, sums the base-level coverage into RESOLUTION-bp bins and writes it as one compressed parquet row. add_chrom_sparse(..., resolution=32) is the knob that sets the bin size.

In [5]:

Copied!





num_rows_per_file = max(10_000 // n_meta_cells, 5)

generator = GenomeChunkDatasetGenerator(
    output_dir=f"{PARQUET_DIR}/",
    genome=dataset.genome,
    window_size=WINDOW_SIZE,
    step_size=STEP_SIZE,
    num_rows_per_file=num_rows_per_file,
)
generator.add_chrom_sparse(prefix=PREFIX, dataset_dir=CSC_DIR, resolution=RESOLUTION)
generator.generate()
del generator
num_rows_per_file = max(10_000 // n_meta_cells, 5)

generator = GenomeChunkDatasetGenerator(
    output_dir=f"{PARQUET_DIR}/",
    genome=dataset.genome,
    window_size=WINDOW_SIZE,
    step_size=STEP_SIZE,
    num_rows_per_file=num_rows_per_file,
)
generator.add_chrom_sparse(prefix=PREFIX, dataset_dir=CSC_DIR, resolution=RESOLUTION)
generator.generate()
del generator

Processing ChromiumPBMC.MetaCell chr1

Processing ChromiumPBMC.MetaCell chr10

Processing ChromiumPBMC.MetaCell chr11

Processing ChromiumPBMC.MetaCell chr12

Processing ChromiumPBMC.MetaCell chr13

Processing ChromiumPBMC.MetaCell chr14

Processing ChromiumPBMC.MetaCell chr15

Processing ChromiumPBMC.MetaCell chr16

Processing ChromiumPBMC.MetaCell chr17

Processing ChromiumPBMC.MetaCell chr18

Processing ChromiumPBMC.MetaCell chr19

Processing ChromiumPBMC.MetaCell chr2

Processing ChromiumPBMC.MetaCell chr20

Processing ChromiumPBMC.MetaCell chr21

Processing ChromiumPBMC.MetaCell chr22

Processing ChromiumPBMC.MetaCell chr3

Processing ChromiumPBMC.MetaCell chr4

Processing ChromiumPBMC.MetaCell chr5

Processing ChromiumPBMC.MetaCell chr6

Processing ChromiumPBMC.MetaCell chr7

Processing ChromiumPBMC.MetaCell chr8

Processing ChromiumPBMC.MetaCell chr9

Processing ChromiumPBMC.MetaCell chrX

Processing ChromiumPBMC.MetaCell chrY

Creating dataset for chromosome chr1.

2026-07-12 20:18:06,715	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:06,716	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr10.

2026-07-12 20:18:07,791	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:07,791	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr11.

2026-07-12 20:18:08,425	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:08,426	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr12.

2026-07-12 20:18:09,101	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:09,102	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr13.

2026-07-12 20:18:09,571	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:09,572	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr14.

2026-07-12 20:18:10,059	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:10,060	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr15.

2026-07-12 20:18:10,569	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:10,570	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr16.

2026-07-12 20:18:11,115	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:11,115	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr17.

2026-07-12 20:18:11,792	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:11,793	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr18.

2026-07-12 20:18:12,300	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:12,300	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr19.

2026-07-12 20:18:12,869	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:12,869	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr2.

2026-07-12 20:18:13,662	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:13,662	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr20.

2026-07-12 20:18:14,250	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:14,250	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr21.

2026-07-12 20:18:14,659	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:14,659	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr22.

2026-07-12 20:18:15,108	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:15,108	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr3.

2026-07-12 20:18:15,767	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:15,768	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr4.

2026-07-12 20:18:16,408	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:16,409	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr5.

2026-07-12 20:18:16,996	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:16,997	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr6.

2026-07-12 20:18:17,640	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:17,641	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr7.

2026-07-12 20:18:18,306	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:18,307	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr8.

2026-07-12 20:18:18,890	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:18,891	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chr9.

2026-07-12 20:18:19,471	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:19,471	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chrX.

2026-07-12 20:18:19,967	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:19,967	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

2026-07-12 20:18:20,285	INFO streaming_executor.py:108 -- Starting execution of Dataset. Full logs are in /tmp/ray/session_2026-07-12_20-08-03_897060_1845820/logs/ray-data

2026-07-12 20:18:20,285	INFO streaming_executor.py:109 -- Execution plan of Dataset: InputDataBuffer[Input] -> TaskPoolMapOperator[Write]

Creating dataset for chromosome chrY.

Other resolutions. The same builder produces any resolution. For the 1 bp store used elsewhere, set WINDOW_SIZE = 100_000, STEP_SIZE = 90_000, RESOLUTION = 1 and rerun Step 3 into a -MetaCell-1bp directory. We build only 32 bp here because that is what the paper's model uses.

What's inside the parquet dataset¶

A parquet_db is a directory of parquet files, partitioned by chromosome, plus three small sidecar files. Each row of a parquet file is one genomic window ("meta-region") holding a gzip-compressed CSR matrix of shape (metacells, bins). The matrix is split across four binary columns following the grammar "{prefix}:{part}+{dtype}" — the CSR indices, indptr, data, and shape — alongside a region string like chr1:0-600000.

The sidecars are: config.joblib (the build parameters), row_names.joblib (the metacell order of every matrix's rows), and parquet_row_regions.feather (the region → file lookup used for random access).

In [6]:

Copied!





root = Path(PARQUET_DIR)
print("top level:", sorted(p.name for p in root.iterdir())[:12])

# the metacell order shared by every stored matrix's rows
row_names = joblib.load(root / "row_names.joblib")
metacells = row_names[PREFIX]
print(f"{len(metacells)} metacell rows, e.g. {list(metacells[:3])}")

# schema of one parquet file
example_file = next((root / "chr1").glob("*.parquet"))
print(pq.read_schema(example_file))
root = Path(PARQUET_DIR)
print("top level:", sorted(p.name for p in root.iterdir())[:12])

# the metacell order shared by every stored matrix's rows
row_names = joblib.load(root / "row_names.joblib")
metacells = row_names[PREFIX]
print(f"{len(metacells)} metacell rows, e.g. {list(metacells[:3])}")

# schema of one parquet file
example_file = next((root / "chr1").glob("*.parquet"))
print(pq.read_schema(example_file))

top level: ['chr1', 'chr10', 'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19', 'chr2']
619 metacell rows, e.g. ['group0+SEACell-0', 'group0+SEACell-1', 'group0+SEACell-2']
ChromiumPBMC.MetaCell:indices+uint32: binary
ChromiumPBMC.MetaCell:indptr+uint32: binary
ChromiumPBMC.MetaCell:data+float32: binary
ChromiumPBMC.MetaCell:shape+uint32: binary
region: string

In [7]:

Copied!

# region -> file lookup that powers random access
regions = pd.read_feather(root / "parquet_row_regions.feather")
regions.head()
# region -> file lookup that powers random access
regions = pd.read_feather(root / "parquet_row_regions.feather")
regions.head()

Out[7]:

	region	parquet_file	row_id_in_parquet
0	chr4:108000000-108600000	chr4/64_000036_000000.parquet	0
1	chr4:108200000-108800000	chr4/64_000036_000000.parquet	1
2	chr4:108400000-109000000	chr4/64_000036_000000.parquet	2
3	chr4:108600000-109200000	chr4/64_000036_000000.parquet	3
4	chr4:108800000-109400000	chr4/64_000036_000000.parquet	4

Demo 1 — Decode one meta-region into a matrix¶

At the lowest level, reading the dataset means decompressing the four CSR arrays of a row and rebuilding the sparse matrix. bolero.tl.dataset.sc_transforms.compressed_bytes_to_array is the one-line helper (np.frombuffer(gzip.decompress(b), dtype)); the data loader uses exactly this.

In [8]:

Copied!





table = pq.read_table(example_file)

# pick the first meta-region that actually has coverage (telomeric windows are empty)
row_i = next(
    i
    for i in range(table.num_rows)
    if compressed_bytes_to_array(
        table[f"{PREFIX}:data+float32"][i].as_py(), "float32"
    ).size
)
row = table.slice(row_i, 1).to_pydict()

indices = compressed_bytes_to_array(row[f"{PREFIX}:indices+uint32"][0], "uint32")
indptr = compressed_bytes_to_array(row[f"{PREFIX}:indptr+uint32"][0], "uint32")
data = compressed_bytes_to_array(row[f"{PREFIX}:data+float32"][0], "float32")
shape = compressed_bytes_to_array(row[f"{PREFIX}:shape+uint32"][0], "uint32")

mat = csr_matrix((data, indices, indptr), shape=tuple(shape))
print(f"{row['region'][0]} -> {mat.shape} {mat.dtype}, {mat.nnz:,} nonzero entries")
# rows = metacells (row_names order); columns = 32 bp bins (18750 = 600000 / 32)
table = pq.read_table(example_file)

# pick the first meta-region that actually has coverage (telomeric windows are empty)
row_i = next(
    i
    for i in range(table.num_rows)
    if compressed_bytes_to_array(
        table[f"{PREFIX}:data+float32"][i].as_py(), "float32"
    ).size
)
row = table.slice(row_i, 1).to_pydict()

indices = compressed_bytes_to_array(row[f"{PREFIX}:indices+uint32"][0], "uint32")
indptr = compressed_bytes_to_array(row[f"{PREFIX}:indptr+uint32"][0], "uint32")
data = compressed_bytes_to_array(row[f"{PREFIX}:data+float32"][0], "float32")
shape = compressed_bytes_to_array(row[f"{PREFIX}:shape+uint32"][0], "uint32")

mat = csr_matrix((data, indices, indptr), shape=tuple(shape))
print(f"{row['region'][0]} -> {mat.shape} {mat.dtype}, {mat.nnz:,} nonzero entries")
# rows = metacells (row_names order); columns = 32 bp bins (18750 = 600000 / 32)

chr1:97200000-97800000 -> (619, 18750) float32, 76,355 nonzero entries

Demo 2 — Random access with `GenomeParquetDB`¶

You rarely decode rows by hand. GenomeParquetDB wraps the parquet_db with a DuckDB region index: give it any genomic region and it finds the right file(s), decompresses, and slices out the exact window at the dataset's resolution. iter_batches yields dense metacell × bin arrays — the shape Bolero's data loader (bolero.tl.dataset.ray_dataset.RayGenomeChunkDataset) feeds to the model.

In [9]:

Copied!





db = GenomeParquetDB(dataset_dir=PARQUET_DIR)

# a 100 kb window inside a region we know is covered (reuse the one found above)
chrom, span = row["region"][0].split(":")
start = int(span.split("-")[0]) + 100_000
region = f"{chrom}:{start}-{start + 100_000}"

result = db.query_regions([region])[0]
coverage = result[PREFIX]  # CSR: metacells x bins for this window
print(
    region,
    "->",
    coverage.shape,
    f"({100_000} bp / {RESOLUTION} = {100_000 // RESOLUTION} bins)",
)

batch = next(db.iter_batches([region], batch_size=1))
print({k: np.asarray(v).shape for k, v in batch.items() if k != "region"})
db = GenomeParquetDB(dataset_dir=PARQUET_DIR)

# a 100 kb window inside a region we know is covered (reuse the one found above)
chrom, span = row["region"][0].split(":")
start = int(span.split("-")[0]) + 100_000
region = f"{chrom}:{start}-{start + 100_000}"

result = db.query_regions([region])[0]
coverage = result[PREFIX]  # CSR: metacells x bins for this window
print(
    region,
    "->",
    coverage.shape,
    f"({100_000} bp / {RESOLUTION} = {100_000 // RESOLUTION} bins)",
)

batch = next(db.iter_batches([region], batch_size=1))
print({k: np.asarray(v).shape for k, v in batch.items() if k != "region"})

chr1:97300000-97400000 -> (619, 3125) (100000 bp / 32 = 3125 bins)
{'ChromiumPBMC.MetaCell': (1, 619, 3125)}

Housekeeping¶

ChromiumPBMC_csc_dataset/ is only an intermediate — the parquet dataset is self-contained, so the CSC store can be deleted once you're satisfied with the build. This notebook leaves it in place; remove it manually when you no longer need it, e.g. shutil.rmtree(CSC_DIR).

ChromiumPBMC-MetaCell-32bp/ now holds Bolero's base-resolution training signal for this dataset: a metacell × 32 bp coverage matrix for every genomic window, random-accessible by region. This random-access layer is exactly what the training and inference data loaders (RayGenomeChunkDataset, GenericGenomeDataManager) are built on.