Skills › Scientific Agent Skills › pathway-enrichment › run_enrichment.py
run_enrichment.py
266 lines · skills/pathway-enrichment/scripts/run_enrichment.py at 5f2bc03
1#!/usr/bin/env python3
2"""Run ORA or classic permutation preranked GSEA with GSEApy 1.3.1.
4Inputs must already use the organism and identifier namespace of the gene sets.
5No identifier conversion or orthology mapping is performed. CSV/TSV gene lists
6have a header; rank files have two comma/tab-separated columns, header optional.
7Local GMT + explicit ORA background or a rank file supports offline analysis.
9Examples:
10 python run_enrichment.py ora --genes hits.txt --background tested.txt \
11 --libraries pathways.gmt --outdir results
12 python run_enrichment.py gsea --rnk genes.rnk --libraries pathways.gmt \
13 --seed 123 --outdir results
15Requires gseapy==1.3.1, numpy, pandas, matplotlib. Online Enrichr submits query
16lists (and an optional background) to the public service.
17"""
18from __future__ import annotations
20import argparse
21import hashlib
22import json
23import sys
24import warnings
25from datetime import datetime, timezone
26from pathlib import Path
28import numpy as np
29import pandas as pd
31try:
32 import gseapy as gp
33except ImportError:
34 sys.exit("[FAIL] Install gseapy==1.3.1 in an isolated environment.")
36DEFAULT_LIBRARIES = ["MSigDB_Hallmark_2020", "GO_Biological_Process_2026"]
37_MISSING = {"", "nan", "none", "na", "n/a"}
40def _clean_index(index, organism: str):
41 """Trim IDs; casing is not an identifier or orthology mapping algorithm."""
42 return pd.Index(["" if pd.isna(g) else str(g).strip() for g in index])
45def _clean_symbols(genes, organism: str):
46 return list(dict.fromkeys(g for g in _clean_index(genes, organism)
47 if g.lower() not in _MISSING))
50def _read_gene_list(path: Path):
51 """CSV/TSV requires a header; text files contain one identifier per line."""
52 if path.suffix.lower() in {".csv", ".tsv"}:
53 sep = "\t" if path.suffix.lower() == ".tsv" else ","
54 return pd.read_csv(path, sep=sep, dtype=str, keep_default_na=False).iloc[:, 0].tolist()
55 return path.read_text(encoding="utf-8").splitlines()
58def _numeric(values):
59 # Missing estimates may be excluded; malformed numbers must never disappear.
60 values = pd.Series(values, copy=True)
61 missing = values.astype(str).str.strip().str.lower().isin(_MISSING) | values.isna()
62 return pd.to_numeric(values.mask(missing), errors="raise")
65def _prepare_rank(rnk: pd.Series, organism: str) -> pd.Series:
66 rnk = rnk.copy()
67 rnk.index = _clean_index(rnk.index, organism)
68 if any(g.lower() in _MISSING for g in rnk.index):
69 raise ValueError("Ranked genes must have nonempty identifiers; repair missing IDs.")
70 rnk = pd.Series(_numeric(rnk.to_numpy()).to_numpy(), index=rnk.index, dtype=float)
71 missing = int(rnk.isna().sum())
72 if missing:
73 warnings.warn(f"Excluded {missing} genes with missing ranking statistics.")
74 rnk = rnk.dropna()
75 if not np.isfinite(rnk.to_numpy()).all():
76 raise ValueError("Ranking statistics must be finite; repair infinity upstream.")
77 if rnk.index.has_duplicates:
78 raise ValueError("Duplicate gene IDs in ranking; resolve mappings/probes upstream.")
79 if len(rnk) < 2 or not rnk.ne(0).any():
80 raise ValueError("Ranking needs at least two genes and a nonzero score.")
81 ties = int(rnk.duplicated().sum())
82 if ties:
83 warnings.warn(f"{ties} tied scores; preserving input order within ties. Assess tie sensitivity.")
84 rnk = rnk.sort_values(ascending=False, kind="stable")
85 rnk.attrs["missing_statistics_excluded"] = missing
86 return rnk
89def _build_rank_from_deseq2(path: Path, organism: str) -> pd.Series:
90 """Use a signed Wald stat, else signed -log10(raw p), clipping p=0 to 1e-300."""
91 df = pd.read_csv(path, index_col=0, dtype=str, keep_default_na=False)
92 cols = {c.lower(): c for c in df.columns}
93 if "stat" in cols:
94 rnk = df[cols["stat"]]
95 metric = "supplied stat (must be signed Wald, not unsigned LRT)"
96 elif "log2foldchange" in cols and "pvalue" in cols:
97 lfc = _numeric(df[cols["log2foldchange"]])
98 pval = _numeric(df[cols["pvalue"]])
99 if ((pval.dropna() < 0) | (pval.dropna() > 1)).any():
100 raise ValueError("Raw p-values must lie in [0, 1].")
101 if not np.isfinite(lfc.dropna().to_numpy()).all():
102 raise ValueError("Fold changes must be finite.")
103 if pval.eq(0).any():
104 warnings.warn("Zero p-values clipped to 1e-300; prefer the original signed statistic.")
105 rnk = np.sign(lfc) * -np.log10(pval.clip(lower=1e-300))
106 metric = "sign(log2FoldChange) * -log10(raw pvalue clipped at 1e-300)"
107 else:
108 raise ValueError(f"Need signed 'stat' or 'log2FoldChange' + 'pvalue'; found {list(df.columns)}")
109 rnk = _prepare_rank(rnk, organism)
110 rnk.attrs["ranking_metric"] = metric
111 return rnk
114def _read_rnk(path: Path, organism: str) -> pd.Series:
115 lines = path.read_text(encoding="utf-8").splitlines()
116 if not lines:
117 raise ValueError("Ranked file is empty.")
118 first_line = lines[0]
119 sep = "\t" if "\t" in first_line else ","
120 df = pd.read_csv(path, sep=sep, header=None, dtype=str, keep_default_na=False)
121 if df.shape[1] != 2:
122 raise ValueError("Ranked file must have exactly two columns: gene,score")
123 if df.iloc[0, 0].strip().lower() in {"gene", "gene_id", "symbol"} and df.iloc[0, 1].strip().lower() in {"score", "stat", "rank"}:
124 df = df.iloc[1:]
125 return _prepare_rank(pd.Series(df[1].to_numpy(), index=df[0]), organism)
128def _validate_libraries(libraries, organism):
129 if not libraries:
130 if organism != "human":
131 raise ValueError("Nonhuman analyses require explicit --libraries matching the organism.")
132 libraries = DEFAULT_LIBRARIES.copy()
133 if len(libraries) != len(set(libraries)):
134 raise ValueError("Duplicate libraries are not allowed.")
135 remote = []
136 local_names = []
137 for lib in libraries:
138 if Path(lib).is_file():
139 if Path(lib).suffix.lower() != ".gmt":
140 raise ValueError(f"Local library must be a GMT file: {lib}")
141 local_names.append(Path(lib).name)
142 elif lib.lower().endswith(".gmt"):
143 raise ValueError(f"GMT file does not exist: {lib}")
144 else:
145 remote.append(lib)
146 if len(local_names) != len(set(local_names)):
147 raise ValueError("Local GMT basenames must be unique to preserve library provenance.")
148 if remote:
149 available = set(gp.get_library_name(organism=organism))
150 invalid = set(remote) - available
151 if invalid:
152 raise ValueError(f"Unavailable libraries for {organism}: {sorted(invalid)}")
153 return libraries
156def _dotplot(df: pd.DataFrame, column: str, title: str, outpath: Path):
157 if df.empty:
158 print("[OK] No significant terms; no significance plot written.")
159 return False
160 try:
161 import matplotlib.pyplot as plt
162 ax = gp.dotplot(df, column=column, title=title, top_term=15, cutoff=1.0)
163 fig = ax.get_figure()
164 fig.savefig(outpath, dpi=200, bbox_inches="tight")
165 plt.close(fig)
166 return True
167 except Exception as exc:
168 print(f"[WARN] Dotplot skipped: {exc}")
169 return False
172def run_ora(args):
173 genes = _clean_symbols(_read_gene_list(Path(args.genes)), args.organism)
174 if not genes:
175 raise ValueError("The hit list is empty after cleanup.")
176 background = None
177 if args.background:
178 background = _clean_symbols(_read_gene_list(Path(args.background)), args.organism)
179 if not background or not set(genes).issubset(background):
180 raise ValueError("Background must be nonempty and contain every query gene.")
181 elif any(Path(lib).is_file() for lib in args.libraries):
182 raise ValueError("Local ORA requires --background with the tested gene universe.")
183 else:
184 warnings.warn("Using Enrichr's service background, not an assay-specific universe.")
185 if background and args.organism not in {"human", "mouse"} and any(not Path(lib).is_file() for lib in args.libraries):
186 raise ValueError("Speedrichr background routing is not organism-specific; use local GMT for this organism.")
187 enr = gp.enrichr(gene_list=genes, gene_sets=args.libraries,
188 organism=args.organism, background=background, outdir=None)
189 if not isinstance(enr.results, pd.DataFrame):
190 raise ValueError("No overlapping terms returned; inspect namespace and library coverage.")
191 res = enr.results.copy()
192 sig = res[res["Adjusted P-value"] < args.fdr].sort_values("Adjusted P-value")
193 args.input_summary = {"query_genes": len(genes), "background_genes": len(background) if background else None,
194 "background_mode": "explicit" if background else "service-default",
195 "correction": "provider/local BH per library; inspect returned term family"}
196 return res, sig, "Adjusted P-value"
199def run_gsea(args):
200 rnk = (_build_rank_from_deseq2(Path(args.deseq2), args.organism) if args.deseq2
201 else _read_rnk(Path(args.rnk), args.organism))
202 pre = gp.prerank(rnk=rnk, gene_sets=args.libraries, organism=args.organism,
203 min_size=args.min_size, max_size=args.max_size,
204 permutation_num=args.permutations, seed=args.seed, threads=args.threads,
205 method="permutation", ascending=None, outdir=None)
206 res = pre.res2d.copy()
207 res["FDR q-val"] = pd.to_numeric(res["FDR q-val"], errors="raise")
208 sig = res[res["FDR q-val"] < args.fdr].sort_values("FDR q-val")
209 args.input_summary = {"ranked_genes": len(rnk), "tied_scores": int(rnk.duplicated().sum()),
210 "ranking_metric": rnk.attrs.get("ranking_metric", "supplied score"),
211 "missing_statistics_excluded": rnk.attrs["missing_statistics_excluded"],
212 "method": "permutation", "permutation_type": "gene_set", "weight": 1.0,
213 "correction": "GSEA permutation FDR per library prefix"}
214 return res, sig, "FDR q-val"
217def main():
218 p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
219 sub = p.add_subparsers(dest="method", required=True)
220 common = argparse.ArgumentParser(add_help=False)
221 common.add_argument("--libraries", nargs="+", help="Enrichr names or local GMT paths; human default: Hallmark 2020 + GO BP 2026.")
222 common.add_argument("--organism", default="human", choices=["human", "mouse", "fly", "yeast", "worm", "fish"])
223 common.add_argument("--outdir", default="enrichment_results")
224 common.add_argument("--fdr", type=float, default=0.05)
225 ora = sub.add_parser("ora", parents=[common], help="Over-representation analysis.")
226 ora.add_argument("--genes", required=True, help="Text gene list or CSV/TSV with header.")
227 ora.add_argument("--background", help="Tested gene universe; required with local GMT.")
228 gsea = sub.add_parser("gsea", parents=[common], help="Classic permutation preranked GSEA.")
229 source = gsea.add_mutually_exclusive_group(required=True)
230 source.add_argument("--deseq2", help="CSV with gene index and signed Wald stat (not LRT stat).")
231 source.add_argument("--rnk", help="Two columns: gene,score; comma or tab separated.")
232 gsea.add_argument("--min-size", type=int, default=15)
233 gsea.add_argument("--max-size", type=int, default=500)
234 gsea.add_argument("--permutations", type=int, default=1000)
235 gsea.add_argument("--seed", type=int, default=123)
236 gsea.add_argument("--threads", type=int, default=4)
237 args = p.parse_args()
238 try:
239 if not 0 < args.fdr <= 1:
240 raise ValueError("--fdr must be in (0, 1].")
241 if args.method == "gsea" and (args.min_size < 1 or args.max_size < args.min_size or args.permutations < 1 or args.threads < 1):
242 raise ValueError("Require 1 <= min-size <= max-size, permutations >= 1 and threads >= 1.")
243 args.libraries = _validate_libraries(args.libraries, args.organism)
244 outdir = Path(args.outdir)
245 outdir.mkdir(parents=True, exist_ok=True)
246 # Avoid stale figures/tables looking like outputs of this run.
247 if any(outdir.glob(f"{args.method}_*")):
248 raise ValueError("Output directory already contains results for this method; choose a fresh directory.")
249 res, sig, col = run_ora(args) if args.method == "ora" else run_gsea(args)
250 res.to_csv(outdir / f"{args.method}_results.csv", index=False)
251 sig.to_csv(outdir / f"{args.method}_significant.csv", index=False)
252 plotted = _dotplot(sig, col, args.method.upper(), outdir / f"{args.method}_dotplot.png")
253 inputs = [getattr(args, k, None) for k in ("genes", "background", "deseq2", "rnk")] + args.libraries
254 hashes = {str(f): hashlib.sha256(Path(f).read_bytes()).hexdigest() for f in inputs if f and Path(f).is_file()}
255 metadata = {"timestamp_utc": datetime.now(timezone.utc).isoformat(), "gseapy_version": gp.__version__,
256 "arguments": vars(args), "input_sha256": hashes, "results": len(res),
257 "significant": len(sig), "dotplot_written": plotted,
258 "identifier_policy": "trim only; preserve case; ranked duplicate IDs rejected"}
259 (outdir / f"{args.method}_metadata.json").write_text(json.dumps(metadata, indent=2)+"\n", encoding="utf-8")
260 print(f"[OK] {len(sig)}/{len(res)} terms pass {col} < {args.fdr}; outputs: {outdir}")
261 except (ValueError, OSError, IndexError) as exc:
262 p.error(str(exc))
265if __name__ == "__main__":
266 main()