Gseapy Gene Enrichment

pip install gseapy

# Verify
python -c "import gseapy; print(gseapy.__version__)"
# 1.1.3

decisions:
  - id: D1
    param: analysis_mode
    kind: required
    source: upstream
    ask: "Test only the genes that passed a significance cut (ORA), or the whole ranked list without a cut (GSEA prerank)?"
    default: "prerank when the upstream stage yields a full ranking; ORA when it yields a gene list"

  - id: D2
    param: organism
    kind: derived
    source: upstream
    ask: "Which species are these genes from?"
    default: "carried from the upstream dataset"

  - id: D3
    param: gene_sets
    kind: required
    source: user
    depends_on: [D2]
    ask: "Which collections should the genes be tested against - GO terms, pathways, disease or perturbation signatures?"
    default: null

  - id: D4
    param: rnk
    kind: optional
    source: user
    depends_on: [D1]
    ask: "What should order the genes - effect size, signed significance, or the test statistic?"
    default: "log2FoldChange"
    skip_if: "ORA mode - ranking is not used"

  - id: D5
    param: cutoff
    kind: optional
    source: user
    ask: "How strong must a term's evidence be before it is reported?"
    default: 0.05

  - id: D6
    param: min_size, max_size
    kind: optional_conditional
    source: user
    ask: "Should very small or very broad gene sets be excluded as uninformative?"
    default: "15-500 genes"

  - id: D7
    param: permutation_num
    kind: optional_conditional
    source: user
    ask: "More permutations give finer p-values at proportionally more runtime - is the default enough?"
    default: 1000

  - id: D8
    param: weighted_score_type
    kind: optional
    source: user
    ask: "Should GSEA use the standard weighted enrichment score, or a different weighting for this ranking?"
    default: 1

  - id: D9
    param: seed, threads
    kind: never_ask
    source: data
    reason: "A fixed seed and available CPU cores affect reproducibility and runtime, not the enrichment model."
    default: "seed fixed, threads min(4, available_cores)"

import gseapy as gp

# ORA: test a gene list against GO Biological Process
gene_list = ["TP53", "BRCA1", "CDK2", "CCND1", "MYC", "EGFR", "KRAS", "PTEN"]

enr = gp.enrichr(gene_list=gene_list,
                 gene_sets=["GO_Biological_Process_2023"],
                 organism="human",
                 outdir=None)
print(enr.results.head(5)[["Term", "P-value", "Adjusted P-value", "Genes"]])

import gseapy as gp
import pandas as pd

# Gene list from DESeq2 (significant upregulated genes)
sig_genes = ["TP53", "BRCA1", "CDK2", "CCND1", "MYC", "EGFR",
             "KRAS", "PTEN", "RB1", "AKT1", "PIK3CA", "MDM2"]

# Run ORA against multiple databases
enr = gp.enrichr(
    gene_list=sig_genes,
    gene_sets=[
        "GO_Biological_Process_2023",
        "KEGG_2021_Human",
        "MSigDB_Hallmark_2020",
        "Reactome_2022",
    ],
    organism="human",
    outdir="enrichr_results/",
    cutoff=0.05,
)

# Display top results
results = enr.results
print(f"Enriched terms: {len(results[results['Adjusted P-value'] < 0.05])}")
print(results[results["Adjusted P-value"] < 0.05].sort_values("Adjusted P-value")
      .head(10)[["Gene_set", "Term", "Adjusted P-value", "Combined Score"]])

import gseapy as gp

# List all available gene set libraries
libraries = gp.get_library_name(organism="human")
print(f"Available databases: {len(libraries)}")
print("Selected databases:")
for lib in sorted(libraries):
    if any(kw in lib for kw in ["GO_Bio", "KEGG", "Hallmark", "Reactome"]):
        print(f"  {lib}")

# Mouse databases
mouse_libs = gp.get_library_name(organism="mouse")
print(f"\nMouse databases: {len(mouse_libs)}")

import gseapy as gp
import pandas as pd
import numpy as np

# Load DESeq2 results (or create example ranked list)
# deseq_results = pd.read_csv("deseq2_results.tsv", sep="\t", index_col=0)
# ranked = deseq_results["log2FoldChange"].dropna().sort_values(ascending=False)

