Skip to content
PIASO

Human PBMC scRNA-seq end to end (cytome)

The same analysis as the AnnData version, with the matrix on disk instead of in memory. The calls are nearly identical; what changes is where the data lives and where each result is written.

Read the AnnData page for the reasoning — why both tails of the count distribution are filtered, why a cluster’s markers are its most informative QC metric. This page covers what is different about doing it on a file.

import numpy as np
import pandas as pd
import cytome
import piaso
import cosg


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

1. Convert once

h5 = piaso.data.fetch_dataset("pbmc_multiome_san2")   # ~88 MB, cached


ds = cytome.from_10x_h5(h5, "san2.cytome", sample_name="SAN2", force=True)
ds.n_cells, ds.n_genes, sorted(ds.modalities)
(4360, 36601, ['ATAC', 'RNA'])

Use the dataset from_10x_h5 returns. It hands back an open dataset. A common pattern is to ignore the return value and call cytome.open(path) on the next line — that leaves two writers on one SQLite file, and the next step fails with database is locked. One handle, from the function that made the file.

This is a multiome, so both modalities come across. Everything below passes modality="RNA"; the ATAC side stays on the file, untouched and available.

import os
os.path.getsize("san2.cytome") / 1e6
140.0     # MB, both modalities

2. QC by streaming

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

The metrics are computed in cell batches and written into the cells table — they are on the file, not in a variable, and reopening gives them back.

ds.cells.to_pandas() pulls that table into memory. It is one row per cell and a handful of columns, so it is small; the matrix is what stays on disk.

piaso.pl.scatter(ds, x="n_counts", y="n_genes", color="pct_counts_mt",
                 logx=True, logy=True, marginals=True)
QC scatter 
piaso.pp.scrublet(ds, expected_doublet_rate=0.06, random_state=0)

3. Normalize, embed, cluster

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

save_layer=True writes the normalized matrix onto the file so later steps read it instead of recomputing. That costs disk: 140 MB → 201 MB here. save_layer=False keeps the file small and normalizes on the fly — the right choice when disk is the constraint and the run is one-shot.

piaso.pl.embedding(ds, basis="X_umap", color="leiden", legend_loc="both")
First-pass clusters 
piaso.pl.plot_features_violin(
    ds, ["n_genes", "n_counts", "pct_counts_mt", "pct_counts_ribo",
         "scrublet_score"], groupby="leiden")
QC metrics per cluster

4. The same artefact check

On a cytome, cosg.cosg returns a dict of arrays — there is no uns to write into — so the columns have to be lined up with groups_order by hand:

markers = cosg.cosg(ds, groupby="leiden", modality="RNA", n_genes_user=25,
                    layer="infog")
order = list(markers["groups_order"])
names = pd.DataFrame(markers["names"], columns=order)
scores = pd.DataFrame(markers["scores"], columns=order)
cells = ds.cells.to_pandas()
qc = cells.groupby("leiden", observed=True).agg(
    n_cells=("n_genes", "size"),
    median_genes=("n_genes", "median"),
    median_pct_mt=("pct_counts_mt", "median"),
    median_scrublet=("scrublet_score", "median"),
    frac_high_scrublet=("scrublet_score", lambda v: float((v > 0.2).mean())))
qc = qc.loc[[c for c in order if c in qc.index]]     # match COSG's column order
qc["top_cosg"] = [float(scores[c].iloc[0]) for c in qc.index]
qc["mt_in_top10"] = [sum(str(g).startswith("MT-") for g in names[c].head(10))
                     for c in qc.index]


drop = list(qc.index[(qc["median_scrublet"] > 0.15)
                     | (qc["frac_high_scrublet"] > 0.5)
                     | (qc["mt_in_top10"] >= 3)])

Reindexing qc by order is not cosmetic. The groupby sorts cluster labels as strings and COSG returns its own order; zipping them without aligning attaches each cluster’s markers to a different cluster’s statistics, and the result looks entirely plausible.

15 clusters on the first pass, and the check finds two:

leidennmedian mt%top COSGMT in top 10scrublet > 0.2markers
1121515.90.11601%PID1, LYZ, SASH1
923715.40.69101%TCF7L2, FCGR3A
1230113.10.030464%MTRNR2L12, MT-…
48611.30.035049%ZNF227, …

