{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. QC and Filtering\n",
    "\n",
    "**Pipeline phase:** Phase 1 — Data Ingestion and QC\n",
    "\n",
    "**What you will learn:**\n",
    "- Calculate per-spot QC metrics (total counts, n_genes, % mitochondrial)\n",
    "- Filter low-quality spots and genes\n",
    "- Generate QC plots (2x2 grid, violin, scatter)\n",
    "- Run sample-level QC classification (pass/fail with MAD outlier detection)\n",
    "- Produce an HTML QC report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import anndata as ad\n",
    "import sc_tools.qc as qc\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data\n",
    "\n",
    "We create a small dataset with a handful of mitochondrial genes (prefixed `MT-`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "n_spots, n_genes = 300, 150\n",
    "\n",
    "X = np.random.negative_binomial(3, 0.4, size=(n_spots, n_genes)).astype(float)\n",
    "\n",
    "# Add some mitochondrial genes\n",
    "mt_genes = [f\"MT-{g}\" for g in [\"CO1\", \"CO2\", \"ND1\", \"ND2\", \"ATP6\"]]\n",
    "other_genes = [f\"GENE_{i}\" for i in range(n_genes - len(mt_genes))]\n",
    "gene_names = mt_genes + other_genes\n",
    "\n",
    "adata = ad.AnnData(\n",
    "    X=X,\n",
    "    obs={\"sample\": [f\"sample_{i % 3}\" for i in range(n_spots)],\n",
    "         \"library_id\": [f\"lib_{i % 3}\" for i in range(n_spots)]},\n",
    ")\n",
    "adata.var_names = gene_names\n",
    "adata.obs_names = [f\"spot_{i}\" for i in range(n_spots)]\n",
    "\n",
    "print(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate QC metrics\n",
    "\n",
    "`calculate_qc_metrics` wraps scanpy's `pp.calculate_qc_metrics` and adds\n",
    "optional MT/HB percentage columns when patterns are provided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qc.calculate_qc_metrics(\n",
    "    adata,\n",
    "    mt_pattern=\"^MT-\",\n",
    "    hb_pattern=\"^HB\",\n",
    ")\n",
    "\n",
    "print(\"QC columns:\", [c for c in adata.obs.columns if \"count\" in c or \"gene\" in c or \"pct\" in c])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QC plots\n",
    "\n",
    "### 2x2 grid (pre-filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = qc.qc_2x2_grid(adata)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Violin metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = qc.qc_violin_metrics(adata, groupby=\"sample\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter: total counts vs n_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = qc.qc_scatter_counts_genes(adata)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter cells and genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Before: {adata.shape}\")\n",
    "\n",
    "qc.filter_cells(adata, min_genes=5)\n",
    "qc.filter_genes(adata, min_cells=3)\n",
    "\n",
    "print(f\"After:  {adata.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample-level QC\n",
    "\n",
    "Compute per-sample aggregate metrics and classify samples as pass/fail\n",
    "using absolute thresholds combined with adaptive MAD outlier detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_metrics = qc.compute_sample_metrics(adata, sample_col=\"sample\")\n",
    "print(sample_metrics[[\"n_spots\", \"median_counts\", \"median_genes\"]].to_string())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = qc.classify_samples(\n",
    "    sample_metrics,\n",
    "    modality=\"visium\",\n",
    "    mad_multiplier=3.0,\n",
    ")\n",
    "print(\"Pass/fail summary:\")\n",
    "print(results[[\"sample\", \"qc_pass\", \"fail_reasons\"]].to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HTML QC report\n",
    "\n",
    "In a real project, generate the HTML report and save to `figures/QC/qc_report.html`:\n",
    "\n",
    "```python\n",
    "qc.generate_qc_report(\n",
    "    adata,\n",
    "    sample_metrics=sample_metrics,\n",
    "    classification=results,\n",
    "    output_path=\"figures/QC/qc_report.html\",\n",
    ")\n",
    "```\n",
    "\n",
    "Or run the CLI script that produces all QC figures and the report in one step:\n",
    "\n",
    "```bash\n",
    "python scripts/run_qc_report.py \\\n",
    "    --input results/adata.raw.p1.h5ad \\\n",
    "    --output-dir figures/QC \\\n",
    "    --modality visium \\\n",
    "    --apply-filter\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sc_tools",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}