# Example ranked gene list (gene → log2FC)
np.random.seed(42)
gene_names = [f"GENE_{i}" for i in range(1000)]
log2fc = np.random.normal(0, 2, 1000)
ranked = pd.Series(log2fc, index=gene_names).sort_values(ascending=False)

# Run preranked GSEA against MSigDB Hallmarks
pre_res = gp.prerank(
    rnk=ranked,
    gene_sets="MSigDB_Hallmark_2020",
    threads=4,
    min_size=15,
    max_size=500,
    permutation_num=1000,
    outdir="gsea_results/prerank/",
    seed=42,
    verbose=True,
)

# View results
res_df = pre_res.res2d
sig = res_df[res_df["FDR q-val"] < 0.25]
print(f"Significant gene sets (FDR < 0.25): {len(sig)}")
print(sig.sort_values("NES", ascending=False)[["Term", "NES", "NOM p-val", "FDR q-val"]].head(10))

from gseapy import gseaplot

# pre_res: prerank result from Step 3
top_term = pre_res.res2d.sort_values("NES", ascending=False).index[0]
gseaplot(
    rank_metric=pre_res.ranking,
    term=top_term,
    **pre_res.results[top_term],
    ofname="figures/gsea_running_score.png",
)
print(f"Saved figures/gsea_running_score.png ({top_term})")

import gseapy as gp

# enr.results is the enrichment table: Term, Overlap, Adjusted P-value, Genes, ...
enr = gp.enrichr(
    gene_list=["TP53", "BRCA1", "CDK2", "CCND1", "MYC", "EGFR",
               "KRAS", "PTEN", "RB1", "AKT1", "PIK3CA", "MDM2",
               "BCL2", "CDKN1A", "E2F1", "CCNE1"],
    gene_sets=["KEGG_2021_Human"],
    organism="human",
    outdir=None,
    cutoff=0.05,
)
enr.results.to_csv("enrichment.csv", index=False)
print(f"ORA terms: {len(enr.results)} -> enrichment.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "GSEA dot plot" recipe.

import gseapy as gp
import pandas as pd

# From DESeq2 output loaded into Python
# deseq_df = pd.read_csv("deseq2_results.tsv", sep="\t", index_col=0)
# deseq_df = deseq_df.dropna(subset=["log2FoldChange", "padj"])

# Simulate DESeq2 output
import numpy as np
np.random.seed(0)
n = 500
deseq_df = pd.DataFrame({
    "log2FoldChange": np.random.normal(0, 1.5, n),
    "padj": np.random.uniform(0, 1, n),
}, index=[f"GENE{i}" for i in range(n)])

# Significant up/down gene lists for ORA
up_genes = deseq_df[(deseq_df["padj"] < 0.05) & (deseq_df["log2FoldChange"] > 1)].index.tolist()
dn_genes = deseq_df[(deseq_df["padj"] < 0.05) & (deseq_df["log2FoldChange"] < -1)].index.tolist()
print(f"Upregulated: {len(up_genes)}, Downregulated: {len(dn_genes)}")

# ORA on upregulated genes
if up_genes:
    enr_up = gp.enrichr(gene_list=up_genes,
                         gene_sets=["GO_Biological_Process_2023", "KEGG_2021_Human"],
                         organism="human", outdir=None)
    sig_up = enr_up.results[enr_up.results["Adjusted P-value"] < 0.05]
    print(f"Enriched terms (upregulated): {len(sig_up)}")
    print(sig_up.sort_values("Adjusted P-value").head(5)[["Term", "Adjusted P-value"]])

# Preranked GSEA on full ranked list
ranked = deseq_df["log2FoldChange"].sort_values(ascending=False)
pre = gp.prerank(rnk=ranked, gene_sets="MSigDB_Hallmark_2020",
                 threads=4, permutation_num=500, outdir="gsea_out/", seed=42)
print(pre.res2d[pre.res2d["FDR q-val"] < 0.25].sort_values("NES", ascending=False)
      .head(5)[["Term", "NES", "FDR q-val"]])