This is an independent confirmation of the AnnData result. Different backend, different clustering (15 clusters rather than 19), and the same mitochondrial population turns up — four MT genes in its top ten, a best-marker score of 0.030, and 64% of its cells above the doublet threshold. And again the two highest-mt% clusters are not the bad ones: cluster 11 tops the column at 15.9% and cluster 9 is a clean non-classical monocyte population at 15.4%.

5. Filter writes a new file

hi_c = np.percentile(cells["n_counts"], 99)
mt_cut = np.percentile(cells["pct_counts_mt"], 98)


keep = (~cells["leiden"].astype(str).isin(drop)
        & (cells["n_counts"] >= 500) & (cells["n_counts"] <= hi_c)
        & (cells["n_genes"] >= 250)
        & (cells["pct_counts_mt"] <= mt_cut)
        & ~(cells["is_doublet"].astype(bool) | (cells["scrublet_score"] > 0.3)))


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

A cytome is a file, so subsetting it is a copy. inplace=False with an output path writes a new one and leaves the original intact; inplace=True, output=None replaces the file atomically. Both cell-level and cluster-level criteria go into one mask, so the copy happens once.

DROP: ['4', '12']
keeping 3515 of 4360   (counts <= 15,603, mt <= 21.4%)

Then the same four steps again, for the same reason as on the AnnData path — the gene selection and the SVD were fitted with the removed cells present:

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")
UMAP after cleanup

3,515 cells and 13 clusters, from 4,360 and 15.

markers = cosg.cosg(ds, groupby="leiden", modality="RNA", n_genes_user=25,
                    layer="infog")
order = list(markers["groups_order"])
names = pd.DataFrame(markers["names"], columns=order)


top3 = []
for c in order:
    top3 += [g for g in map(str, names[c][:3]) if g not in top3]


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

dotplot reads only the plotted genes, in one streaming pass, so its cost is set by how many genes you ask for rather than by how many cells there are.

6. Annotation, written back to the file

markers_df, marker_sets = piaso.tl.getMarkers(
    study="AllenHumanImmuneHealthAtlas_L2", as_dict=True)


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

Scoring streams the matrix and the labels land in the cells table. Nothing has to be exported for the annotation to persist — it is on the file the moment the call returns.

piaso.pl.embedding(ds, basis="X_umap", color="CellTypes")
Predicted cell types 
piaso.pl.stackedBarplot(ds, groupby="CellTypes", splitby="leiden")
Cluster composition per cell type 
gm = piaso.pp.calculateGroupMetrics(ds, groupby="CellTypes")
piaso.pl.plotGroupMetrics(gm, data=ds, groupby="CellTypes")
QC metrics per cell type

The plotting functions take the dataset wherever they took an AnnData. They stream what they need.

7. GDR

piaso.tl.runGDR(ds, batch_key=None, groupby="CellTypes", n_gene=30, mu=10,
                modality="RNA", layer="infog")


piaso.tl.neighbors(ds, use_rep="X_gdr", n_neighbors=15)
piaso.tl.umap(ds, use_rep="X_gdr", key_added="X_umap_gdr")
piaso.pl.embedding(ds, basis="X_umap_gdr", color=["CellTypes", "leiden"])
GDR embedding

runGDR writes X_gdr into the cytome and returns None, mirroring the way the AnnData path writes into .obsm in place. Pass write_to_cytome=False for the tuple return instead.

The parameter is layer=, not cytome_layer=. Several cytome-aware functions took cytome_layer in earlier releases and runGDR no longer does; passing it raises rather than being ignored, which is the behaviour you want.

8. The file is the result

ds.list_embeddings()
list(ds.cells.to_pandas().columns)[:10]
['RNA_svd', 'X_umap', 'X_gdr', 'X_umap_gdr']

Every embedding, label and QC column computed above is on disk. The file went 140 MB (conversion) -> 201 MB (with the INFOG layer) -> 158 MB after filtering, and that last file holds the counts, the normalized layer, four embeddings, the cluster labels, the predicted cell types and every QC metric — for both modalities. Reopening the file gives it all back with no recomputation and no export step. To hand it to an AnnData-based tool:

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

AnnData or cytome

AnnDatacytome
where the matrix livesmemorydisk
memory usescales with cells × genesset by batch size
resultsin the object, lost unless writtenon the file, immediately
subsettinga view or a copy in memorywrites a new file
second modalitya second objectalready there
best foranything that fits comfortablylarge data, or many sessions on one dataset

At 4,360 cells this sample fits in memory easily and the AnnData path is simpler. The point of running it here is that the calls are the same ones, so the step up to data that does not fit is a change of backend rather than a rewrite.

Next