Mouse brain scRNA-seq end to end (cytome)

The same analysis as Mouse brain scRNA-seq end to end (AnnData), on the same data, with one difference: the matrix is never loaded into memory. PIASO reads it from the .cytome file in chunks and writes results back onto the file.

Compare the two side by side — the calls are nearly identical. That is the point: switching backends is not a rewrite.

import cytome
import piaso
import cosg

piaso.settings.set_figure_params(style="cell")

1. Convert Cell Ranger output to a cytome

cytome.from_10x_h5 writes the file once; every later step reads from it.

fetch_dataset downloads the file once into ~/.piaso/data/datasets/, checks it against a published md5, and returns the path. Running it again on a complete file is instant; a download interrupted half way is re-fetched rather than silently used.

piaso.data.list_datasets()                              # what is available
h5 = piaso.data.fetch_dataset("mouse_brain_10k_gemx")   # ~65 MB, cached
ds_path = "mouse_brain_10k.cytome"

cytome.from_10x_h5(h5, ds_path, sample_name="mouse_brain_10k", force=True)
ds = cytome.open(ds_path)
ds

The conversion reads the matrix in cell batches rather than loading it, so peak memory is set by the batch size and not by the size of the file. batch_size is chosen automatically from the file’s own density; pass it explicitly to override.

To use your own Cell Ranger output instead, point from_10x_h5 at it — nothing else in this tutorial changes.

The dataset is a file on disk. Opening it reads metadata only — not the matrix.

ds.n_cells, ds.n_genes, sorted(ds.modalities)   # metadata only, no matrix read

2. Quality control

The same functions as the AnnData path, pointed at the dataset. Metrics are computed by streaming and stored in the cells table.

piaso.pp.calculateCellMetrics(ds, modality="RNA",
                              prefix_vars={"mt": "mt-", "ribo": ["Rps", "Rpl"]})
ds.cells.to_pandas()[["n_counts", "n_genes", "pct_counts_mt",
                      "pct_counts_ribo"]].describe()

piaso.pl.scatter(ds, x="n_counts", y="n_genes", color="pct_counts_ribo",
                 logx=True, logy=True, marginals=True)

3. Doublets

scrublet streams the matrix and writes scrublet_score and is_doublet back into the cells table:

piaso.pp.scrublet(ds, expected_doublet_rate=0.06, random_state=0)

piaso.pl.scatter(ds, x="n_counts", y="scrublet_score", color="is_doublet",
                 logx=True)

The automatic threshold calls 34 of 11,357; a > 0.3 cut marks 511. As on the AnnData path, the fixed cut is the one to use — and the cells it marks are removed together with any low-quality clusters, after the first clustering.

4. Normalization and feature selection

INFOG streams the counts and stores the normalized layer on the file, so the next step reads RNA_infog rather than recomputing it.

piaso.tl.infog(ds, modality="RNA", n_top_genes=3000, save_layer=True)

save_layer=True writes the full normalized matrix. That makes the file much larger — 92 MB to 448 MB here — and every later step reads it instead of recomputing. save_layer=False keeps the file small and normalizes on the fly.

5. Embedding and clusters

piaso.tl.runSVD(ds, modality="RNA", layer="infog",
                n_components=50, key_added="X_svd")

piaso.tl.neighbors(ds, use_rep="X_svd", n_neighbors=15)
piaso.tl.leiden(ds, resolution=1.0, key_added="leiden")
piaso.tl.umap(ds, use_rep="X_svd")

Everything lands on the file, not in memory:

ds.list_embeddings()

piaso.pl.embedding(ds, basis="X_umap", color="leiden", legend_loc="both")

40 clusters on the full 11,357 cells.

piaso.pl.plot_features_violin(
    ds, ["n_genes", "n_counts", "pct_counts_mt", "pct_counts_ribo",
         "scrublet_score"],
    groupby="leiden")

Three clusters stand out, all for the same reason: median scrublet_score 0.350 with 88–94% of their cells above 0.2. They are doublet clusters, and because this path clustered before removing the flagged cells, they are larger and more obvious here than on the AnnData path.

Note the score is coarsely quantized — 143 distinct values across 11,357 cells. That is Scrublet, not the cytome backend: the score is the fraction of a cell’s k nearest neighbours that are simulated doublets, so it can only take k+1 values.

6. Marker genes

cosg.cosg takes the dataset directly and reads it in chunks. On a cytome it returns a dict of arrays rather than writing into uns, since there is no uns to write into:

