tooluniverse-gene-enrichment

What this skill does

## COMPUTE, DON'T DESCRIBE
When analysis requires computation (statistics, data processing, scoring, enrichment), write and run Python code via Bash. Don't describe what you would do — execute it and report actual results. Use ToolUniverse tools to retrieve data, then Python (pandas, scipy, statsmodels, matplotlib) to analyze it.

# Gene Enrichment and Pathway Analysis

## RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:
- `*_executed.ipynb` → read with `tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}'` and cite its cell outputs as the authoritative answer
- Pre-computed enrichment files (CSV/TSV named `*enrich*`, `*go*`, `*kegg*`, `*reactome*`, `*ego*`, `*_simplified.csv`) → read directly
- Canonical analysis scripts (`analysis.R`, `run_*.py`, `find_*.R`, `*.Rmd`) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if **none** of the above exist. Re-running enrichment from raw DEG lists produces different numbers than the published answer due to subtle filter differences upstream, and is much slower.

---

## PRIMARY SCRIPTS — use these FIRST

Three deterministic CLI scripts cover the bulk of enrichment questions.
Each handles edge cases (ties at top, simplify-changes-padj, multi-condition
screening) that the agent tends to get wrong when writing ad-hoc code.
**Always write outputs to `/tmp/...` — never into the data folder.**

### 1. `scripts/gseapy_enrichment_runner.py` — gseapy enrichr / prerank

**When to use**: the question references `gseapy`, `enrichr`, "Enrichr library", or any GO BP/MF/CC, KEGG, Reactome, WikiPathways, MSigDB enrichment via the gseapy package.

```bash
python skills/tooluniverse-gene-enrichment/scripts/gseapy_enrichment_runner.py \
    --gene-list /tmp/sig_symbols.txt \
    --library GO_Biological_Process_2021,Reactome_2022 \
    --organism Human \
    --top 5 \
    --candidate "negative regulation of epithelial cell proliferation" \
    --workdir /tmp/gseapy_run
```

What it reports (parseable lines):
- `# TOP_BY_ADJ_PVALUE: <term>` — what `df.sort_values('Adjusted P-value').iloc[0]` returns (this is what published notebooks usually print)
- `# TIES_AT_TOP: n=K` — number of terms tied at the lowest Adjusted P-value
- `# TOP_TIE_BROKEN: <term>` — deterministic tie-break (adj_p, raw_p, overlap desc, alphabetic)
- `# TOPN_BY_ADJ_PVALUE:` — full top N listing
- `# CANDIDATE_RANK '<term>': rank=R adj_p=...` — for any `--candidate` substring you pass
- `# SUBSTRING_COUNT_TOPN '<sub>': K` — for `--count-substring` queries (e.g., "how many top-20 terms contain 'Oxidative'")

Pass `--mode prerank --ranked-list /tmp/lfc.tsv` for GSEA preranked.

### 2. `scripts/enrichgo_runner.py` — clusterProfiler::enrichGO + simplify

**When to use**: the question references `enrichGO`, `clusterProfiler`, `simplify`, `simplify(cutoff=0.7)`, or the data folder contains an `analysis.R` / `find_*.R` that uses these. This is the canonical R workflow — gseapy does NOT reproduce it faithfully because `simplify` changes the multiple-testing denominator and thus the p.adjust values for surviving terms.

```bash
python skills/tooluniverse-gene-enrichment/scripts/enrichgo_runner.py \
    --gene-list /tmp/sig_ensembl.txt \
    --background /tmp/bg_ensembl.txt \
    --keytype ENSEMBL \
    --ontology BP \
    --simplify-cutoff 0.7 \
    --candidate "regulation of T cell activation" \
    --candidate "potassium ion transmembrane transport" \
    --workdir /tmp/enrichgo_run
```

What it reports:
- `# TOP10_RAW:` — top 10 from `as.data.frame(ego)` (BEFORE simplify; raw p.adjust)
- `# TOP10_SIMPLIFIED:` — top 10 from `as.data.frame(simplify(ego, cutoff=0.7))` (AFTER simplify; p.adjust differs)
- `# CANDIDATE '<term>': raw_rank=R raw_padj=... simp_rank=R simp_padj=...` — both pre- and post-simplify ranks for each candidate. `simp_rank=NA (collapsed by simplify)` means the term was redundant with a more-significant parent/sibling and was dropped.

