Mouse brain scRNA-seq end to end (AnnData)

From a Cell Ranger .h5 to clusters, cell types and marker genes. Everything here runs on a plain pip install piaso-tools — no scanpy, no scrublet.

The companion tutorial, Mouse brain scRNA-seq end to end (cytome), does the same analysis without loading the matrix into memory. The calls are almost identical; the difference is where the data lives.

For a human sample — where the mitochondrial thresholds actually bite — see Human PBMC scRNA-seq end to end (AnnData).

import numpy as np
import piaso
import cosg

# Journal preset: square single-column panels at Cell's width and base font.
# 'nature', 'nature_methods', 'science' and others are also available.
piaso.settings.set_figure_params(style="cell")

Every plot below is shown inline. To write one to disk, pass save=:

# piaso.pl.embedding(adata, basis="X_umap", color="leiden", save="umap.png")

piaso.settings.figdir sets where those go.

1. Get the data

piaso.data keeps a small registry of example datasets and caches downloads, so the same call works on a laptop and on a cluster.

piaso.data.list_datasets()

We use a 10x Genomics mouse brain sample: 11,357 nuclei, 33,696 genes, Cell Ranger filtered feature-barcode matrix (~65 MB).

adata = piaso.data.load_dataset("mouse_brain_10k_gemx")
adata

load_dataset downloads the file once into ~/.piaso/data/datasets/, verifies it against a published md5, then opens it. Running it again is instant; a download interrupted half way is re-fetched rather than silently used. To download without opening — useful on a cluster login node — use piaso.data.fetch_dataset(...), which returns the path.

load_dataset reads Cell Ranger HDF5 with piaso.pp.read_10x_h5. If you have your own file, that function is the entry point:

# adata = piaso.pp.read_10x("path/to/filtered_feature_bc_matrix.h5")
# adata = piaso.pp.read_10x("path/to/filtered_feature_bc_matrix/")   # MTX dir

.X holds raw UMI counts — that is what the next step expects.

2. Quality control

calculateCellMetrics computes per-cell totals. prefix_vars adds the percentage of counts in any gene-name prefix — mitochondrial genes are mt- in mouse, MT- in human.

piaso.pp.calculateCellMetrics(adata, prefix_vars={"mt": "mt-",
                                                  "ribo": ["Rps", "Rpl"]})
adata.obs[["n_counts", "n_genes", "pct_counts_mt", "pct_counts_ribo"]].describe()

prefix_vars takes any number of prefix groups. Mitochondrial genes start with mt- in mouse and MT- in human; ribosomal protein genes start with Rps/Rpl and RPS/RPL. Both are worth having: a cell that is mostly mitochondrial or mostly ribosomal reads is usually a cell that has lost its cytoplasm.

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

One plot, three questions: the depth–complexity relation, the marginal distribution of each, and where the mitochondrial reads sit. Cells falling below the main band with high pct_counts_mt are the ones to cut.

This is a nuclei sample, so both fractions are near zero — median pct_counts_mt 0.011% and median pct_counts_ribo 0.20% — and neither can do the work it does for whole cells. The gene count carries the signal here instead. That is worth knowing before you set a mitochondrial threshold that would remove nothing. (On whole cells it is completely different: the human PBMC walkthrough has a median of 12%, and a whole cluster there turns out to be mitochondrial reads.)

Filter both tails, not only the low one

The convention is to drop cells with too few counts. Too many is also a signal — a barcode with several times the median amount of RNA is often two nuclei in one droplet — and on this sample the upper tail is where the cells are:

hi_c = np.percentile(adata.obs["n_counts"], 99)          # 69,373
ribo_cut = np.percentile(adata.obs["pct_counts_ribo"], 99)   # 0.83%

keep = ((adata.obs["n_counts"] >= 500)
        & (adata.obs["n_counts"] <= hi_c)
        & (adata.obs["n_genes"] >= 250)
        & (adata.obs["pct_counts_mt"] <= 5.0)
        & (adata.obs["pct_counts_ribo"] <= ribo_cut))
adata = adata[keep].copy()

n_counts < 500                     0
n_counts > 69,373 (99th pct)     114
n_genes < 250                      3
pct_counts_mt > 5.0               27
pct_counts_ribo > 0.83 (99th)    114
                     union      233 of 11,357  (2.1%)

Not one cell fails the low-count floor. Cell Ranger’s filtered matrix has already dropped empty droplets, so the usual first filter does nothing here — and the 114 cells at the top of the count distribution, which no conventional threshold touches, are the ones worth removing. The mitochondrial cut at 5% takes 27 cells; on this sample it is close to a no-op, which is the honest outcome for a nuclei prep, not a reason to lower it until it bites.

