{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bed36601",
   "metadata": {},
   "source": [
    "# 1.3 Million brain cells on a laptop GPU\n",
    "\n",
    "Mirrors the [rapids-singlecell 1M-brain tutorial](https://rapids-singlecell.readthedocs.io/en/latest/notebooks/brain_1M.html) — 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.)\n",
    "\n",
    "> Requires the 10x 1.3M-neuron file at `data/external/1M_neurons.h5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5525a9f8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:23:50.182260Z",
     "iopub.status.busy": "2026-06-23T18:23:50.181974Z",
     "iopub.status.idle": "2026-06-23T18:23:51.279222Z",
     "shell.execute_reply": "2026-06-23T18:23:51.278645Z"
    }
   },
   "outputs": [],
   "source": [
    "import warnings; warnings.filterwarnings('ignore')\n",
    "import time, gc, numpy as np, scanpy as sc\n",
    "import metalsinglecell as msc\n",
    "t0 = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac79fca",
   "metadata": {},
   "source": [
    "## Load and subset to 1M cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b34f1e00",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:23:51.280622Z",
     "iopub.status.busy": "2026-06-23T18:23:51.280456Z",
     "iopub.status.idle": "2026-06-23T18:25:36.462842Z",
     "shell.execute_reply": "2026-06-23T18:25:36.458866Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000000, 27998)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "def _find(rel):              # resolve path from notebook dir up to the repo root\n",
    "    d = os.getcwd()\n",
    "    for _ in range(6):\n",
    "        if os.path.exists(os.path.join(d, rel)): return os.path.join(d, rel)\n",
    "        d = os.path.dirname(d)\n",
    "    raise FileNotFoundError(rel + ' (place the 10x 1.3M-neuron .h5 at data/external/)')\n",
    "adata = sc.read_10x_h5(_find('data/external/1M_neurons.h5'))\n",
    "adata.var_names_make_unique()\n",
    "adata = adata[:1_000_000].copy()\n",
    "adata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1152244d",
   "metadata": {},
   "source": [
    "## QC and filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "453df773",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:25:36.560853Z",
     "iopub.status.busy": "2026-06-23T18:25:36.560333Z",
     "iopub.status.idle": "2026-06-23T18:26:48.024570Z",
     "shell.execute_reply": "2026-06-23T18:26:48.008577Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(986434, 22597)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata.var['mt'] = adata.var_names.str.lower().str.startswith('mt-')\n",
    "sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True, percent_top=None)\n",
    "msc.pp.filter_genes(adata, min_cells=3)\n",
    "adata = adata[(adata.obs.n_genes_by_counts > 500) & (adata.obs.pct_counts_mt < 20)].copy()\n",
    "gc.collect(); adata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb8f0780",
   "metadata": {},
   "source": [
    "## Normalize, log1p, highly variable genes (Seurat v3 on counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1b32768a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:26:48.050881Z",
     "iopub.status.busy": "2026-06-23T18:26:48.050735Z",
     "iopub.status.idle": "2026-06-23T18:30:43.591388Z",
     "shell.execute_reply": "2026-06-23T18:30:43.589703Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(986434, 5000)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata.layers['counts'] = adata.X.copy()\n",
    "msc.pp.normalize_total(adata, target_sum=1e4)\n",
    "msc.pp.log1p(adata)\n",
    "msc.pp.highly_variable_genes(adata, n_top_genes=5000, flavor='seurat_v3', layer='counts')\n",
    "adata = adata[:, adata.var.highly_variable].copy()\n",
    "del adata.layers['counts']; gc.collect(); adata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fae5046",
   "metadata": {},
   "source": [
    "## Sparse PCA, neighbors, UMAP, Leiden\n",
    "We run PCA directly on the sparse log-normalized HVG matrix (no densifying `scale`/`regress_out`, which would need ~20 GB at this scale)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "45c29ff3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:30:43.615165Z",
     "iopub.status.busy": "2026-06-23T18:30:43.615025Z",
     "iopub.status.idle": "2026-06-23T18:30:54.859463Z",
     "shell.execute_reply": "2026-06-23T18:30:54.854762Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(986434, 50)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msc.pp.pca(adata, n_comps=50, use_highly_variable=False)\n",
    "adata.obsm['X_pca'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a490224d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:30:54.864915Z",
     "iopub.status.busy": "2026-06-23T18:30:54.864710Z",
     "iopub.status.idle": "2026-06-23T18:32:13.284687Z",
     "shell.execute_reply": "2026-06-23T18:32:13.283855Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msc.pp.neighbors(adata, n_neighbors=15)\n",
    "msc.tl.leiden(adata, resolution=1.0, backend='gpu')\n",
    "adata.obs.leiden.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "08904c84",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:32:13.288065Z",
     "iopub.status.busy": "2026-06-23T18:32:13.287946Z",
     "iopub.status.idle": "2026-06-23T18:35:28.371122Z",
     "shell.execute_reply": "2026-06-23T18:35:28.370647Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(986434, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msc.tl.umap(adata, min_dist=0.3)\n",
    "adata.obsm['X_umap'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7b080459",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:35:28.373662Z",
     "iopub.status.busy": "2026-06-23T18:35:28.373566Z",
     "iopub.status.idle": "2026-06-23T18:35:28.378112Z",
     "shell.execute_reply": "2026-06-23T18:35:28.377720Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "986,434 cells · 31 Leiden clusters · total wall time 697s\n"
     ]
    }
   ],
   "source": [
    "print(f'{adata.n_obs:,} cells · {adata.obs.leiden.nunique()} Leiden clusters · total wall time {time.time() - t0:.0f}s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0355b7e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-23T18:35:28.381151Z",
     "iopub.status.busy": "2026-06-23T18:35:28.381072Z",
     "iopub.status.idle": "2026-06-23T18:35:28.790189Z",
     "shell.execute_reply": "2026-06-23T18:35:28.789547Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['leiden'], legend_loc=None, size=2, frameon=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fad2dc0f",
   "metadata": {},
   "source": [
    "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}