When a question says "in the simplified results" or "after simplify", read **simp_padj**. When it just says "the most enriched" without mentioning simplify, default to the simplified frame anyway IF the canonical `analysis.R` calls `simplify`.

Requires R packages `clusterProfiler`, `org.Hs.eg.db` (or `org.Mm.eg.db` for mouse). Install via `Rscript skills/evals/install_r_packages.R` if missing.

### 3. `scripts/condition_enrichment_screen.py` — per-condition enrichment

**When to use**: the question asks "what fraction/percentage of conditions/screens/timepoints/groups had significant enrichment of <category>", or you have an N-by-many gene table and need per-condition enrichment.

```bash
# Per-condition gene-list files:
python skills/tooluniverse-gene-enrichment/scripts/condition_enrichment_screen.py \
    --condition-genes acute=/tmp/acute_sig.txt \
    --condition-genes round1=/tmp/r1_sig.txt \
    --condition-genes round2=/tmp/r2_sig.txt \
    --condition-genes round3=/tmp/r3_sig.txt \
    --library /path/to/local_pathways.gmt \
    --background /tmp/expressed.txt \
    --keyword immune --keyword cytokine --keyword interferon \
    --workdir /tmp/cond_screen
```

Or pass a single 2-col TSV (`condition<TAB>gene`) via `--conditions-tsv`.

What it reports:
- Per condition: `n_genes`, `sig_terms` (Adj P < cutoff), `sig_terms_keyword` (sig terms whose Term contains any --keyword)
- `# n_with_any_sig=N pct_with_any_sig=N%` — the fraction with any significant term
- `# n_with_keyword_sig=N pct_with_keyword_sig=N%` — the fraction whose sig terms include a category keyword

Notes:
- The `--library` can be either an Enrichr library name (online) or a path to a local `.gmt` file. **Prefer the local GMT if the data folder ships one** (avoids rate-limits and exactly reproduces published results).
- Use `--exclude-condition <label>` for "control" / "baseline" conditions that the question wants excluded from the denominator.
- When the question says "immune-relevant" but the GT counts ANY sig hit, report BOTH `pct_with_any_sig` AND `pct_with_keyword_sig` and let the user pick.

### Why these scripts exist (debugging notes)

Enrichment top-hits depend critically on three things:
1. **Upstream DEG filter** (padj only? padj+|LFC|>0.5? +baseMean>10? lfc-shrunk?). The "right" filter is whatever the canonical notebook used. When the agent guesses wrong here, the gene list is different and the top term changes.
2. **Library snapshot** — Enrichr libraries get republished. `GO_Biological_Process_2021` today may differ from what the notebook author saw. There is NO good fix; report the candidate's rank and let the user judge.
3. **Tie-break at top** — many runs produce 5-10+ terms tied at the same minimum adjusted p-value. `df.sort_values(...).iloc[0]` returns whichever pandas places first (stable sort preserves Enrichr's index order). Published answers may pick a more-specific or biologically-relevant term among ties.

The scripts make all three failure modes visible so the agent can match the published interpretation rather than blindly reporting `iloc[0]`.

### When `# TIES_AT_TOP: n=N` is large (warning sign)

If `gseapy_enrichment_runner.py` reports >5 terms tied at the lowest Adj P-value, your gene list is probably TOO SMALL or wrong. Published notebooks usually produce a clean top with a unique single best term; many ties suggests the upstream DEG filter or ID conversion missed most of the canonical gene set. Re-check:
- Did you apply the SAME filter the notebook used? (padj only vs padj+|LFC|>thr vs +baseMean>10)
- Is your gene-ID space the same? (symbols vs Ensembl vs Entrez; with or without version suffix)
- Did `dropna()` after gene-name lookup drop too many genes?
Re-run after fixing and the ties at top should drop sharply.

### DEG filter default — use ONLY what the question names

When the question describes the input gene list, apply ONLY the thresholds it
names. Do NOT silently add `|LFC| > x`,