scatter defaults to magma_r for a continuous colour: sequential, perceptually uniform and colourblind-safe. cmap='Spectral_r' or 'RdBu_r' gives the diverging look.

3. Doublets

piaso.pp.scrublet is PIASO’s own implementation of the Scrublet algorithm — no extra install. It simulates doublets from the observed cells and scores each barcode by how close it sits to that simulated population.

piaso.pp.scrublet(adata, expected_doublet_rate=0.06, random_state=0)
adata.obs["is_doublet"].sum(), adata.n_obs

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

On this sample the automatic threshold calls 28 of 11,357 barcodes (0.25%), far below the 6% prior — the prior is what the simulation uses, not what it must find.

That automatic threshold is conservative. It looks for a dip between two modes in the score distribution, and on a clean sample there is no clear dip to find. A fixed cut on the score itself is the usual practice, and it removes an order of magnitude more:

adata.obs["doublet_call"] = (adata.obs["is_doublet"]
                             | (adata.obs["scrublet_score"] > 0.3))
adata.obs["doublet_call"].sum()

476 cells here, against 28 from the automatic call. 0.3 is a convention, not a constant — read it off the scatter above, and off the per-cluster violins in step 6.

For a multi-sample object, pass library_key= so each library is scored against its own background:

# piaso.pp.scrublet(adata, library_key="sample_id")

adata = adata[~adata.obs["doublet_call"]].copy()
piaso.pp.filter_cells(adata, min_counts=500, min_features=250)
adata.shape

4. Normalization and feature selection

INFOG normalizes and selects informative genes in one step. It reads raw UMI counts: from adata.X by default, or from a layer if you pass layer=.

piaso.tl.infog(adata, n_top_genes=3000)

It wrote adata.layers['infog'] and adata.var['highly_variable'].

5. Dimensionality reduction

Pass layer="infog" explicitly. runSVD defaults to adata.X, which would run the SVD on the raw counts and silently ignore step 4.

piaso.tl.runSVD(adata, layer="infog", n_components=50, key_added="X_svd")

6. Neighbours, clusters, UMAP

piaso.tl.neighbors(adata, use_rep="X_svd", n_neighbors=15)

piaso.tl.leiden(adata, resolution=1.0, key_added="leiden")
adata.obs["leiden"].value_counts().head()

piaso.tl.umap(adata, use_rep="X_svd")

39 clusters at resolution=1.0. That is a property of the data, not a setting to fix — brain samples split finely because the neuronal subtypes really are distinct.

legend_loc decides where the cluster names go. On the data while you are working out what each cluster is:

piaso.pl.embedding(adata, basis="X_umap", color="leiden", legend_loc="on_data")

Beside the panel when you want to read a colour off it:

piaso.pl.embedding(adata, basis="X_umap", color="leiden", legend_loc="right")

Or both at once:

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

7. Which clusters are artefacts?

Per-cell filtering does not catch everything. A cluster can be made entirely of doublets, or of empty droplets carrying only ambient RNA, and still survive: every one of its cells looks unremarkable next to the whole sample. It only stands out next to the other clusters.

Start with the QC metrics per cluster:

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

And the marker genes, three per cluster:

cosg.cosg(adata, groupby="leiden", key_added="cosg", n_genes_user=25,
          layer="infog")

top3 = []
for c in adata.obs["leiden"].cat.categories:
    top3 += [g for g in adata.uns["cosg"]["names"][c][:3] if g not in top3]

piaso.pl.dotplot(adata, top3, groupby="leiden", standard_scale="var")

Reading the dot plot

A good cluster’s own markers are both darker and larger there than anywhere else — a high mean and a high fraction of cells expressing.

Colour alone is not enough, and this is where a low-quality cluster is easy to misread. Normalization divides by each cell’s total, so a cluster with a low total count can show a high normalized mean for a gene that is not particularly its own. What gives it away is the dot size: the same gene is expressed in a large fraction of the other clusters too. A cluster whose top three genes are large dots everywhere has no specificity of its own, which is what empty droplets with ambient RNA, and other low-quality clusters, look like.

Why not a simple gene-count floor

A rule like “median genes below half the sample median” looks reasonable and is wrong on this data: glia genuinely have far fewer detected genes than neurons, so it flags real biology. Compare each cluster to the spread of clusters, and require a second, independent signal before dropping anything:

