404 lines
10 KiB
Markdown
404 lines
10 KiB
Markdown
---
|
|
title: "Enrichment Analysis with gseapy"
|
|
task: ""
|
|
lineage_type: import
|
|
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-rnaseq-deseq2/references/enrichment_analysis.md
|
|
upstream_sha: e2520a96
|
|
imported_at: 2026-06-26
|
|
prompt_class: prompt
|
|
upstream_changes: accepted
|
|
author: upstream
|
|
validated: false
|
|
---
|
|
|
|
# Enrichment Analysis with gseapy
|
|
|
|
Complete guide to pathway and GO enrichment analysis.
|
|
|
|
## Basic Over-Representation Analysis (ORA)
|
|
|
|
```python
|
|
import gseapy as gp
|
|
|
|
# Prepare gene list (from DESeq2 results)
|
|
sig_genes = results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 0.5)]
|
|
gene_list = sig_genes.index.tolist()
|
|
|
|
# Run enrichment
|
|
enr = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets='GO_Biological_Process_2023',
|
|
background=None, # or provide background gene list
|
|
outdir=None, # Don't save files
|
|
cutoff=0.05,
|
|
no_plot=True,
|
|
verbose=False
|
|
)
|
|
|
|
# Access results
|
|
results_df = enr.results
|
|
print(results_df.head(10))
|
|
```
|
|
|
|
## Gene Set Library Selection
|
|
|
|
### Gene Ontology (GO)
|
|
|
|
```python
|
|
# Human/Mouse GO Biological Process (most recent)
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')
|
|
|
|
# GO Molecular Function
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Molecular_Function_2021')
|
|
|
|
# GO Cellular Component
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Cellular_Component_2021')
|
|
```
|
|
|
|
### KEGG Pathways
|
|
|
|
```python
|
|
# Human KEGG
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')
|
|
|
|
# Mouse KEGG
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2019_Mouse')
|
|
```
|
|
|
|
### Reactome
|
|
|
|
```python
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='Reactome_2022')
|
|
```
|
|
|
|
### WikiPathways
|
|
|
|
```python
|
|
# Human
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='WikiPathways_2019_Human')
|
|
|
|
# Mouse
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='WikiPathways_2019_Mouse')
|
|
```
|
|
|
|
### Other Libraries
|
|
|
|
```python
|
|
# MSigDB Hallmark gene sets
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='MSigDB_Hallmark_2020')
|
|
|
|
# GWAS Catalog
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GWAS_Catalog_2019')
|
|
|
|
# BioCarta
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='BioCarta_2016')
|
|
```
|
|
|
|
## Using Background Gene Sets
|
|
|
|
```python
|
|
# Background = all genes tested in DESeq2
|
|
background = results.dropna(subset=['padj']).index.tolist()
|
|
|
|
enr = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets='GO_Biological_Process_2023',
|
|
background=background, # Provide background
|
|
outdir=None,
|
|
cutoff=0.05
|
|
)
|
|
```
|
|
|
|
## Extract Specific Results
|
|
|
|
```python
|
|
def extract_enrichment_answer(enr_results, term_query=None, metric='Adjusted P-value'):
|
|
"""Extract specific enrichment result.
|
|
|
|
Args:
|
|
enr_results: gseapy enrichment results DataFrame
|
|
term_query: String to search in Term column (case-insensitive)
|
|
metric: Column to return ('Adjusted P-value', 'Odds Ratio', 'P-value', etc.)
|
|
|
|
Returns:
|
|
Value or DataFrame of matches
|
|
"""
|
|
if term_query:
|
|
# Case-insensitive search
|
|
mask = enr_results['Term'].str.lower().str.contains(term_query.lower())
|
|
matches = enr_results[mask]
|
|
if len(matches) == 1:
|
|
return matches.iloc[0][metric]
|
|
elif len(matches) > 1:
|
|
return matches[['Term', metric]]
|
|
else:
|
|
return None
|
|
|
|
# Return top result
|
|
return enr_results.sort_values(metric).head(1)
|
|
|
|
# Usage
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')
|
|
pval = extract_enrichment_answer(enr.results, term_query='ABC transporters', metric='Adjusted P-value')
|
|
print(f"ABC transporters adjusted p-value: {pval}")
|
|
```
|
|
|
|
## Extract Gene Count in Pathway
|
|
|
|
```python
|
|
# Enrichment results have 'Overlap' column (e.g., "11/42")
|
|
# Also 'Genes' column with semicolon-separated gene list
|
|
|
|
pathway_row = enr.results[enr.results['Term'].str.contains('ABC transporters')].iloc[0]
|
|
|
|
# Number of DEGs in pathway
|
|
overlap_str = pathway_row['Overlap'] # e.g., "11/42"
|
|
n_overlap = int(overlap_str.split('/')[0]) # 11
|
|
|
|
# Pathway size
|
|
n_pathway = int(overlap_str.split('/')[1]) # 42
|
|
|
|
# Gene list
|
|
genes_in_pathway = pathway_row['Genes'].split(';')
|
|
print(f"{n_overlap} genes contribute to this pathway:")
|
|
print(genes_in_pathway)
|
|
```
|
|
|
|
## Multi-Library Enrichment
|
|
|
|
```python
|
|
# Run enrichment on multiple libraries
|
|
libraries = [
|
|
'GO_Biological_Process_2023',
|
|
'KEGG_2021_Human',
|
|
'Reactome_2022'
|
|
]
|
|
|
|
all_results = {}
|
|
for lib in libraries:
|
|
enr = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets=lib,
|
|
outdir=None,
|
|
cutoff=0.05,
|
|
no_plot=True,
|
|
verbose=False
|
|
)
|
|
all_results[lib] = enr.results
|
|
|
|
# Combine top results
|
|
combined = []
|
|
for lib, res in all_results.items():
|
|
top5 = res.head(5).copy()
|
|
top5['Library'] = lib
|
|
combined.append(top5)
|
|
combined_df = pd.concat(combined, ignore_index=True)
|
|
```
|
|
|
|
## GO Term Simplification
|
|
|
|
```python
|
|
def simplify_go_terms(enr_results, similarity_threshold=0.7):
|
|
"""Simplify GO terms by removing highly similar terms.
|
|
|
|
Approximation of R clusterProfiler::simplify().
|
|
Uses Jaccard similarity on gene sets.
|
|
"""
|
|
if len(enr_results) == 0:
|
|
return enr_results
|
|
|
|
# Parse gene sets from Genes column
|
|
terms = enr_results.sort_values('Adjusted P-value').copy()
|
|
gene_sets = {}
|
|
for _, row in terms.iterrows():
|
|
genes = set(row['Genes'].split(';'))
|
|
gene_sets[row['Term']] = genes
|
|
|
|
# Compute Jaccard similarity between terms
|
|
keep = []
|
|
removed = set()
|
|
|
|
for i, (term_i, genes_i) in enumerate(gene_sets.items()):
|
|
if term_i in removed:
|
|
continue
|
|
keep.append(term_i)
|
|
|
|
for term_j, genes_j in list(gene_sets.items())[i+1:]:
|
|
if term_j in removed:
|
|
continue
|
|
# Jaccard similarity
|
|
intersection = len(genes_i & genes_j)
|
|
union = len(genes_i | genes_j)
|
|
if union > 0:
|
|
similarity = intersection / union
|
|
if similarity > similarity_threshold:
|
|
removed.add(term_j) # Remove the less significant term
|
|
|
|
return terms[terms['Term'].isin(keep)]
|
|
|
|
# Usage
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')
|
|
simplified = simplify_go_terms(enr.results, similarity_threshold=0.7)
|
|
print(f"Original: {len(enr.results)} terms")
|
|
print(f"Simplified: {len(simplified)} terms")
|
|
```
|
|
|
|
## Gene Set Enrichment Analysis (GSEA)
|
|
|
|
```python
|
|
# GSEA requires ranked gene list (not just significant genes)
|
|
# Rank by -log10(pvalue) * sign(log2FC)
|
|
|
|
results_ranked = results.dropna(subset=['pvalue', 'log2FoldChange'])
|
|
results_ranked['rank'] = -np.log10(results_ranked['pvalue']) * np.sign(results_ranked['log2FoldChange'])
|
|
results_ranked = results_ranked.sort_values('rank', ascending=False)
|
|
|
|
# Create rank dictionary
|
|
rank_dict = dict(zip(results_ranked.index, results_ranked['rank']))
|
|
|
|
# Run GSEA
|
|
gsea_res = gp.prerank(
|
|
rnk=rank_dict,
|
|
gene_sets='KEGG_2021_Human',
|
|
outdir=None,
|
|
permutation_num=1000,
|
|
no_plot=True,
|
|
verbose=False
|
|
)
|
|
|
|
# Access results
|
|
gsea_df = gsea_res.res2d
|
|
print(gsea_df[gsea_df['FDR q-val'] < 0.05])
|
|
```
|
|
|
|
## Organism-Specific Libraries
|
|
|
|
### Human
|
|
|
|
```python
|
|
libraries_human = [
|
|
'GO_Biological_Process_2023',
|
|
'GO_Molecular_Function_2021',
|
|
'GO_Cellular_Component_2021',
|
|
'KEGG_2021_Human',
|
|
'Reactome_2022',
|
|
'WikiPathways_2019_Human',
|
|
'MSigDB_Hallmark_2020',
|
|
'BioCarta_2016'
|
|
]
|
|
```
|
|
|
|
### Mouse
|
|
|
|
```python
|
|
libraries_mouse = [
|
|
'GO_Biological_Process_2023',
|
|
'GO_Molecular_Function_2021',
|
|
'GO_Cellular_Component_2021',
|
|
'KEGG_2019_Mouse',
|
|
'WikiPathways_2019_Mouse'
|
|
]
|
|
```
|
|
|
|
### Other Organisms
|
|
|
|
For other organisms, use custom gene sets:
|
|
|
|
```python
|
|
# Load custom GMT file
|
|
enr = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets='/path/to/custom.gmt',
|
|
background=background
|
|
)
|
|
```
|
|
|
|
## Complete Example: DEG to Enrichment
|
|
|
|
```python
|
|
import pandas as pd
|
|
import gseapy as gp
|
|
from pydeseq2.dds import DeseqDataSet
|
|
from pydeseq2.ds import DeseqStats
|
|
|
|
# Run DESeq2 (from pydeseq2_workflow.md)
|
|
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", quiet=True)
|
|
dds.deseq2()
|
|
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
|
|
stat_res.run_wald_test()
|
|
results = stat_res.results_df
|
|
|
|
# Filter DEGs
|
|
sig_genes = results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 0.5)]
|
|
gene_list = sig_genes.index.tolist()
|
|
|
|
# Run GO enrichment
|
|
enr_go = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets='GO_Biological_Process_2023',
|
|
background=results.dropna(subset=['padj']).index.tolist(),
|
|
outdir=None,
|
|
cutoff=0.05,
|
|
no_plot=True,
|
|
verbose=False
|
|
)
|
|
|
|
# Run KEGG enrichment
|
|
enr_kegg = gp.enrich(
|
|
gene_list=gene_list,
|
|
gene_sets='KEGG_2021_Human',
|
|
background=results.dropna(subset=['padj']).index.tolist(),
|
|
outdir=None,
|
|
cutoff=0.05,
|
|
no_plot=True,
|
|
verbose=False
|
|
)
|
|
|
|
# Display top results
|
|
print("\nTop 5 GO terms:")
|
|
print(enr_go.results[['Term', 'Adjusted P-value', 'Overlap']].head(5))
|
|
|
|
print("\nTop 5 KEGG pathways:")
|
|
print(enr_kegg.results[['Term', 'Adjusted P-value', 'Overlap']].head(5))
|
|
|
|
# Answer specific question
|
|
if 'immune response' in question.lower():
|
|
immune_result = enr_go.results[enr_go.results['Term'].str.contains('immune', case=False)]
|
|
if len(immune_result) > 0:
|
|
answer = immune_result.iloc[0]['Adjusted P-value']
|
|
print(f"\nImmune response adjusted p-value: {answer}")
|
|
```
|
|
|
|
## Common Enrichment Patterns
|
|
|
|
### Pattern 1: Extract adjusted p-value for specific pathway
|
|
|
|
```python
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')
|
|
pathway = enr.results[enr.results['Term'].str.contains('ABC transporters')]
|
|
answer = round(pathway.iloc[0]['Adjusted P-value'], 4)
|
|
```
|
|
|
|
### Pattern 2: Count significant pathways
|
|
|
|
```python
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023', cutoff=0.05)
|
|
answer = len(enr.results[enr.results['Adjusted P-value'] < 0.05])
|
|
```
|
|
|
|
### Pattern 3: Gene count in pathway
|
|
|
|
```python
|
|
pathway = enr.results[enr.results['Term'].str.contains('ribosome')]
|
|
overlap = pathway.iloc[0]['Overlap'] # e.g., "25/150"
|
|
answer = int(overlap.split('/')[0]) # 25
|
|
```
|
|
|
|
### Pattern 4: Simplify GO terms
|
|
|
|
```python
|
|
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')
|
|
simplified = simplify_go_terms(enr.results, similarity_threshold=0.7)
|
|
answer = len(simplified)
|
|
```
|