diff --git a/CHANGELOG.md b/CHANGELOG.md
index 49e4533..d06f530 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,12 @@ and this project adheres to [Semantic Versioning][].
 [keep a changelog]: https://keepachangelog.com/en/1.0.0/
 [semantic versioning]: https://semver.org/spec/v2.0.0.html
 
+## 2.2.0
+
+### Added
+- Bundled a Claude Code [Agent Skill](https://docs.claude.com/en/docs/claude-code/skills) that teaches coding agents to use decoupler correctly (method selection, calling conventions, output locations, priors, datasets, and benchmarking)
+- Added the `decoupler-install-skills` console script to install the skill into `~/.claude/skills/decoupler/` (use `--force` to refresh after upgrading, or `--print-path` to point `CLAUDE_SKILLS_PATH` at the bundled copy in place)
+
 ## 2.1.6
 
 ### Bugfixes
diff --git a/README.md b/README.md
index c3dc8e2..cef6b3a 100644
--- a/README.md
+++ b/README.md
@@ -79,6 +79,20 @@ mamba create -n=dcp conda-forge::decoupler-py
 pip install git+https://github.com/scverse/decoupler.git@main
 ```
 
+## Claude Code Skill
+
+decoupler ships an [Agent Skill](https://docs.claude.com/en/docs/claude-code/skills) for [Claude Code](https://claude.com/claude-code) that teaches the agent to use the library correctly. Install it once into your personal skills directory:
+
+```bash
+decoupler-install-skills
+```
+
+This copies the skill to `~/.claude/skills/decoupler/`, making it available in **every** project (re-run with `--force` after upgrading). Once installed, just ask Claude Code to do decoupler tasks ("find the most active TFs in this cell cluster", "score hallmark gene sets with ora", "build a consensus across methods") and the skill is consulted automatically. To avoid copying into your home directory, point Claude Code at the bundled copy instead:
+
+```bash
+export CLAUDE_SKILLS_PATH="$(decoupler-install-skills --print-path)"
+```
+
 ## Release notes
 
 See the [changelog][].
diff --git a/pyproject.toml b/pyproject.toml
index 4d543f3..27fb84a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ requires = [ "hatchling" ]
 
 [project]
 name = "decoupler"
-version = "2.1.6"
+version = "2.2.0"
 description = "Python package to perform enrichment analysis from omics data."
 readme = "README.md"
 license = { file = "LICENSE" }
@@ -38,6 +38,7 @@ dependencies = [
 urls.Documentation = "https://decoupler.readthedocs.io/"
 urls.Homepage = "https://github.com/scverse/decoupler"
 urls.Source = "https://github.com/scverse/decoupler"
+scripts.decoupler-install-skills = "decoupler._skills.install:main"
 
 [dependency-groups]
 dev = [
diff --git a/src/decoupler/_skills/__init__.py b/src/decoupler/_skills/__init__.py
new file mode 100644
index 0000000..34656be
--- /dev/null
+++ b/src/decoupler/_skills/__init__.py
@@ -0,0 +1,6 @@
+"""Bundled Claude Code Agent Skill and its installer.
+
+This subpackage keeps the decoupler Agent Skill (the ``data/`` directory) and the
+command that installs it (``install.py``) together and out of the main package
+namespace. See ``install.py`` for usage.
+"""
diff --git a/src/decoupler/_skills/data/SKILL.md b/src/decoupler/_skills/data/SKILL.md
new file mode 100644
index 0000000..4f038b9
--- /dev/null
+++ b/src/decoupler/_skills/data/SKILL.md
@@ -0,0 +1,85 @@
+---
+name: decoupler
+description: Use for any task involving the decoupler library — inferring biological activity/enrichment scores from omics data (bulk, single-cell, spatial). Triggers on estimating transcription factor (TF) activity, pathway activity, or gene-set enrichment from an AnnData/DataFrame; running ulm, mlm, ora, gsea, gsva, aucell, viper, zscore, waggr, mdt, udt; multi-method consensus; fetching prior knowledge (CollecTRI, DoRothEA, PROGENy, MSigDB hallmarks, OmniPath); pseudobulk and expression filtering; using bundled example datasets (pbmc3k, covid5k, erygast1k, hsctgfb, msvisium, knocktf, toy) and whether they are raw counts or normalized; ranking/plotting activity scores; or benchmarking methods against ground truth. This is a router skill — read the relevant file under references/ before writing decoupler code, because outputs land in different places and several defaults are non-obvious.
+---
+
+# decoupler
+
+`decoupler` estimates biological activities (transcription factors, pathways, gene
+sets) from omics data. You give it a **data matrix** (observations × features, e.g.
+cells × genes) and a **prior-knowledge network** (which features belong to which
+biological program), and it returns an **enrichment score per observation per
+program** — all methods sharing one unified interface. Part of scverse; works
+directly on `AnnData`, `pandas.DataFrame`, or raw matrices.
+
+This skill is a **router**. Each topic below has a detailed reference file with exact
+signatures and footguns. **Read the relevant reference before writing code** — do not
+rely on memory of the API, because results are written to different places depending
+on the input type and several methods have non-obvious defaults.
+
+## Two cross-cutting concepts (read these first if unsure)
+
+- **[The unified method call](references/calling-convention.md)** — every method is
+  called the same way: `dc.mt.<method>(data, net, tmin=5, ...)`. Covers the accepted
+  input types (AnnData / DataFrame / `[matrix, obs, var]` **list**), the shared
+  arguments (`tmin`, `raw`, `layer`, `empty`, `bsize`), and the big footguns: input
+  must be **normalized (e.g. log1p), not raw counts**, and `net` must be a
+  **long-format** `source`/`target`/`weight` table.
+
+- **[Where the results go](references/io-and-outputs.md)** — the return type depends
+  on the input type. **AnnData → written in place** into `.obsm["score_<method>"]`
+  (and `.obsm["padj_<method>"]` if the method tests), returns `None`. **DataFrame /
+  list → returns an `(es, pv)` tuple** (`pv` is `None` for non-testing methods). Read
+  scores back out of an AnnData with `dc.pp.get_obsm`. **This is the #1 thing agents
+  get wrong.**
+
+## Task → reference file
+
+| If the task is… | Read |
+|---|---|
+| **starting from scratch** — set up an end-to-end run (load data → get a net → score → rank → plot) | [references/getting-started.md](references/getting-started.md) |
+| how to **call a method**, what data/net formats are accepted, `tmin`/`raw`/`layer` semantics | [references/calling-convention.md](references/calling-convention.md) |
+| **where scores land** and how to read them back (`.obsm` vs `(es, pv)` tuple, `get_obsm`) | [references/io-and-outputs.md](references/io-and-outputs.md) |
+| **choosing and configuring a method** — ulm/mlm/ora/gsea/gsva/aucell/viper/zscore/waggr/mdt/udt + their kwargs | [references/methods.md](references/methods.md) |
+| getting **prior knowledge** — TF regulons (CollecTRI, DoRothEA), pathways (PROGENy), gene sets (MSigDB hallmarks), OmniPath, `.gmt`, organism translation | [references/priors.md](references/priors.md) |
+| **preprocessing** — pseudobulk from single cell, expression/sample filtering, layers, spatial kNN | [references/preprocessing.md](references/preprocessing.md) |
+| running **many methods at once** and combining them into a **consensus** | [references/multi-method.md](references/multi-method.md) |
+| **ranking** activity scores by group/ordering and **plotting** them (barplot, dotplot, volcano, network) or overlaying them on a **UMAP / spatial H&E** via scanpy | [references/ranking-and-plotting.md](references/ranking-and-plotting.md) |
+| **benchmarking** methods/nets against a perturbation ground truth | [references/benchmarking.md](references/benchmarking.md) |
+| using a **bundled example dataset** (`dc.ds.*`) — and whether its `.X` is raw counts or already normalized | [references/datasets.md](references/datasets.md) |
+
+## Quick orientation (submodule map)
+
+decoupler exposes everything through namespaced submodules — canonical import is
+`import decoupler as dc`:
+
+- `dc.mt` — **methods**: the enrichment/activity algorithms (`ulm`, `mlm`, `ora`,
+  `gsea`, `gsva`, `aucell`, `viper`, `zscore`, `waggr`, `mdt`, `udt`) plus
+  `dc.mt.decouple` (run several) and `dc.mt.consensus`. `dc.mt.show()` lists them.
+- `dc.op` — **prior knowledge** (OmniPath): `collectri`, `dorothea`, `progeny`,
+  `hallmark`, `resource`, `translate`.
+- `dc.pp` — **preprocessing**: `pseudobulk`, `filter_by_expr`, `filter_by_prop`,
+  `filter_samples`, `get_obsm`, `extract`, `swap_layer`, `knn`, `read_gmt`.
+- `dc.ds` — **datasets**: `toy` (used in every example: `adata, net = dc.ds.toy()`),
+  `pbmc3k`, `covid5k`, and other real examples. **They differ in normalization state**
+  (some raw counts, some already log-normalized) — see
+  [datasets.md](references/datasets.md) before scoring one.
+- `dc.tl` — **tools**: `rankby_group`, `rankby_obsm`, `rankby_order`.
+- `dc.pl` — **plotting**: `barplot`, `dotplot`, `volcano`, `network`, `obsm`, etc.
+- `dc.bm` — **benchmarking**: `benchmark`, `metric`, `pl`.
+
+## Conventions used throughout
+
+- **Data** = observations × features (cells/samples × genes), **normalized** (log1p),
+  not raw counts. **Net** = long-format DataFrame: `source` (program), `target`
+  (feature), optional signed `weight`.
+- Output score matrices are observations × sources. AnnData results live in `.obsm`
+  under `score_<method>` / `padj_<method>`.
+- `tmin` (default 5) silently drops sources with fewer than `tmin` targets present in
+  the data — lower it for small/toy nets (examples use `tmin=3`).
+- `verbose=True` on any method prints the pruning/run log — use it to see how many
+  sources survived `tmin`.
+- **Plotting**: use `dc.pl.*` for summary plots and scanpy `sc.pl.umap` / `sc.pl.spatial`
+  (on a `dc.pp.get_obsm` score AnnData) for embedding/tissue overlays. Never reimplement
+  a plot from matplotlib primitives — see
+  [ranking-and-plotting.md](references/ranking-and-plotting.md).
diff --git a/src/decoupler/_skills/data/references/benchmarking.md b/src/decoupler/_skills/data/references/benchmarking.md
new file mode 100644
index 0000000..4c399b6
--- /dev/null
+++ b/src/decoupler/_skills/data/references/benchmarking.md
@@ -0,0 +1,59 @@
+# Benchmarking methods and nets
+
+`dc.bm` evaluates how well a method + net recover a known ground truth — typically a
+perturbation dataset where you know which program was up/down in each sample. Use it
+to **choose a method** or **validate a net**, not for routine scoring.
+
+## `dc.bm.benchmark`
+
+```python
+dc.bm.benchmark(
+    adata,                 # AnnData with a perturbation ground truth in .obs
+    net,                   # a net, or a dict of {name: net} to compare several
+    metrics=None,          # e.g. ["auroc", "auprc"] (None = defaults)
+    groupby=None,          # obs column(s) to stratify the evaluation
+    runby="expr",          # evaluate per-observation ("expr") or per-source
+    sfilt=False,           # filter sources to those with a ground-truth label
+    thr=0.1, emin=5,
+    kws_decouple=None,     # kwargs forwarded to dc.mt.decouple (methods, args, tmin, ...)
+    verbose=False,
+)
+```
+
+- Internally runs `dc.mt.decouple`, so pass method selection and per-method args via
+  `kws_decouple` (e.g. `kws_decouple={"methods": ["ulm", "mlm"], "tmin": 3}`).
+- Pass a **dict of nets** to benchmark several priors head to head.
+- The ground truth (which source is perturbed, and its sign) must be encoded in
+  `adata.obs` — see the decoupler benchmarking tutorial for the expected columns.
+
+## Metrics and plots
+
+- `dc.bm.metric` — the scoring functions: `auc`, `fscore`, `qrank`, `hmean`,
+  `dict_metric`.
+- `dc.bm.pl` — plots of benchmark results: `auc`, `bar`, `fscore`, `qrank`, `summary`.
+
+```python
+res = dc.bm.benchmark(adata, {"collectri": net_a, "dorothea": net_b},
+                      kws_decouple={"methods": ["ulm", "mlm"]})
+dc.bm.pl.summary(res)   # compare methods × nets
+```
+
+## When to use
+
+- Deciding **which method** to trust for your data/organism.
+- Comparing **candidate nets** (e.g. CollecTRI vs DoRothEA levels).
+- Sanity-checking a **custom net** against a perturbation benchmark
+  (`dc.ds.knocktf` and similar provide ground-truth datasets).
+
+## Footguns
+
+- Benchmarking is only meaningful with a **real ground truth** — a toy net won't tell
+  you anything.
+- It re-runs scoring internally; keep `kws_decouple` consistent with how you intend to
+  score for real, or the comparison won't transfer.
+
+## Related references
+
+[multi-method.md](multi-method.md) (`decouple`, which `benchmark` wraps),
+[methods.md](methods.md) (candidate methods),
+[priors.md](priors.md) (candidate nets).
diff --git a/src/decoupler/_skills/data/references/calling-convention.md b/src/decoupler/_skills/data/references/calling-convention.md
new file mode 100644
index 0000000..c435986
--- /dev/null
+++ b/src/decoupler/_skills/data/references/calling-convention.md
@@ -0,0 +1,91 @@
+# The unified method call
+
+Every enrichment method in `dc.mt` is a callable object with the **same signature**.
+Learn it once and it applies to `ulm`, `mlm`, `ora`, `gsea`, `gsva`, `aucell`,
+`viper`, `zscore`, `waggr`, `mdt`, `udt` (and `dc.mt.decouple`).
+
+```python
+import decoupler as dc
+
+adata, net = dc.ds.toy()          # AnnData (obs × var) + long-format net
+dc.mt.ulm(adata, net, tmin=3)     # writes adata.obsm["score_ulm"], ["padj_ulm"]
+```
+
+## Shared signature
+
+```python
+dc.mt.<method>(
+    data,              # AnnData | pandas.DataFrame | [matrix, obs_names, var_names]
+    net,               # long-format DataFrame: source, target, [weight]
+    tmin=5,            # drop sources with < tmin targets present in the data
+    raw=False,         # use adata.raw instead of adata.X
+    layer=None,        # use adata.layers[layer] instead of adata.X (passed via kwargs)
+    empty=True,        # drop all-zero observations/features before scoring
+    bsize=250_000,     # batch size (rows) for sparse / backed matrices
+    verbose=False,
+    **kwargs,          # method-specific args — see references/methods.md
+)
+```
+
+`layer=` is accepted even though it is not in the explicit signature — it is forwarded
+through `**kwargs` to the runner. Use `layer=` **or** `raw=`, not both.
+
+## Accepted input types (`data`)
+
+| Form | What it is | Returns |
+|---|---|---|
+| `AnnData` | scored on `.X` (or `.raw`/`layer`) | **`None`** — results written to `.obsm` in place |
+| `pandas.DataFrame` | observations (rows) × features (cols) | `(es, pv)` tuple of DataFrames |
+| `list` `[matrix, obs_names, var_names]` | numpy/sparse matrix + name arrays | `(es, pv)` tuple of DataFrames |
+
+See [io-and-outputs.md](io-and-outputs.md) for the return details.
+
+### Footgun — the matrix form must be a `list`, not a `tuple`
+
+```python
+# WRONG — a tuple raises "mat must be a list of [matrix, samples, features]"
+es, pv = dc.mt.ulm((X, obs, var), net)
+
+# RIGHT — wrap it in a list
+es, pv = dc.mt.ulm([X, obs, var], net)
+```
+
+## The three footguns that silently give wrong results
+
+1. **Input must be normalized, not raw counts.** decoupler expects
+   log-normalized expression — the same matrix you would use for differential
+   expression. **When normalizing, use `sc.pp.normalize_total(adata, target_sum=1e4)`**
+   followed by `sc.pp.log1p(adata)`. Passing raw counts runs without error but produces
+   meaningless scores. If your normalized data is in a layer, pass
+   `layer="lognorm"`; if in `.raw`, pass `raw=True`. **Check before normalizing:** if
+   `.X` is integer-valued with a large max it is raw counts; if it is small floats
+   (often with `adata.uns["log1p"]`) it is already normalized — do not normalize twice.
+   The bundled datasets differ; see [datasets.md](datasets.md) (e.g. `pbmc3k` is already
+   log-normalized, `covid5k` is raw counts).
+
+2. **`net` must be long-format** with columns `source`, `target`, and (for weighted
+   methods) `weight` — one row per (program, feature) edge. A wide/matrix net or wrong
+   column names will fail or misbehave. See [priors.md](priors.md).
+
+3. **`tmin` silently prunes.** Sources with fewer than `tmin` (default 5) targets
+   present in `data` are dropped *before* scoring — they simply won't appear in the
+   output. On small or toy nets this can drop everything; the examples use `tmin=3`.
+   Run with `verbose=True` to see `<method> - X sources with < tmin targets` and how
+   many survived.
+
+## Listing and introspecting methods
+
+```python
+dc.mt.show()          # DataFrame: name, desc, stype, weight, test, limits, reference
+```
+
+- `weight` — whether the method uses signed edge weights (`True` for ulm/mlm/viper/…).
+- `test` — whether it produces p-values (`padj`). `False` for aucell/gsva/mdt/udt.
+- `stype` / `limits` — value type and output range (e.g. aucell `(0, 1)`).
+
+## Related references
+
+[io-and-outputs.md](io-and-outputs.md) (return types, reading `.obsm` back),
+[methods.md](methods.md) (per-method kwargs and how to choose one),
+[priors.md](priors.md) (building/fetching `net`),
+[preprocessing.md](preprocessing.md) (getting a normalized matrix / layers).
diff --git a/src/decoupler/_skills/data/references/datasets.md b/src/decoupler/_skills/data/references/datasets.md
new file mode 100644
index 0000000..5c9c05d
--- /dev/null
+++ b/src/decoupler/_skills/data/references/datasets.md
@@ -0,0 +1,88 @@
+# Bundled example datasets (`dc.ds`)
+
+decoupler ships loaders for a toy dataset and several real ones. **They are not all in
+the same state** — some `.X` are raw counts, some are already log-normalized, one is a
+log2FC contrast statistic. Since methods expect **normalized** input (see
+[calling-convention.md](calling-convention.md)), you must know which is which or you
+will either skip a needed normalization or re-normalize already-normalized data.
+
+## Rule: inspect `.X` before assuming counts
+
+```python
+import numpy as np, scipy.sparse as sps
+X = adata.X.toarray() if sps.issparse(adata.X) else np.asarray(adata.X)
+X.max(), np.allclose(X[:50], np.round(X[:50])), "log1p" in adata.uns
+```
+
+- **Integer-valued and large max (hundreds–millions)** → raw counts → normalize first.
+- **Float, small max (roughly ≤ tens), often `adata.uns["log1p"]` present** → already
+  log-normalized → **score directly, do not normalize again**.
+- **Signed values (negatives)** → a contrast statistic (e.g. log2FC), not expression →
+  score directly.
+
+## The datasets
+
+| Loader | What it is | `.X` state | Before scoring | Organism |
+|---|---|---|---|---|
+| `dc.ds.toy()` | synthetic AnnData **+ matching net** | synthetic continuous | score directly (`tmin=3`) | — |
+| `dc.ds.toy_bench()` | toy + benchmark ground truth | synthetic continuous | feed to `dc.bm` | — |
+| `dc.ds.pbmc3k()` | 10x PBMCs (~2.6k cells) | **log-normalized** (max ~6) | **score directly** | human |
+| `dc.ds.erygast1k()` | mouse erythroid gastrulation (~800 cells) | **log-normalized** (`uns['log1p']`) | score directly | **mouse** |
+| `dc.ds.covid5k()` | COVID-19 PBMCs (~5k cells) | **raw counts** (max ~1e4) | **normalize + log1p** | human |
+| `dc.ds.hsctgfb()` | bulk RNA-seq, 6 HSC samples | **raw counts** | **filter + normalize** (bulk) | human |
+| `dc.ds.msvisium()` | Visium MS brain slide (spatial) | **raw counts** | **normalize** | human |
+| `dc.ds.knocktf()` | KnockTF TF-perturbation contrasts | **log2FC** (−20…20) | score directly (benchmark) | human |
+
+### Ready to score directly (do NOT normalize)
+
+```python
+adata = dc.ds.pbmc3k()            # already log-normalized
+net = dc.op.collectri(organism="human")
+dc.mt.ulm(adata, net)             # -> adata.obsm["score_ulm"], ["padj_ulm"]
+```
+
+- `erygast1k` is the same, but **mouse** — fetch a mouse net (`organism="mouse"`) or
+  translate a human one (`dc.op.translate`, see [priors.md](priors.md)). Its `obs` has
+  `celltype`, `stage`.
+- `pbmc3k.obs` has `celltype` and `leiden`; `obsm` has PCA/UMAP/tSNE already.
+
+### Raw counts → normalize first
+
+```python
+import scanpy as sc
+adata = dc.ds.covid5k()                    # raw integer counts
+sc.pp.normalize_total(adata, target_sum=1e4)
+sc.pp.log1p(adata)
+dc.mt.ulm(adata, dc.op.collectri())
+```
+
+- `covid5k.obs`: `individual`, `sex`, `disease`, `celltype`.
+- `hsctgfb` is **bulk** (6 samples, `obs`: `condition`, `sample_id`) — filter lowly
+  expressed genes with `dc.pp.filter_by_expr` then normalize; a natural fit for a
+  bulk/DE + activity workflow ([preprocessing.md](preprocessing.md)).
+- `msvisium` is **spatial** — normalize, then optionally build a neighbor graph with
+  `dc.pp.knn(adata, key="spatial")`. It carries `obsm["spatial"]`, `uns["spatial"]`
+  (H&E images), and `obs["niches"]`. To plot activity on the H&E, use `dc.pp.get_obsm` +
+  `sc.pl.spatial` (see [ranking-and-plotting.md](ranking-and-plotting.md)).
+
+### Benchmark inputs (not for routine scoring)
+
+- `dc.ds.knocktf()` — `.X` is already log2FCs; `obs["source"]` is the perturbed TF and
+  `obs["logFC"]` its effect. It is the ground truth for `dc.bm.benchmark`, not a matrix
+  you normalize. `thr_fc=-1` (default) keeps clearly down-regulated perturbations.
+- `dc.ds.toy_bench()` — toy data with ground truth in `obs["source"]` / `class` /
+  `type_p` for exercising the benchmark pipeline ([benchmarking.md](benchmarking.md)).
+
+## Gene identifiers
+
+Real datasets use gene **symbols** in `var_names` (matching the OmniPath priors). If you
+bring data with Ensembl IDs, convert with `dc.ds.ensmbl_to_symbol(...)` so `net`
+`target`s line up — mismatched namespaces get silently pruned by `tmin`
+(see the identifier footgun in [priors.md](priors.md)).
+
+## Related references
+
+[calling-convention.md](calling-convention.md) (normalized-input requirement),
+[preprocessing.md](preprocessing.md) (normalize / filter / pseudobulk / spatial kNN),
+[priors.md](priors.md) (matching organism and gene identifiers),
+[benchmarking.md](benchmarking.md) (`knocktf`, `toy_bench`).
diff --git a/src/decoupler/_skills/data/references/getting-started.md b/src/decoupler/_skills/data/references/getting-started.md
new file mode 100644
index 0000000..a2d72b9
--- /dev/null
+++ b/src/decoupler/_skills/data/references/getting-started.md
@@ -0,0 +1,82 @@
+# Getting started (end-to-end)
+
+A complete decoupler run has four steps: **get a normalized data matrix**, **get a
+prior net**, **score**, then **inspect/rank/plot**. Below are the two most common
+end-to-end shapes.
+
+## Minimal, self-contained (toy data)
+
+```python
+import decoupler as dc
+
+# 1. Data (obs × var, normalized) + 2. a prior net (source/target/weight)
+adata, net = dc.ds.toy()
+
+# 3. Score — results written into adata.obsm in place
+dc.mt.ulm(adata, net, tmin=3)
+
+# 4. Inspect: scores + adjusted p-values, both obs × sources
+adata.obsm["score_ulm"]
+adata.obsm["padj_ulm"]
+```
+
+`dc.ds.toy()` returns a small AnnData and a matching weighted net — use it for any
+quick check or when reproducing behavior.
+
+## Single-cell TF activity (realistic shape)
+
+```python
+import scanpy as sc
+import decoupler as dc
+
+# adata: raw counts in a single-cell AnnData
+sc.pp.normalize_total(adata, target_sum=1e4)   # when normalizing, use target_sum=1e4
+sc.pp.log1p(adata)                      # decoupler needs normalized input, not counts
+# NOTE: only normalize raw counts. Data that is already log-normalized (e.g.
+# dc.ds.pbmc3k()) must NOT be normalized again — see references/datasets.md.
+
+net = dc.op.collectri(organism="human")  # signed TF regulons
+
+dc.mt.ulm(adata, net)                    # -> adata.obsm["score_ulm"], ["padj_ulm"]
+
+# Turn scores into an AnnData for downstream scanpy/plotting
+acts = dc.pp.get_obsm(adata, key="score_ulm")   # obs × TFs
+
+# Which TFs are most active per cell type?
+df = dc.tl.rankby_group(acts, groupby="cell_type")
+# barplot/dotplot take DataFrames (not `acts`) — see ranking-and-plotting.md
+```
+
+## Pathway activity or gene-set enrichment
+
+Swap the net (and, for unweighted gene sets, the method):
+
+```python
+# Pathways (weighted) — ulm is fine
+prog = dc.op.progeny(organism="human", top=500)
+dc.mt.ulm(adata, prog)                   # adata.obsm["score_ulm"] = pathway activities
+
+# Hallmark gene sets (unweighted) — use a gene-set method
+hall = dc.op.hallmark(organism="human")
+dc.mt.ora(adata, hall)                   # adata.obsm["score_ora"], ["padj_ora"]
+```
+
+## Bulk / pseudobulk
+
+For bulk RNA-seq or pseudobulk from single cell, aggregate and filter first (see
+[preprocessing.md](preprocessing.md)), then score exactly the same way:
+
+```python
+pdata = dc.pp.pseudobulk(adata, sample_col="sample", groups_col="cell_type")
+dc.pp.filter_by_expr(pdata)              # drop lowly-expressed genes
+# ... normalize pdata (e.g. CPM + log1p) ...
+dc.mt.ulm(pdata, net)
+```
+
+## Where to go next
+
+- Method choice and kwargs → [methods.md](methods.md)
+- Getting/formatting the net → [priors.md](priors.md)
+- Return types and reading scores back → [io-and-outputs.md](io-and-outputs.md)
+- Running many methods + consensus → [multi-method.md](multi-method.md)
+- Ranking + plotting → [ranking-and-plotting.md](ranking-and-plotting.md)
diff --git a/src/decoupler/_skills/data/references/io-and-outputs.md b/src/decoupler/_skills/data/references/io-and-outputs.md
new file mode 100644
index 0000000..e9e6cc8
--- /dev/null
+++ b/src/decoupler/_skills/data/references/io-and-outputs.md
@@ -0,0 +1,80 @@
+# Where the results go (input type decides the return)
+
+The **return type of every `dc.mt` method depends on the input type**. This is the
+single most common source of confusion — an AnnData call looks like it returned
+nothing, but it wrote the results into the object.
+
+## AnnData in → results written in place, returns `None`
+
+```python
+adata, net = dc.ds.toy()
+out = dc.mt.ulm(adata, net, tmin=3)
+assert out is None                       # nothing is returned
+
+adata.obsm["score_ulm"]                  # DataFrame: obs × sources (activity scores)
+adata.obsm["padj_ulm"]                   # DataFrame: obs × sources (FDR-adj p-values)
+```
+
+- Results are stored under `.obsm` with keys `score_<method>` and (only if the method
+  tests) `padj_<method>`. Method names: `score_ulm`, `padj_ulm`, `score_ora`, etc.
+- **Do not** write `adata = dc.mt.ulm(adata, net)` — that overwrites `adata` with
+  `None` and loses your data.
+- Non-testing methods (`aucell`, `gsva`, `mdt`, `udt`) write only `score_<method>`,
+  no `padj_<method>`.
+- Edge case: if `empty=True` dropped observations, decoupler returns a repaired copy
+  of the AnnData instead of `None`; capture it if you pass data with empty rows.
+
+### Reading scores back out as an AnnData
+
+`dc.pp.get_obsm` pulls an `.obsm` score matrix into a new AnnData whose `.X` is the
+scores and whose `var` are the sources. It **preserves `obs`, `uns`, and `obsm`** from
+the original, so the result is ready for scanpy — you can `sc.pl.umap` / `sc.pl.spatial`
+it directly, coloring by a source (`var_name`) or an `obs` column. This is the intended
+way to visualize scores; see
+[ranking-and-plotting.md](ranking-and-plotting.md).
+
+```python
+acts = dc.pp.get_obsm(adata, key="score_ulm")   # AnnData: obs × sources
+acts.var_names                                   # the programs (e.g. TFs)
+# acts.obs / acts.uns / acts.obsm are carried over from adata (e.g. spatial images)
+```
+
+## DataFrame or list in → returns an `(es, pv)` tuple
+
+```python
+import pandas as pd
+df = pd.DataFrame(adata.X, index=adata.obs_names, columns=adata.var_names)
+
+es, pv = dc.mt.ulm(df, net, tmin=3)   # es: obs × sources scores; pv: obs × sources padj
+```
+
+- `es` (enrichment scores) and `pv` (adjusted p-values) are both `obs × sources`
+  DataFrames.
+- For **non-testing** methods, `pv is None`:
+  ```python
+  es, pv = dc.mt.aucell(df, net, tmin=3)
+  assert pv is None
+  ```
+- The `list` form `[matrix, obs_names, var_names]` returns the same `(es, pv)` tuple.
+
+## Multi-method returns
+
+- `dc.mt.decouple(...)` returns a **dict** keyed `score_<m>` / `padj_<m>` for DataFrame
+  input, or writes those keys into `.obsm` for AnnData input (returns `None`).
+- `dc.mt.consensus(...)` reads all `score_*` entries and returns/writes a single
+  consensus score. See [multi-method.md](multi-method.md).
+
+## Quick decision table
+
+| You passed… | You get back… | Where scores are |
+|---|---|---|
+| `AnnData` | `None` | `adata.obsm["score_<m>"]`, `["padj_<m>"]` |
+| `DataFrame` / `list` | `(es, pv)` | the returned DataFrames (`pv=None` if no test) |
+| `decouple(AnnData)` | `None` | multiple `adata.obsm[...]` keys |
+| `decouple(DataFrame)` | `dict` | `dict["score_<m>"]`, `dict["padj_<m>"]` |
+
+## Related references
+
+[calling-convention.md](calling-convention.md) (input types and shared args),
+[ranking-and-plotting.md](ranking-and-plotting.md) (using the scores downstream),
+[multi-method.md](multi-method.md) (`decouple` / `consensus` outputs).
diff --git a/src/decoupler/_skills/data/references/methods.md b/src/decoupler/_skills/data/references/methods.md
new file mode 100644
index 0000000..0cba852
--- /dev/null
+++ b/src/decoupler/_skills/data/references/methods.md
@@ -0,0 +1,67 @@
+# Choosing and configuring a method
+
+All methods share the call in [calling-convention.md](calling-convention.md). This
+file covers **which method to pick** and each method's **extra kwargs**. List them
+live with `dc.mt.show()`.
+
+## The methods at a glance
+
+| Method | Kind | Uses weights? | p-values? | Output range | Typical use |
+|---|---|:--:|:--:|---|---|
+| `ulm` | Univariate Linear Model | yes (signed) | yes | (-inf, inf) | **Default** for TF/pathway activity; fast, robust |
+| `mlm` | Multivariate Linear Model | yes (signed) | yes | (-inf, inf) | Accounts for correlated regulons; needs full-rank net |
+| `viper` | VIPER regulon analysis | yes (signed) | yes | (-inf, inf) | TF activity with pleiotropy/target-overlap correction |
+| `waggr` | Weighted Aggregate | yes (signed) | yes | (-inf, inf) | wmean/wsum/median aggregate; permutation p-values |
+| `zscore` | Z-score | yes (signed) | yes | (-inf, inf) | Simple standardized weighted mean |
+| `gsea` | Gene Set Enrichment Analysis | no | yes | (-inf, inf) | Ranked enrichment of gene sets |
+| `ora` | Over-Representation Analysis | no | yes | (-inf, inf) | Fisher test on top/bottom features |
+| `aucell` | AUCell | no | no | (0, 1) | Rank-based gene-set activity per cell |
+| `gsva` | Gene Set Variation Analysis | no | no | (-1, 1) | Sample-relative gene-set variation |
+| `mdt` | Multivariate Decision Tree | yes | no | (0, 1) | Random-forest feature importance as activity |
+| `udt` | Univariate Decision Tree | yes | no | (0, 1) | Single-tree importance as activity |
+
+Rules of thumb:
+- **Weighted net (has signed `weight`, e.g. CollecTRI/DoRothEA/PROGENy)** → prefer
+  `ulm` (default), `mlm`, `viper`, or `waggr`.
+- **Unweighted gene sets (e.g. MSigDB hallmarks, `.gmt`)** → `ora`, `gsea`, `aucell`,
+  or `gsva`.
+- **Unsure / want robustness** → run several and take a consensus
+  ([multi-method.md](multi-method.md)).
+
+## Method-specific kwargs
+
+Pass these as extra keyword args to the call, e.g. `dc.mt.ora(adata, net, n_up=300)`.
+
+- **`ulm`, `mlm`** — `tval=True`: return the t-value (default) vs the raw coefficient.
+- **`ora`** — `n_up=None` (n top features as the "up" set; default derives from data),
+  `n_bm=0` (n bottom features), `n_bg=20000` (background size for the Fisher test),
+  `ha_corr=0.5` (Haldane–Anscombe continuity correction).
+- **`gsea`** — `times=1000` (permutations), `seed=42`. **Reproducibility:** p-values
+  come from permutations; fix `seed` (and keep `times` constant) for stable results.
+- **`gsva`** — `kcdf="gaussian"`, `maxdiff=True`, `absrnk=False`, `tau=1`.
+- **`aucell`** — `n_up=None`: size of the ranking cutoff (top-N genes per cell;
+  default derives from the data).
+- **`viper`** — `pleiotropy=True`, `reg_sign=0.05`, `n_targets=10`, `penalty=20`
+  (pleiotropy correction for shared targets between regulons).
+- **`zscore`** — `flavor="RoKAI"` (alternative standardization flavor).
+- **`waggr`** — `fun="wmean"` (or `"wsum"`, `"median"`), `times=1000`, `seed=42`
+  (permutation p-values; fix `seed` for reproducibility).
+- **`mdt`, `udt`** — no extra kwargs.
+
+## Notes / footguns
+
+- **`ulm` is the recommended default** for TF and pathway activity in most workflows.
+- **`mlm`** solves a single multivariate regression; if the net is rank-deficient
+  (highly collinear regulons) it can be unstable — prefer `ulm`/`viper` then.
+- **`mlm` p-values are not FDR-adjusted** the way other testing methods are (it is
+  the one method the runner skips BH correction for); interpret `padj_mlm` accordingly.
+- **Permutation methods (`gsea`, `waggr`)** are stochastic — set `seed` and hold
+  `times` fixed, or scores/p-values will drift between runs.
+- **Non-testing methods (`aucell`, `gsva`, `mdt`, `udt`)** produce no `padj_<method>`.
+
+## Related references
+
+[calling-convention.md](calling-convention.md) (shared args, input formats),
+[multi-method.md](multi-method.md) (`decouple` / `consensus`),
+[priors.md](priors.md) (weighted vs gene-set nets),
+[io-and-outputs.md](io-and-outputs.md) (where scores land).
diff --git a/src/decoupler/_skills/data/references/multi-method.md b/src/decoupler/_skills/data/references/multi-method.md
new file mode 100644
index 0000000..c7c007f
--- /dev/null
+++ b/src/decoupler/_skills/data/references/multi-method.md
@@ -0,0 +1,72 @@
+# Running many methods + consensus
+
+When you don't want to commit to one method, run several and combine them. Two tools:
+`dc.mt.decouple` (run N methods) and `dc.mt.consensus` (aggregate their scores).
+
+## `dc.mt.decouple` — run several methods at once
+
+```python
+dc.mt.decouple(
+    data, net,
+    methods="all",     # or a list, e.g. ["ulm", "ora", "aucell"]
+    args=None,         # per-method kwargs: {"ora": {"n_up": 300}, "gsea": {"seed": 0}}
+    cons=False,        # also compute a consensus score
+    **kwargs,          # shared kwargs forwarded to every method (e.g. tmin=3)
+)
+```
+
+- **AnnData input** → writes `score_<m>` / `padj_<m>` for each method into `.obsm`,
+  returns `None`.
+- **DataFrame input** → returns a **dict** keyed `score_<m>` / `padj_<m>`.
+- `methods="all"` runs every applicable method; pass a list to restrict.
+- `args` gives per-method overrides; top-level `**kwargs` (like `tmin`) apply to all.
+
+```python
+adata, net = dc.ds.toy()
+dc.mt.decouple(adata, net, methods=["ulm", "ora"], tmin=3)
+adata.obsm["score_ulm"], adata.obsm["score_ora"]
+```
+
+## `dc.mt.consensus` — combine method scores
+
+```python
+dc.mt.consensus(result, verbose=False)
+```
+
+Computes a single consensus activity by aggregating the per-method scores (a signed
+z-score mean across methods), with its own p-values.
+
+- Pass an **AnnData** that already has `score_*` in `.obsm` → writes
+  `score_consensus` / `padj_consensus`, returns `None`.
+- Pass the **dict** returned by `decouple` (DataFrame input) → returns
+  `(es, pv)`.
+- Shortcut: `dc.mt.decouple(..., cons=True)` runs the methods and adds the consensus
+  in one call.
+
+```python
+# One-shot: run several methods and get the consensus
+dc.mt.decouple(adata, net, methods=["ulm", "mlm", "ora"], cons=True, tmin=3)
+adata.obsm["score_consensus"]
+```
+
+## When to use this
+
+- **Robustness / reduce method bias** — consensus across methods is a common default
+  for TF activity when you don't want to defend a single choice.
+- **Comparing methods** — `decouple` gives you every method's scores side by side in
+  one object.
+
+## Footguns
+
+- Consensus needs **at least two** method scores present; running it after a single
+  method is meaningless.
+- Mixing methods with very different ranges is fine (consensus standardizes), but keep
+  the **same net** across methods so sources line up.
+- Per-method reproducibility still applies — set `seed` for `gsea`/`waggr` via
+  `args={"gsea": {"seed": 0}}`.
+
+## Related references
+
+[methods.md](methods.md) (per-method kwargs, which to include),
+[io-and-outputs.md](io-and-outputs.md) (dict vs `.obsm` returns),
+[benchmarking.md](benchmarking.md) (deciding which methods win on your data).
diff --git a/src/decoupler/_skills/data/references/preprocessing.md b/src/decoupler/_skills/data/references/preprocessing.md
new file mode 100644
index 0000000..444820f
--- /dev/null
+++ b/src/decoupler/_skills/data/references/preprocessing.md
@@ -0,0 +1,66 @@
+# Preprocessing (`dc.pp`)
+
+Helpers to get your data into the shape decoupler expects: pseudobulk aggregation,
+expression/sample filtering, layer handling, and spatial neighborhoods. Reminder:
+methods need **normalized** input, not raw counts (see
+[calling-convention.md](calling-convention.md)).
+
+## Pseudobulk (single cell → sample-level)
+
+```python
+pdata = dc.pp.pseudobulk(
+    adata,
+    sample_col="sample",        # obs column identifying each sample/replicate
+    groups_col="cell_type",     # obs column to aggregate within (or None)
+    mode="sum",                 # "sum" (default), "mean", or a callable/dict
+    layer=None, raw=False,      # where to read counts from
+)
+```
+
+Returns a new AnnData of pseudobulk profiles (one per sample × group). Aggregate
+**raw counts** here, then filter and normalize before scoring.
+
+## Filtering
+
+- `dc.pp.filter_by_expr(adata, group=None, min_count=10, min_total_count=15, ...)` —
+  edgeR-style filter of lowly-expressed genes on count data (for bulk/pseudobulk).
+  In-place by default (`inplace=True`); pass `inplace=False` to get a mask.
+- `dc.pp.filter_by_prop(adata, min_prop=0.2, min_smpls=2)` — keep genes expressed in a
+  minimum proportion of cells across a minimum number of samples.
+- `dc.pp.filter_samples(adata, min_cells=10, min_counts=1000)` — drop low-quality
+  pseudobulk samples.
+
+Typical bulk/pseudobulk order: `pseudobulk` → `filter_by_expr` → normalize (e.g.
+CPM + log1p) → `dc.mt.<method>`.
+
+## Layers and reading scores back
+
+- `dc.pp.get_obsm(adata, key="score_ulm")` — pull an `.obsm` score matrix into a new
+  AnnData (`.X` = scores, `var` = sources). The bridge to scanpy/plotting.
+- `dc.pp.swap_layer(adata, key, X_key="X", inplace=False)` — move a layer into `.X`
+  (or vice versa) so you can control which matrix a method reads.
+- `dc.pp.extract(data, layer=None, raw=False)` — the low-level extractor methods use
+  internally; rarely needed directly, but useful to see exactly what matrix/obs/var a
+  method will consume.
+
+## Spatial
+
+- `dc.pp.knn(adata, key="spatial", bw=100, max_nn=100)` — build a spatial
+  neighbor graph (used before spatially-smoothed analyses).
+- `dc.pp.bin_order(adata, order, nbins=100)` — bin observations along a continuous
+  ordering (e.g. pseudotime) for trajectory-style summaries.
+
+## Footguns
+
+- **Pseudobulk on counts, score on normalized.** `pseudobulk` should aggregate raw
+  counts; normalize the resulting `pdata` before `dc.mt.*`.
+- **`filter_by_expr` expects counts**, not log-normalized values — run it before
+  normalization.
+- **`inplace=True` is the default** for the filters — they mutate `adata` and return
+  `None`. Pass `inplace=False` if you want a boolean mask instead.
+
+## Related references
+
+[getting-started.md](getting-started.md) (where preprocessing fits in the flow),
+[calling-convention.md](calling-convention.md) (normalization requirement, `layer`/`raw`),
+[io-and-outputs.md](io-and-outputs.md) (`get_obsm`).
diff --git a/src/decoupler/_skills/data/references/priors.md b/src/decoupler/_skills/data/references/priors.md
new file mode 100644
index 0000000..8e66f23
--- /dev/null
+++ b/src/decoupler/_skills/data/references/priors.md
@@ -0,0 +1,79 @@
+# Prior knowledge (the `net`)
+
+Every method needs a **`net`**: a long-format `pandas.DataFrame` describing which
+features (genes) belong to which biological program.
+
+## Net format
+
+| column | meaning | required |
+|---|---|---|
+| `source` | the program (TF, pathway, gene set) | yes |
+| `target` | the feature (gene) | yes |
+| `weight` | signed importance of the edge | only for weighted methods |
+
+One row per edge. Weighted methods (`ulm`, `mlm`, `viper`, `zscore`, `waggr`, `mdt`,
+`udt`) use `weight` (sign encodes activation/repression); gene-set methods (`ora`,
+`gsea`, `aucell`, `gsva`) ignore it.
+
+```python
+import decoupler as dc
+net = dc.op.collectri(organism="human")
+net.head()      # columns: source, target, weight
+```
+
+## Fetching curated priors from `dc.op` (OmniPath)
+
+All default to `organism="human"` and `license="academic"`.
+
+- **TF regulons**
+  - `dc.op.collectri(organism="human", remove_complexes=False)` — CollecTRI TF→gene
+    regulons with signed weights. Recommended default for TF activity.
+  - `dc.op.dorothea(organism="human", levels=None)` — DoRothEA regulons; `levels`
+    selects confidence levels (e.g. `["A", "B", "C"]`).
+- **Pathways**
+  - `dc.op.progeny(organism="human", top=inf, thr_padj=0.05)` — PROGENy responsive
+    genes; `top` limits genes per pathway (e.g. `top=500` is common).
+- **Gene sets**
+  - `dc.op.hallmark(organism="human")` — MSigDB hallmark gene sets (no weights → use
+    with `ora`/`gsea`/`aucell`/`gsva`).
+- **Anything else on OmniPath**
+  - `dc.op.resource(name, organism="human")` — fetch an arbitrary resource by name.
+  - `dc.op.show_resources()` — list available resource names.
+
+## Loading your own gene sets from a `.gmt`
+
+```python
+net = dc.pp.read_gmt("my_sets.gmt")   # -> long-format source/target (no weights)
+```
+
+Use with a gene-set method (`ora`, `gsea`, `aucell`, `gsva`).
+
+## Cross-species: translate a net to another organism
+
+`dc.op` priors accept `organism=` directly; to translate an existing human net to
+another species by orthology use `dc.op.translate`:
+
+```python
+mouse_net = dc.op.translate(net, target_organism="mouse", min_evidence=3)
+dc.op.show_organisms()   # list supported organisms
+```
+
+- Translates the `source`/`target` gene symbols (control columns via `columns=`).
+- `min_evidence` and `one_to_many` control ortholog stringency.
+
+## Footguns
+
+- **Match gene identifiers.** The net's `target` symbols must match your data's
+  `var_names` (both symbols, or both Ensembl IDs). Mismatched namespaces yield empty
+  results after `tmin` pruning. `dc.ds.ensmbl_to_symbol` helps convert.
+- **Right net for the method.** Weighted methods on an unweighted net treat all
+  weights as 1 (loses activation/repression sign); gene-set methods on a weighted net
+  ignore the weights. Pick the pair deliberately — see [methods.md](methods.md).
+- **`organism` must match your data.** Don't score mouse data with a human net;
+  fetch with `organism="mouse"` or `dc.op.translate`.
+
+## Related references
+
+[methods.md](methods.md) (weighted vs gene-set methods),
+[calling-convention.md](calling-convention.md) (how `net` is consumed, `tmin`),
+[getting-started.md](getting-started.md) (end-to-end example).
diff --git a/src/decoupler/_skills/data/references/ranking-and-plotting.md b/src/decoupler/_skills/data/references/ranking-and-plotting.md
new file mode 100644
index 0000000..251df37
--- /dev/null
+++ b/src/decoupler/_skills/data/references/ranking-and-plotting.md
@@ -0,0 +1,131 @@
+# Ranking and plotting activity scores
+
+After scoring, you usually want to (a) find which programs are characteristic of which
+groups, and (b) visualize them. Ranking lives in `dc.tl`, plotting in `dc.pl`.
+
+Most of these operate on the **score matrix as an AnnData** — get one with
+`dc.pp.get_obsm` (see [io-and-outputs.md](io-and-outputs.md)):
+
+```python
+dc.mt.ulm(adata, net, tmin=3)
+acts = dc.pp.get_obsm(adata, key="score_ulm")   # obs × sources AnnData
+```
+
+## Ranking (`dc.tl`)
+
+- `dc.tl.rankby_group(acts, groupby, reference="rest", method="t-test_overestim_var")`
+  → DataFrame ranking sources per group (which TFs/pathways mark each cell type).
+- `dc.tl.rankby_obsm(adata, key, uns_key="rank_obsm")` → rank the sources in an
+  `.obsm` score matrix across obs groups; can store into `.uns`.
+- `dc.tl.rankby_order(adata, order, stat="dcor")` → associate sources with a
+  continuous ordering (e.g. pseudotime) via a dependence statistic.
+
+```python
+df = dc.tl.rankby_group(acts, groupby="cell_type")   # columns include group, name, stat, ...
+```
+
+## Plotting (`dc.pl`)
+
+Each returns a matplotlib Figure (or `None`). Common ones:
+
+- `dc.pl.barplot(data, name, top=25)` — top scoring sources for one group/contrast.
+  `data` is a **wide `DataFrame`** (rows = groups/contrasts in `.index`, cols = sources)
+  and `name` is one **row label**, **not** the `get_obsm` AnnData. You get such a frame
+  for free when scoring a contrast/pseudobulk matrix (DataFrame input → wide `scores`),
+  or from per-cell scores via a group-mean aggregation (see the example below).
+- `dc.pl.dotplot(df, x, y, c, s, top=10)` — dots for a **single contrast**: one dot per
+  source, positioned at `x`, labelled by `y`, colored by `c`, sized by `s`. **`x`, `c`,
+  and `s` must be numeric columns** — `x` is sorted by `|x|` to pick the `top` rows, so a
+  categorical `x` errors. It is *not* a sources × groups grid; filter a `rankby_group`
+  result to one group first.
+- `dc.pl.volcano(data, x, y, ...)` — volcano of a differential result; pass `net` +
+  `name` to highlight a program's target genes.
+- `dc.pl.network(net, sources=5, targets=10, score=...)` — draw a source→target graph
+  (needs `igraph`); overlay `score`/`data` to color nodes.
+- `dc.pl.obsm(adata, key="rank_obsm")` — heatmap of a ranked `.obsm` matrix.
+- `dc.pl.source_targets(...)`, `dc.pl.order(...)`, `dc.pl.leading_edge(...)` — target-
+  level and ordering/GSEA-style views.
+
+```python
+# rank cols: group, reference, name, stat, meanchange, pval, padj
+df = dc.tl.rankby_group(acts, groupby="cell_type")
+
+# dotplot: ONE group; x/c/s numeric (x sorted by |x| for `top`), y = source label
+import numpy as np
+sub = df[df["group"] == "B cells"].copy()
+sub["sig"] = -np.log10(sub["padj"].clip(lower=2.22e-16))
+dc.pl.dotplot(sub, x="stat", y="name", c="meanchange", s="sig", top=10)
+
+# barplot: wide DataFrame (rows=groups, cols=sources), name = a row label
+mean_acts = acts.to_df().groupby(acts.obs["cell_type"], observed=True).mean()
+dc.pl.barplot(mean_acts, name="B cells", top=15)
+```
+
+Column names for `dotplot`/`volcano` must match the columns of the DataFrame you pass
+(inspect `df.columns` from `rankby_group` first), and `dotplot`'s `x`/`c`/`s` columns
+must be **numeric** — a categorical `x` (e.g. `x="group"`) raises "category cannot
+perform absolute".
+
+## Visualizing scores with scanpy (UMAP / spatial / H&E)
+
+`dc.pl.*` has **no UMAP or spatial plot**. To overlay a program's activity on an
+embedding or on tissue, do **not** hand-roll the plot and do **not** copy a single score
+into `adata.obs`. Instead, turn the scores into a scanpy-ready AnnData with
+`dc.pp.get_obsm` and let scanpy do it. `get_obsm` puts the score matrix in `.X` (so each
+**source/TF becomes a `var_name`**) while **preserving `obs`, `uns` (H&E images), and
+`obsm` (spatial coords)** — exactly what `sc.pl.umap` / `sc.pl.spatial` need.
+
+```python
+import scanpy as sc
+
+dc.mt.ulm(adata, net)                          # adata.obsm["score_ulm"]
+acts = dc.pp.get_obsm(adata, "score_ulm")      # scores in .X; keeps obs/uns/obsm
+
+# embedding overlay: "RFXAP" is a var_name -> colored from .X (its activity)
+sc.pl.umap(acts, color=["RFXAP"], cmap="RdBu_r")
+
+# spatial / H&E (e.g. dc.ds.msvisium): mix var names and obs columns in one call
+sc.pl.spatial(acts, color=["RFXAP", "niches"], cmap="RdBu_r", size=1.5, wspace=0.3)
+```
+
+`color=["RFXAP", "niches"]` colors the first panel by RFXAP activity (from `.X`) and the
+second by the `niches` obs column — scanpy builds the multi-panel figure itself.
+
+### Spatial/scanpy footguns
+
+- **Do NOT pass `ax=` (from your own `plt.subplots`) into `sc.pl.spatial`.** It manages
+  its own figure and colorbar, and passing `ax=` silently breaks `vmin`/`vmax`/`vcenter`
+  (the color scale collapses, e.g. to ±0.1, and spots may not render). For multiple
+  panels pass a **`color=[...]` list** and use `wspace=`/`ncols=` for layout.
+- **Activities are signed and often asymmetric** (e.g. −2 to +9). Use a diverging cmap
+  (`RdBu_r`). If you want an explicit scale, set a **symmetric robust limit** rather than
+  `vcenter=`:
+  ```python
+  import numpy as np
+  lim = float(np.percentile(np.abs(acts[:, "RFXAP"].X), 99))
+  sc.pl.spatial(acts, color="RFXAP", cmap="RdBu_r", vmin=-lim, vmax=lim, size=1.5)
+  ```
+- `dc.ds.msvisium` is **raw counts** — normalize before scoring (see
+  [datasets.md](datasets.md)).
+
+## Footguns
+
+- **Use existing plotting functions — never reimplement a plot.** For summary views use
+  `dc.pl.*`; for embeddings/tissue use scanpy `sc.pl.umap` / `sc.pl.spatial` (or
+  `squidpy.pl.spatial_scatter`). Do not hand-blit the H&E image + a scatter of spots; if
+  a scanpy call misbehaves, fix the call (usually a stray `ax=`), don't replace it.
+- `dc.pl.*` functions take **DataFrames**, not the raw `get_obsm` AnnData. But the two
+  take *different* frames: `barplot` wants a **wide** group×source matrix (aggregate
+  per-cell scores with `acts.to_df().groupby(...).mean()`, or use a wide `scores` from a
+  contrast run), while `dotplot` wants a **long** rank table (a `rankby_group` result,
+  filtered to one group).
+- `dc.pl.network` requires `igraph`; install the `full` extra if missing.
+- For `dotplot`, the `x`/`y`/`c`/`s` arguments are **column names** — passing values or
+  wrong names errors; check `df.columns` from the ranking step. `x`/`c`/`s` must be
+  **numeric** (`x` is sorted by `|x|` internally), so a categorical column there fails.
+
+## Related references
+
+[io-and-outputs.md](io-and-outputs.md) (`get_obsm`, score layout),
+[getting-started.md](getting-started.md) (full flow),
+[methods.md](methods.md) (what the scores mean per method).
diff --git a/src/decoupler/_skills/install.py b/src/decoupler/_skills/install.py
new file mode 100644
index 0000000..42dc5a0
--- /dev/null
+++ b/src/decoupler/_skills/install.py
@@ -0,0 +1,143 @@
+"""Install the bundled decoupler Agent Skill for Claude Code.
+
+This module exposes the ``decoupler-install-skills`` console script. It copies
+the Agent Skill that ships inside the installed package into a location where
+Claude Code discovers skills, so that a coding agent gains decoupler-specific
+guidance (API contracts and footguns) in every project, not just this repository.
+
+Claude Code does not scan Python ``site-packages`` for skills; it only looks in a
+small set of fixed roots. The personal root, ``~/.claude/skills/``, makes a skill
+available across all of a user's projects, so that is the default destination. The
+skill itself is one self-contained directory -- a ``SKILL.md`` router plus a
+``references/`` folder of detailed topic files that are read on demand -- bundled
+here under ``data/``, so installation is simply copying that directory.
+
+The copy is deliberately opt-in via this command rather than being performed on
+import or install, because silently writing into a user's home configuration would
+be surprising and invasive.
+
+Run ``decoupler-install-skills`` to install, or ``--print-path`` to see where the
+bundled source lives (useful with the ``CLAUDE_SKILLS_PATH`` environment variable
+if you would rather point Claude Code at the package in place than copy it).
+"""
+
+import argparse
+import shutil
+import sys
+from pathlib import Path
+
+SKILL_NAME = "decoupler"
+
+
+def _bundled_skill_dir() -> Path:
+    """Return the path to the Agent Skill bundled inside the package.
+
+    Returns
+    -------
+    skill_dir : pathlib.Path
+        The ``data/`` directory containing the bundled ``SKILL.md`` and
+        ``references/``.
+    """
+    return Path(__file__).resolve().parent / "data"
+
+
+def _default_dest() -> Path:
+    """Return the default install destination in the user's personal skill root.
+
+    Returns
+    -------
+    dest : pathlib.Path
+        ``~/.claude/skills/decoupler``.
+    """
+    return Path.home() / ".claude" / "skills" / SKILL_NAME
+
+
+def install_skill(dest: Path | None = None, force: bool = False) -> Path:
+    """Copy the bundled decoupler skill into a Claude Code skills directory.
+
+    The bundled skill directory is copied wholesale to ``dest``. If ``dest``
+    already exists it is left untouched unless ``force`` is set, in which case it
+    is removed and replaced so that an upgrade reflects the installed package
+    version exactly.
+
+    Parameters
+    ----------
+    dest : pathlib.Path or None, optional
+        The directory to install the skill into. If None, defaults to
+        ``~/.claude/skills/decoupler``. Default is None.
+    force : bool, optional
+        Whether to overwrite an existing installation at ``dest``. Default is
+        False.
+
+    Returns
+    -------
+    dest : pathlib.Path
+        The directory the skill was installed into.
+    """
+    src = _bundled_skill_dir()
+    if not (src / "SKILL.md").is_file():
+        raise FileNotFoundError(
+            f"Bundled skill not found at {src}. The package may be installed without its skill data."
+        )
+
+    if dest is None:
+        dest = _default_dest()
+
+    if dest.exists():
+        if not force:
+            raise FileExistsError(f"{dest} already exists. Re-run with --force to overwrite.")
+        shutil.rmtree(dest)
+
+    dest.parent.mkdir(parents=True, exist_ok=True)
+    shutil.copytree(src, dest, ignore=shutil.ignore_patterns(".ipynb_checkpoints", "__pycache__"))
+    return dest
+
+
+def main(argv: list[str] | None = None) -> int:
+    """Command-line entry point for ``decoupler-install-skills``.
+
+    Parameters
+    ----------
+    argv : list of str or None, optional
+        Arguments to parse. If None, ``sys.argv`` is used. Default is None.
+
+    Returns
+    -------
+    code : int
+        Process exit code; 0 on success and 1 on a handled error.
+    """
+    parser = argparse.ArgumentParser(
+        prog="decoupler-install-skills",
+        description="Install the decoupler Agent Skill for Claude Code.",
+    )
+    parser.add_argument(
+        "--dest",
+        type=Path,
+        default=None,
+        help="Destination directory (default: ~/.claude/skills/decoupler).",
+    )
+    parser.add_argument("--force", action="store_true", help="Overwrite an existing installation.")
+    parser.add_argument(
+        "--print-path",
+        action="store_true",
+        help="Print the bundled skill directory and exit without installing.",
+    )
+    args = parser.parse_args(argv)
+
+    if args.print_path:
+        print(_bundled_skill_dir())
+        return 0
+
+    try:
+        dest = install_skill(dest=args.dest, force=args.force)
+    except (FileExistsError, FileNotFoundError) as e:
+        print(f"error: {e}", file=sys.stderr)
+        return 1
+
+    print(f"Installed decoupler skill to {dest}")
+    print("It will be available to Claude Code in your next session.")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tests/test_install_skills.py b/tests/test_install_skills.py
new file mode 100644
index 0000000..e77549e
--- /dev/null
+++ b/tests/test_install_skills.py
@@ -0,0 +1,158 @@
+import re
+import subprocess
+import sys
+from pathlib import Path
+
+import pytest
+
+from decoupler._skills import install as install_mod
+from decoupler._skills.install import (
+    _bundled_skill_dir,
+    _default_dest,
+    install_skill,
+)
+
+###
+# Helpers
+
+
+def _link_targets(text):
+    """Return the basenames of every Markdown link to a .md file in `text`."""
+
+    return [m.split("/")[-1] for m in re.findall(r"\]\(([^)]+\.md)\)", text)]
+
+
+def _seed_fake_skill(root):
+    """Create a minimal but well-formed skill directory at `root`."""
+
+    (root / "references").mkdir(parents=True)
+    (root / "SKILL.md").write_text("---\nname: x\ndescription: y\n---\n\nbody\n")
+    (root / "references" / "foo.md").write_text("# foo\n")
+    return root
+
+
+###
+# Bundled-skill integrity (run against the data that actually ships)
+
+
+def test_bundled_skill_dir_is_well_formed():
+    data = _bundled_skill_dir()
+    assert data.is_dir()
+    assert (data / "SKILL.md").is_file()
+
+    refs = list((data / "references").glob("*.md"))
+    assert len(refs) > 0
+
+
+def test_skill_frontmatter():
+    text = (_bundled_skill_dir() / "SKILL.md").read_text()
+    assert text.startswith("---\n")
+
+    frontmatter = text.split("---\n", 2)[1]
+    fields = dict(re.findall(r"^(\w+):\s*(.*)$", frontmatter, flags=re.MULTILINE))
+
+    assert fields.get("name") == "decoupler"
+    # Claude Code caps the description (combined with when_to_use) at 1536 chars.
+    assert 0 < len(fields.get("description", "")) <= 1536
+
+
+def test_internal_links_all_resolve():
+    data = _bundled_skill_dir()
+    refs = data / "references"
+    existing = {p.name for p in refs.glob("*.md")}
+
+    broken = []
+    for f in [data / "SKILL.md", *refs.glob("*.md")]:
+        for target in _link_targets(f.read_text()):
+            if target not in existing:
+                broken.append(f"{f.name} -> {target}")
+
+    assert broken == [], f"broken internal links: {broken}"
+
+
+###
+# install_skill behavior (always into a temporary destination)
+
+
+def test_install_into_fresh_dest(tmp_path):
+    dest = tmp_path / "decoupler"
+    result = install_skill(dest=dest)
+
+    assert result == dest
+    assert (dest / "SKILL.md").is_file()
+
+    bundled = {p.name for p in (_bundled_skill_dir() / "references").glob("*.md")}
+    installed = {p.name for p in (dest / "references").glob("*.md")}
+    assert installed == bundled
+
+
+def test_install_existing_dest_raises_without_force(tmp_path):
+    dest = tmp_path / "decoupler"
+    install_skill(dest=dest)
+
+    with pytest.raises(FileExistsError):
+        install_skill(dest=dest)
+
+
+def test_install_force_replaces_stale_content(tmp_path):
+    dest = tmp_path / "decoupler"
+    install_skill(dest=dest)
+
+    junk = dest / "references" / "stale.md"
+    junk.write_text("remove me")
+
+    install_skill(dest=dest, force=True)
+
+    assert not junk.exists()
+    assert (dest / "SKILL.md").is_file()
+
+
+def test_install_missing_bundle_raises(tmp_path, monkeypatch):
+    empty = tmp_path / "empty"
+    empty.mkdir()
+    monkeypatch.setattr(install_mod, "_bundled_skill_dir", lambda: empty)
+
+    with pytest.raises(FileNotFoundError):
+        install_skill(dest=tmp_path / "decoupler")
+
+
+def test_install_excludes_checkpoints_and_pycache(tmp_path, monkeypatch):
+    src = _seed_fake_skill(tmp_path / "src")
+    (src / ".ipynb_checkpoints").mkdir()
+    (src / ".ipynb_checkpoints" / "SKILL-checkpoint.md").write_text("nope")
+    (src / "__pycache__").mkdir()
+    (src / "__pycache__" / "x.pyc").write_text("nope")
+    monkeypatch.setattr(install_mod, "_bundled_skill_dir", lambda: src)
+
+    dest = install_skill(dest=tmp_path / "decoupler")
+
+    assert not (dest / ".ipynb_checkpoints").exists()
+    assert not (dest / "__pycache__").exists()
+    assert (dest / "references" / "foo.md").is_file()
+
+
+def test_default_dest_is_under_home_and_isolated(tmp_path, monkeypatch):
+    # Patch home so the default-destination path never touches the real ~/.claude.
+    monkeypatch.setattr(Path, "home", lambda: tmp_path)
+
+    result = install_skill()
+
+    assert result == _default_dest()
+    assert result == tmp_path / ".claude" / "skills" / "decoupler"
+    assert (result / "SKILL.md").is_file()
+
+
+###
+# Console-script entry point (opt-in; does not install anything)
+
+
+def test_console_script_print_path():
+    # The entry point is installed next to the running interpreter (bin/Scripts).
+    exe = Path(sys.executable).parent / "decoupler-install-skills"
+    if not exe.exists():
+        pytest.skip("decoupler-install-skills entry point is not installed")
+
+    result = subprocess.run([str(exe), "--print-path"], capture_output=True, text=True)
+
+    assert result.returncode == 0
+    assert result.stdout.strip().endswith("data")