qc = adata.obs.groupby("leiden", observed=True).agg(
    n_cells=("n_genes", "size"),
    median_genes=("n_genes", "median"),
    median_pct_mt=("pct_counts_mt", "median"),
    median_pct_ribo=("pct_counts_ribo", "median"),
    median_scrublet=("scrublet_score", "median"),
    frac_high_scrublet=("scrublet_score", lambda v: float((v > 0.2).mean())),
)

# how specific each cluster's best marker is at all
import pandas as pd
names = pd.DataFrame(adata.uns["cosg"]["names"])
scores = pd.DataFrame(adata.uns["cosg"]["scores"])
qc["top_cosg"] = [float(scores[c].iloc[0]) for c in qc.index]

# and whether the markers are themselves technical
qc["mt_in_top10"] = [sum(g.startswith("mt-") for g in names[c].head(10))
                     for c in qc.index]
qc["ribo_in_top10"] = [sum(g.startswith(("Rps", "Rpl")) for g in names[c].head(10))
                       for c in qc.index]

Does a high-mitochondrial cluster mean a bad cluster?

Worth checking rather than assuming. Sorting the clusters by median pct_mt:

leiden	n	median genes	median mt%	top COSG	top markers
2	106	3,675	0.127	0.882	Bnc2, Slc47a1, Tspan11
37	88	3,834	0.123	0.790	Drd3, Prkar2b, Gpr26
28	91	2,120	0.092	0.920	Flt1, Adgrl4, Ptprb, Slco1a4
39	66	5,244	0.068	0.639	Vipr2, Avp, Dlk1

mt_in_top10 is zero for every one of the 40 clusters. Not a single mitochondrial gene appears among any cluster’s ten best markers, and the whole mt% column spans 0.05–0.13% — a range in which the ordering is noise.

Cluster 28 is the interesting one: it has the fewest detected genes in the object (2,120, well below the sample median) and it would fail any gene-count floor. Its markers are Flt1, Adgrl4, Ptprb, Slco1a4 — endothelial cells, and its top marker scores 0.920, the highest specificity of any cluster here. Endothelial nuclei really do carry less RNA. Dropping it on the gene count would have thrown away the vasculature.

So on this sample, high mitochondrial content identifies nothing, and the check that would matter — are the markers mitochondrial genes? — comes back negative everywhere. It is still worth running, because on whole cells the same check does fire: see §4 of the human PBMC page, where one cluster has five MT genes in its top ten and a top marker score of 0.028.

# position within the spread of clusters, in interquartile units
lo, hi = qc["median_genes"].quantile([0.25, 0.75])
qc["genes_z"] = (qc["median_genes"] - qc["median_genes"].median()) / (hi - lo)

qc["doublet_cluster"] = ((qc["median_scrublet"] > 0.15)
                         | (qc["frac_high_scrublet"] > 0.5))
qc["ambient_cluster"] = (qc["mt_in_top10"] >= 3) | (qc["ribo_in_top10"] >= 3)
qc["low_quality"] = ((qc["genes_z"] < -1.5)
                     & (qc["top_cosg"] < qc["top_cosg"].quantile(0.1)))
qc[qc["low_quality"] | qc["doublet_cluster"] | qc["ambient_cluster"]]

Two clusters fail, both on doublets:

	n	median counts	median genes	genes_z	top COSG	median scrublet	> 0.2	verdict
29	35	23,443	5,292	+0.39	0.010	0.379	1.00	doublet cluster
33	235	32,395	6,083	+0.73	0.033	0.325	0.83	doublet cluster
19	185	8,856	3,099	−0.53	0.479	0.137	0.48	kept, borderline
28	91	5,900	2,120	−1.42	0.920	0.023	0.02	kept — endothelium

Read the two extremes against each other. Cluster 29 has an above-average gene count and every single cell above the doublet threshold, and its best marker scores 0.010 — it has no identity at all. Cluster 28 has by far the fewest genes and a top marker at 0.920. Gene count alone would have kept the first and dropped the second, which is exactly backwards.

Cluster 19 is the honest borderline case: 48% of its cells score above 0.2, just under the 0.5 rule, and its top marker scores 0.479 — respectable but not strong. Its markers (Arhgap15, Ifi207, Ifi202b, Aoah, Mndal) are a coherent immune signature, so it is kept here. A defensible pipeline would flag it and look, rather than let the threshold decide silently.

ambient_cluster fires on nothing in this sample, as §7 showed. Keep the rule anyway — it costs one line and it is the rule that catches a dying-cell cluster on whole-cell data.

adata = adata[~adata.obs["leiden"].isin(
    qc.index[qc["low_quality"] | qc["doublet_cluster"]
             | qc["ambient_cluster"]])].copy()
adata.shape

(10854, 33696)      # from 11,124 after the per-cell filters

8. Re-run the pipeline on the cells that remain

