{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# 3. Preprocessing\n\n**Pipeline phase:** Phase 3\n\n**What you will learn:**\n- Back up raw counts before normalization\n- Normalize and log-transform\n- Select highly variable genes\n- Run PCA, neighbors, UMAP, and Leiden clustering\n- Use the modality-aware `preprocess` recipe\n\n**Note on integration:** The default integration for Visium is scVI.\nThis tutorial uses Harmony (lighter dependency) to stay dependency-free.\nSee the [API docs](../api/index.rst) for all `sc_tools.pp` options."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import anndata as ad\n",
    "import sc_tools.pp as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(7)\n",
    "n_spots, n_genes = 400, 300\n",
    "\n",
    "X = np.random.negative_binomial(4, 0.45, size=(n_spots, n_genes)).astype(float)\n",
    "gene_names = [f\"GENE_{i}\" for i in range(n_genes)]\n",
    "\n",
    "adata = ad.AnnData(\n",
    "    X=X,\n",
    "    obs={\n",
    "        \"sample\": [f\"S{i % 2}\" for i in range(n_spots)],\n",
    "        \"library_id\": [f\"lib_{i % 2}\" for i in range(n_spots)],\n",
    "    },\n",
    ")\n",
    "adata.var_names = gene_names\n",
    "print(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Option A: Step-by-step preprocessing\n",
    "\n",
    "This gives you full control over each step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Backup raw counts before any transformation\n",
    "pp.backup_raw(adata)\n",
    "print(\"adata.raw backed up:\", adata.raw is not None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Optionally filter mitochondrial / ribosomal genes before normalization\n",
    "pp.filter_genes_by_pattern(adata, patterns=[\"^MT-\", \"^RPL\", \"^RPS\"])\n",
    "print(f\"After MT/RP filter: {adata.shape}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. Normalize to 10,000 counts per spot, then log1p transform\n",
    "pp.normalize_total(adata, target_sum=1e4)\n",
    "pp.log_transform(adata)\n",
    "\n",
    "print(\"X max after log1p:\", adata.X.max().round(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. PCA on top-N genes (seurat_v3 HVG by default)\n",
    "pp.pca(adata, n_comps=30, n_top_genes=150)\n",
    "print(\"PCA done. obsm keys:\", list(adata.obsm.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5. Build neighbor graph and compute UMAP\n",
    "pp.neighbors(adata, n_neighbors=15)\n",
    "pp.umap(adata)\n",
    "print(\"UMAP done. Shape:\", adata.obsm[\"X_umap\"].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6. Leiden clustering\n",
    "pp.leiden(adata, resolution=0.5)\n",
    "print(\"Clusters:\", adata.obs[\"leiden\"].value_counts().to_dict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Option B: Recipe (one-call preprocessing)\n",
    "\n",
    "The `preprocess` recipe handles all steps above for a given modality.\n",
    "For IMC data, it uses `arcsinh_transform` instead of `log_transform`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-create fresh adata\n",
    "adata2 = ad.AnnData(X=X.copy())\n",
    "adata2.var_names = gene_names\n",
    "adata2.obs[\"library_id\"] = [f\"lib_{i % 2}\" for i in range(n_spots)]\n",
    "\n",
    "pp.preprocess(\n",
    "    adata2,\n",
    "    modality=\"visium\",\n",
    "    integration=\"harmony\",   # scvi is default; harmony used here for lighter deps\n",
    "    batch_key=\"library_id\",\n",
    "    n_top_genes=150,\n",
    "    n_pcs=30,\n",
    "    resolution=0.5,\n",
    ")\n",
    "\n",
    "print(\"Done. obsm:\", list(adata2.obsm.keys()))\n",
    "print(\"Clusters:\", adata2.obs[\"leiden\"].value_counts().to_dict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## Save Phase 3 checkpoint\n\n```python\nadata.write_h5ad(\"results/adata.normalized.p3.h5ad\")\n```\n\nThe checkpoint must have:\n- `adata.raw` backed up (raw counts)\n- A batch-corrected embedding (e.g. `obsm['X_scvi']` or `obsm['X_pca_harmony']`)\n- `obs['leiden']` cluster labels\n\nSee `Architecture.md §2.2` in the repository for the full requirement table."
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sc_tools",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}