markers = cosg.cosg(ds, groupby="leiden", modality="RNA", n_genes_user=25,
                    layer="infog")
markers["names"].shape          # (n_genes_user, n_clusters)
markers["groups_order"][:5]     # which column is which cluster

top3 = []
for j, clus in enumerate(markers["groups_order"]):
    top3 += [g for g in markers["names"][:, j][:3] if g not in top3]

piaso.pl.dotplot(ds, top3, groupby="leiden", modality="RNA",
                 cytome_layer="infog", standard_scale="var")

dotplot reads the plotted features in a single streaming pass over the file, so the cost is set by how many genes you ask for, not by how many cells there are.

7. Remove the flagged cells, then re-run

filter_cells writes a new file rather than shrinking this one — a cytome is a file, and subsetting it is a copy. Pass inplace=False with an output path; inplace=True, output=None replaces the file atomically instead.

import numpy as np
cells = ds.cells.to_pandas()
keep = (~cells["leiden"].astype(str).isin(["11", "19", "28"])
        & ~(cells["is_doublet"].astype(bool) | (cells["scrublet_score"] > 0.3)))

piaso.pp.filter_cells(ds, mask=np.asarray(keep), inplace=False,
                      output="mouse_brain_10k_clean.cytome", overwrite=True)
ds.close()
ds = cytome.open("mouse_brain_10k_clean.cytome")

Then run the same four steps again, for the same reason as on the AnnData path: the gene selection and the SVD were fitted to the cell set that included the cells you just removed.

piaso.tl.infog(ds, modality="RNA", n_top_genes=3000, save_layer=True)
piaso.tl.runSVD(ds, modality="RNA", layer="infog", n_components=50,
                key_added="X_svd")
piaso.tl.neighbors(ds, use_rep="X_svd", n_neighbors=15)
piaso.tl.leiden(ds, resolution=1.0, key_added="leiden")
piaso.tl.umap(ds, use_rep="X_svd")

10,712 cells and 38 clusters, against 11,357 and 40 before.

piaso.pl.embedding(ds, basis="X_umap", color="leiden", legend_loc="both")

8. Cell types, not cluster numbers

The annotation step is the same call as in the AnnData tutorial, pointed at the dataset. Marker sets still come from an annotated reference held in memory — it is small, and only its markers cross over:

ref = piaso.data.load_dataset("adult_cortex_multiome_rna")
# .X is scaled (and carries NaN); the raw UMIs are in layers['raw'], so name it
piaso.tl.infog(ref, layer="raw", n_top_genes=3000)
cosg.cosg(ref, groupby="CellTypes", key_added="cosg", n_genes_user=30,
          layer="infog")
marker_sets = {ct: list(ref.uns["cosg"]["names"][ct])
               for ct in ref.obs["CellTypes"].cat.categories}
del ref

piaso.tl.predictCellTypeByMarker(ds, marker_gene_set=marker_sets,
                                 modality="RNA", cytome_layer="infog",
                                 score_layer="infog", use_rep="X_svd",
                                 key_added="CellTypes_pred")
ds.cells.to_pandas()["CellTypes_pred"].value_counts().head()

Scoring streams the matrix in two passes and the labels are written back into the cells table, so the annotation is on the file rather than in a variable.

piaso.pl.embedding(ds, basis="X_umap", color="CellTypes_pred")

top_ct = []
for ct in sorted(set(ds.cells.to_pandas()["CellTypes_pred"])):
    top_ct += [g for g in marker_sets[ct] if g not in top_ct][:1]

piaso.pl.dotplot(ds, top_ct, groupby="CellTypes_pred", modality="RNA",
                 cytome_layer="infog", standard_scale="var")

Reference markers by predicted cell type

9. The file is the result

Everything computed above is on disk. Reopening gives it all back, with no recomputation and no export step:

ds.close()

ds = cytome.open(ds_path)
ds.list_embeddings(), list(ds.cells.to_pandas().columns)[:8]

To hand the result to an AnnData-based tool:

# adata = ds.to_anndata(modality="RNA")

When to use which

	AnnData	cytome
Matrix location	memory	disk, read in chunks
Peak memory	scales with cells	set by batch size
Good for	up to ~10⁵ cells	10⁶+ cells, repeated analyses
Results	in the object	persisted on the file

For this 11k-cell dataset either is fine, and AnnData is simpler. The cytome path is what lets the same code run when the dataset is a hundred times larger.