This matters and is easy to skip. infog picked its 3,000 genes from the old cell set, the SVD was fitted to it, and the neighbour graph and UMAP were built on that fit. Keeping them means analysing the clean cells in a space defined partly by the cells you just removed. So run the same four steps again:

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

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

10,849 cells, and 39 clusters become 38 — the removed cluster does not simply vanish, the whole partition shifts a little.

cosg.cosg(adata, groupby="leiden", key_added="cosg", n_genes_user=25,
          layer="infog")
top3 = []
for c in adata.obs["leiden"].cat.categories:
    top3 += [g for g in adata.uns["cosg"]["names"][c][:3] if g not in top3]

piaso.pl.dotplot(adata, top3, groupby="leiden", standard_scale="var")

standard_scale="var" scales each gene to its own 0–1 range across clusters, so a specific but lowly expressed gene is as visible as a loud one. The colour map is Spectral_r; pass cmap="Reds" or "magma_r" for a colourblind-safe sequential scale in a manuscript figure.

9. Cell types, not cluster numbers

Clusters are unnamed. predictCellTypeByMarker scores every cell against a set of marker genes per cell type, assigns the best-scoring type, then smooths the assignment over each cell’s nearest neighbours so single-cell noise does not produce isolated labels.

The marker sets can come from either of two places. They are shown side by side below because they fail differently, and comparing them is itself a check.

Route A — PIASOmarkerDB

PIASOmarkerDB ships curated marker sets from published studies, queryable by species, tissue and study:

piaso.tl.getMarkers(list_studies=True)

For a mouse brain sample, AllenWholeMouseBrain_isocortex is the matched study:

# as_dict=True returns BOTH: the marker table and a
# {cell type: [genes]} dictionary. Unpack them.
markers_df, marker_db = piaso.tl.getMarkers(
    study="AllenWholeMouseBrain_isocortex", as_dict=True)
len(marker_db), list(marker_db)[:4]

piaso.tl.predictCellTypeByMarker(adata, marker_gene_set=marker_db,
                                 score_layer="infog", use_rep="X_svd",
                                 key_added="CellTypes_db")
adata.obs["CellTypes_db"].value_counts().head()

use_rep is the embedding used for smoothing; it defaults to X_gdr, so pass X_svd when the pipeline above is what you ran.

piaso.pl.embedding(adata, basis="X_umap", color="CellTypes_db")

The labels keep the Allen taxonomy’s names (327 Oligo NN, 007 L2/3 IT CTX Glut), which is a feature: they are traceable to the study they came from. Each marker set carries its own citation in the database.

Route B — an annotated reference dataset

When no published set matches, derive one from an annotated dataset of the same tissue with COSG. Here that is the adult mouse cortex multiome, which carries 20 curated types:

ref = piaso.data.load_dataset("adult_cortex_multiome_rna")
ref.obs["CellTypes"].value_counts()

That dataset’s .X is not raw counts — it is scaled (93% of its values are negative, and 349 genes carry NaN). load_dataset warns about this on load and names the layer that does hold raw UMIs. Point infog at that layer:

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_ref = {ct: list(ref.uns["cosg"]["names"][ct])
              for ct in ref.obs["CellTypes"].cat.categories}
del ref                                # 17k × 26k reference, no longer needed

layer= is the better habit: it leaves .X as the author stored it. Copying the counts over .X gives byte-identical results — same layers['infog'] to 0.0, same highly_variable — but it destroys the original. This applies to anything that reads counts: infog, run_TFIDF and score all take layer=.

piaso.tl.predictCellTypeByMarker(adata, marker_gene_set=marker_ref,
                                 score_layer="infog", use_rep="X_svd",
                                 key_added="CellTypes_ref")
piaso.pl.embedding(adata, basis="X_umap", color="CellTypes_ref")

Do the two routes agree?

from sklearn.metrics import adjusted_rand_score
adjusted_rand_score(adata.obs["CellTypes_db"], adata.obs["CellTypes_ref"])

ARI 0.881, and 83.6% of cells fall in the same best-matching pair. The agreement is near-total on the broad classes — oligodendrocytes 4,083/4,083, microglia 942/957, OPCs 941/953 — and breaks down on the fine excitatory subtypes, where the two taxonomies do not draw the same boundaries. Two independent marker sources agreeing on the glia is worth more than either one alone; where they disagree is where the label is a guess.

10. Check the annotation against its own markers

Plot the marker genes that drove the prediction, grouped by the label it produced. A clean diagonal is evidence the annotation holds, rather than a restatement of it.

