1.3 Million brain cells on a laptop GPU

Mirrors the rapids-singlecell 1M-brain tutorial — the same workflow that needs a datacenter GPU there, run here on an Apple M-series laptop via metal-SingleCell. We process 1,000,000 cells end-to-end. (Expect a few minutes; the neighbor graph is the slow step.)

Requires the 10x 1.3M-neuron file at data/external/1M_neurons.h5.

import warnings; warnings.filterwarnings('ignore')
import time, gc, numpy as np, scanpy as sc
import metalsinglecell as msc
t0 = time.time()

Load and subset to 1M cells

import os
def _find(rel):              # resolve path from notebook dir up to the repo root
    d = os.getcwd()
    for _ in range(6):
        if os.path.exists(os.path.join(d, rel)): return os.path.join(d, rel)
        d = os.path.dirname(d)
    raise FileNotFoundError(rel + ' (place the 10x 1.3M-neuron .h5 at data/external/)')
adata = sc.read_10x_h5(_find('data/external/1M_neurons.h5'))
adata.var_names_make_unique()
adata = adata[:1_000_000].copy()
adata.shape

(1000000, 27998)

QC and filtering

adata.var['mt'] = adata.var_names.str.lower().str.startswith('mt-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True, percent_top=None)
msc.pp.filter_genes(adata, min_cells=3)
adata = adata[(adata.obs.n_genes_by_counts > 500) & (adata.obs.pct_counts_mt < 20)].copy()
gc.collect(); adata.shape

(986434, 22597)

Normalize, log1p, highly variable genes (Seurat v3 on counts)

adata.layers['counts'] = adata.X.copy()
msc.pp.normalize_total(adata, target_sum=1e4)
msc.pp.log1p(adata)
msc.pp.highly_variable_genes(adata, n_top_genes=5000, flavor='seurat_v3', layer='counts')
adata = adata[:, adata.var.highly_variable].copy()
del adata.layers['counts']; gc.collect(); adata.shape

(986434, 5000)

Sparse PCA, neighbors, UMAP, Leiden

We run PCA directly on the sparse log-normalized HVG matrix (no densifying scale/regress_out, which would need ~20 GB at this scale).

msc.pp.pca(adata, n_comps=50, use_highly_variable=False)
adata.obsm['X_pca'].shape

(986434, 50)

msc.pp.neighbors(adata, n_neighbors=15)
msc.tl.leiden(adata, resolution=1.0, backend='gpu')
adata.obs.leiden.nunique()

31

msc.tl.umap(adata, min_dist=0.3)
adata.obsm['X_umap'].shape

(986434, 2)

print(f'{adata.n_obs:,} cells · {adata.obs.leiden.nunique()} Leiden clusters · total wall time {time.time() - t0:.0f}s')

986,434 cells · 31 Leiden clusters · total wall time 697s

sc.pl.umap(adata, color=['leiden'], legend_loc=None, size=2, frameon=False)

A full 1M-cell clustering workflow — normalize → HVG → PCA → kNN graph → Leiden → UMAP — on a laptop GPU, with the Metal parallel Leiden doing the clustering.