tooluniverse-gene-enrichment

Gene-set enrichment analysis — GO (Biological Process, Molecular Function, Cellular Component), KEGG, Reactome pathway enrichment via clusterProfiler, gseapy, ORA, GSEA. Use for interpreting DEG lists, screen hit lists, or any gene-list-to-pathways query. Includes simplify-cutoff handling and union-

By mims-harvard · 374 installs

npx skills add mims-harvard/tooluniverse --skill tooluniverse-gene-enrichment

Source repository · Upstream listing

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. 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. 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. 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 , baseMean y , or LFC shrinkage — extra filters shrink the gene list and change overlap counts. Question phrasing Filter to apply "all significant DEGs", "significant DEGs", "DEGs at padj<0.05" padj < 0.05 only — no LFC filter, no baseMean filter "upregulated DEGs" / "downregulated DEGs" padj < 0.05 + sign of log2FoldChange only "DEGs with \ LFC\ 1" or "fold change 2" padj < 0.05 + the stated LFC threshold "after LFC shrinkage" / "apeglm shrunk" Apply lfcShrink() ; otherwise do not Question mentions baseMean or "expressed genes" Apply the named cutoff; otherwise do not Cross check before reporting: count your filtered gene list and state it ( n sig=N in the report). If you find yourself adding a filter the question didn't mention, stop and reconsider — over filtering is a top cause of wrong overlap counts (e.g., reporting 20/64 when the answer is 22/64). Perform comprehensive gene enrichment analysis including Gene Ontology (GO), KEGG, Reactome, WikiPathways, and MSigDB enrichment using both Over Representation Analysis (ORA) and Gene Set Enrichment Analysis (GSEA). Integrates local computation via gseapy with ToolUniverse pathway databases for cross validated, publication ready results. IMPORTANT : Always use English terms in tool calls (gene names, pathway names, organism names), even if the user writes in another language. Only try original language terms as a fallback if English returns no results. Respond in the user's language. Domain Reasoning: Background Selection Enrichment results are only as good as your background. The default background (all annotated genes in the genome) inflates enrichment for tissue specific or context specific gene lists. Always consider: what is the appropriate background for this experiment? For brain RNA seq, use brain expressed genes as background; for a proteomics experiment, use detected proteins. A gene that is never expressed in your system cannot be a true negative control. LOOK UP DON'T GUESS: adjusted p values, gene set overlap counts, and which genes from your input list drive each enriched term. Always retrieve the inputGenes field from enrichment results — do not assume which genes caused a term to be significant. When a term looks surprising, verify by checking which genes overlap. When to Use This Skill Apply when users: Ask about gene enrichment analysis (GO, KEGG, Reactome, etc.) Have a gene list from differential expression, clustering, or any experiment Want to know which biological processes, molecular functions, or cellular components are enriched Need KEGG or Reactome pathway enrichment analysis Ask about GSEA (Gene Set Enrichment Analysis) with ranked gene lists Want over representation analysis (ORA) with Fisher's exact test Need multiple testing correction (Benjamini Hochberg, Bonferroni) Ask about enrichGO, gseapy, clusterProfiler style analyses NOT for (use other skills instead): Network pharmacology / drug repurposing → Use tooluniverse network pharmacology Disease characterization → Use tooluniverse multiomic disease characterization Single gene function lookup → Use tooluniverse disease research Spatial omics analysis → Use tooluniverse spatial omics analysis Protein protein interaction analysis only → Use tooluniverse protein interactions Input Parameters Parameter Required Description Example gene list Yes List of gene symbols, Ensembl IDs, or Entrez IDs ["TP53", "BRCA1", "EGFR"] organism No Organism (default: human). Supported: human, mouse, rat, fly, worm, yeast, zebrafish human analysis type No ORA (default) or GSEA ORA enrichment databases No Which databases to query. Default: all applicable ["GO BP", "GO MF", "GO CC", "KEGG", "Reactome"] gene id type No Input ID type: symbol , ensembl , entrez , uniprot (auto detected if omitted) symbol p value cutoff No Significance threshold (default: 0.05) 0.05 correction method No Multiple testing: BH (Benjamini Hochberg, default), bonferroni , fdr BH background genes No Custom background gene set (default: genome wide) ["GENE1", "GENE2", ...] ranked gene list No For GSEA: gene to score mapping (e.g., log2FC) {"TP53": 2.5, "BRCA1": 1.3, ...} Core Principles 1. Report first approach Create report file FIRST, then populate progressively 2. ID disambiguation FIRST Detect and convert gene IDs before ANY enrichment 3. Multi source validation Run enrichment on at least 2 independent tools, cross validate 4. Exact p values Report raw p values AND adjusted p values with correction method 5. Multiple testing correction ALWAYS apply Benjamini Hochberg unless user specifies otherwise 6. Gene set size filtering Filter by min/max gene set size to avoid trivial/overly broad terms 7. Evidence grading Grade enrichment sources T1 T4 8. Negative results documented "No significant enrichment" is a valid finding 9. Source references Every enrichment result must cite the tool/database/library used 10. Completeness checklist Mandatory section at end showing analysis coverage