top_ct, seen = [], set()
for ct in adata.obs["CellTypes_ref"].cat.categories:
    picked = [g for g in marker_ref[ct] if g in adata.var_names and g not in seen][:3]
    top_ct += picked; seen.update(picked)

piaso.pl.dotplot(adata, top_ct, groupby="CellTypes_ref", standard_scale="var")

Reference markers by predicted cell type

The same check for route A:

PIASOmarkerDB markers by predicted cell type

The broad classes land cleanly. The fine excitatory subtypes are only partly separated, which is what to expect when the reference is cortex and the sample is whole brain — the label is a best match, not a measurement, and the dot plot is where you find that out.

11. Three views of the finished object

Once the clusters have names, three plots answer questions the UMAP cannot.

Per-group QC, as violins. The UMAP shows structure; this shows whether any group is structurally poor:

piaso.pl.violin(adata, ["n_genes", "n_counts", "pct_counts_mt",
                        "pct_counts_ribo"], groupby="leiden")

The same thing as a table, ranked. plotGroupMetrics takes any per-group DataFrame, so it will summarise whatever you can compute — here size and depth side by side, which is how you see that a small cluster is also a shallow one:

summary = adata.obs.groupby("leiden", observed=True).agg(
    n_cells=("n_genes", "size"),
    median_genes=("n_genes", "median"),
    median_counts=("n_counts", "median"),
    median_pct_mt=("pct_counts_mt", "median"),
)
piaso.pl.plotGroupMetrics(summary)

Composition. With more than one sample, split by sample and look for a cluster that is nearly all one library. This is a single sample, so there is nothing to split by — but the same plot against sequencing-depth tertile answers a related question, and it is the one worth asking on one library:

q = adata.obs["n_counts"].quantile([1/3, 2/3]).values
adata.obs["depth_tertile"] = pd.Categorical(
    np.select([adata.obs["n_counts"] <= q[0], adata.obs["n_counts"] <= q[1]],
              ["low", "mid"], "high"),
    categories=["low", "mid", "high"])

piaso.pl.stackedBarplot(adata, groupby="leiden", splitby="depth_tertile")

A cluster made almost entirely of one depth tertile is a depth artefact wearing a cluster’s clothes. On this dataset the tertiles are spread across clusters, which is the answer you want.

12. Per-cell-type QC, after annotation

The cells passed QC and the clusters passed QC. The last check is on the labels: a cell type that sits apart from the others on a technical metric was probably called on that metric.

gm = piaso.pp.calculateGroupMetrics(adata, groupby="CellTypes")
piaso.pl.plotGroupMetrics(gm, data=adata, groupby="CellTypes")

calculateGroupMetrics builds the per-group table (cell count, features detected, mean and median counts, median features per cell); plotGroupMetrics renders it as one panel per metric. Note the API: the plot function takes the table as its first argument, and data=/groupby= only so it can reuse the same cell-type colours as the UMAP.

13. GDR: rebuild the embedding from the cell types

Everything so far ran on X_svd — 50 components of 3,000 genes chosen by variance. Now that the cells have labels, GDR can build an embedding from the genes that distinguish those labels instead:

piaso.tl.runGDR(adata, batch_key=None, groupby="CellTypes",
                n_gene=30, mu=10, layer="infog", score_layer="infog")

piaso.tl.neighbors(adata, use_rep="X_gdr", n_neighbors=15)
piaso.tl.umap(adata, use_rep="X_gdr", key_added="X_umap_gdr")
piaso.pl.embedding(adata, basis="X_umap_gdr", color="CellTypes")

Whether it helped is measurable, not a matter of taste:

from sklearn.metrics import silhouette_score
for key in ("X_svd", "X_gdr"):
    print(key, silhouette_score(adata.obsm[key],
                                adata.obs["CellTypes"].astype(str),
                                sample_size=5000, random_state=0))

X_svd   0.386
X_gdr   0.402

A modest improvement, and modest is the honest description — on a clean, well-separated single sample the SVD embedding was already close to the achievable separation. GDR earns much more where variance-based selection struggles: many samples with batch effects, subtle sub-types, or a reference whose labels you want to project onto new data. Those cases are GDR, multiple samples and projectGDR.

14. Save

# adata.write_h5ad("mouse_brain_processed.h5ad")

Where to go next

Larger data — the cytome version runs the same steps by streaming from disk, so memory does not scale with cell count.
Blood rather than brain tissue — the Human PBMC walkthrough, where the mitochondrial thresholds that do nothing here do real work: median mitochondrial content is 12.0% there against 0.011% here.
Batch integration — piaso.tl.runHarmony, or runGDR(batch_key=...).
GDR in depth — GDR.