Peak Enrichment (genes → TF)

Given a gene set, which TFs are enriched over it.

Overview

Problem. Which TFs co-regulate the genes in my list?

Learning goals

• The Fisher + BH pattern recurs everywhere
• Reverse-infer regulation from public data

Figures

Tutorial

Find ChIP-seq peak enrichment near your genes using the official ChIP-Atlas Enrichment Analysis API.

When to Use This Skill

Use ChIP-Atlas peak enrichment when you need to: - Identify transcription factors binding near your genes from DE analysis or pathway results - Discover chromatin regulators (TFs, histone modifications, chromatin remodelers) enriched near your gene set - Validate regulatory relationships between factors and target genes using public ChIP-seq data - Find cell-type-specific regulators by filtering to specific cell classes - Query 433,000+ ChIP-seq experiments via the official API without manual downloads

Key Concept: Submits your gene list to the ChIP-Atlas API, which performs Fisher's exact test enrichment analysis against all public ChIP-seq experiments. Returns fold enrichment, P-values, and BH-corrected Q-values.

Installation

pip install pandas requests numpy plotnine plotnine-prism

System requirements: Internet connection (API calls to ChIP-Atlas and Ensembl)

Inputs

Outputs

Clarification Questions

1. Input Files (ASK THIS FIRST):

   • Do you have a gene list to analyze?
   • Expected formats: Plain text (one gene per line), CSV with gene column, or DE results file
   • Or use example data? tp53_targets (5 genes, fast test) or immune_response (20 genes)

2. Analysis parameters:

   • Species/genome? Human hg38 (default), hg19, mouse mm10/mm9, rat rn6, fly, worm, yeast
   • Experiment type? "TFs and others" (default), "Histone", "ATAC-Seq"
   • Cell type class? "All cell types" (default), or specific: "Blood", "Neural", "Breast", etc.
   • Peak threshold? 50 (default, balanced), 100 (stringent), 200 (very stringent)
   • TSS window? 5000bp up/down (default), or custom distance_up/distance_down

Standard Workflow

Step 1 - Load data:

# Option 1: Example data
from scripts.load_example_data import load_example_data
data = load_example_data("tp53_targets")
gene_list = data['genes']

# Option 2: Your own genes
# from scripts.load_user_data import load_user_data
# gene_list = load_user_data("my_genes.txt")

Step 2 - Run enrichment analysis:

from scripts.run_enrichment_workflow import run_enrichment_workflow

results = run_enrichment_workflow(
    gene_list=gene_list,
    genome="hg38",
    antigen_class="TFs and others",
    cell_class="All cell types",
    threshold=50,
    output_dir="chipatlas_results"
)

Step 3 - Generate visualizations:

from scripts.generate_all_plots import generate_all_plots
generate_all_plots(results, output_dir="chipatlas_results", top_n=15)

Step 4 - Export results:

from scripts.export_all import export_all
export_all(results, output_dir="chipatlas_results")

Common Issues

Interpretation Guidelines

Q-value (BH-corrected, primary ranking): <0.001 highly significant, <0.01 significant, <0.05 genome-wide significant Fold Enrichment: >10x very strong, 5-10x strong, 2-5x moderate, <2x weak, >100,000x likely sentinel value (zero background overlap) Overlap Rate: >50% core regulon, 20-50% key targets, <20% subset/indirect Threshold: Controls MACS2 peak-calling stringency. Default 50 ≈ p<1e-5. See references/peak_thresholds.md.

Suggested Next Steps

1. Threshold sensitivity check — Re-run at threshold=100 or 200 to test whether top factors remain significant at more stringent peak-calling cutoffs. Offer this to the user as a robustness check (e.g., "Would you like me to re-run at threshold=100 to confirm these findings are robust?").
2. Validate top factors with literature, expression correlation, perturbation data
3. Cell-type-specific analysis with cell_class matching your experimental system
4. Motif analysis of promoter regions for top factor binding motifs
5. Regulatory network construction with top factors and target genes

Related Skills

• gene-correlation-archs4 - Co-expression across 600K RNA-seq samples
• grn-pyscenic - Gene regulatory networks from single-cell data

References

• Zou et al. (2024). ChIP-Atlas 3.0: a data-mining suite to explore chromosome architecture. Nucleic Acids Research. doi:10.1093/nar/gkad884
• Zou et al. (2022). ChIP-Atlas 2021 update. Nucleic Acids Research. doi:10.1093/nar/gkab933
• Oki et al. (2018). ChIP-Atlas: a data-mining suite. EMBO Reports 19(12):e46255. doi:10.15252/embr.201846255
• ChIP-Atlas: https://chip-atlas.org
• API documentation: See references/chipatlas_metadata_format.md
• Enrichment statistics: See references/enrichment_statistics.md

Code preview

scripts/__init__.py

# chip-atlas-peak-enrichment scripts package

scripts/convert_genes_to_regions.py

"""
Convert gene symbols to promoter regions using Ensembl REST API.

Supports batch lookup (POST /lookup/symbol) for all genes in one request,
with sequential fallback for robustness.

Includes rate limiting and error handling for robust API queries.
"""

import requests
import time


# Common gene symbol renames (retired → current HGNC symbol)
# Source: HGNC (genenames.org) approved symbol changes
GENE_ALIASES = {
    # Interleukins/chemokines
    'IL8': 'CXCL8',
    'IL8A': 'CXCL8',
    'GRO1': 'CXCL1',
    'GRO2': 'CXCL2',
    'GRO3': 'CXCL3',
    'SDF1': 'CXCL12',
    'MIG': 'CXCL9',
    'IP10': 'CXCL10',
    'MCP1': 'CCL2',
    'MIP1A': 'CCL3',
    'MIP1B': 'CCL4',
    'RANTES': 'CCL5',
    # Apoptosis/signaling
    'TNFSF6': 'FASLG',
    'APO1': 'FAS',
    'TRAIL': 'TNFSF10',
    # Epigenetic regulators
    'MLL': 'KMT2A',
    'MLL2': 'KMT2D',
    'MLL3': 'KMT2C',
    # Other common renames
    'P53': 'TP53',
    'ERK1': 'MAPK3',
    'ERK2': 'MAPK1',
    'P38': 'MAPK14',
}


def convert_genes_to_regions(gene_list, genome="hg38", upstream=2000, downstream=500):
    """
    Convert gene symbols to promoter regions via Ensembl API.

    Parameters:
    -----------
    gene_list : list of str
        Gene symbols (e.g., ['TP53', 'MYC'])
    genome : str
        Genome assembly ('hg38', 'hg19', 'mm10', 'mm9')
    upstream : int
        Bases upstream of TSS
    downstream : int
        Bases downstream of TSS

    Returns:
    --------
    tuple: (successful_regions, failed_genes)
        successful_regions: list of tuples (chr, start, end, gene, strand)
        failed_genes: list of str (genes that couldn't be mapped)

    Notes:
    ------
    - Rate limited to 0.5s between requests
    - Retries on timeout (max 3 attempts)
    - Skips genes not found in Ensembl
    """

    # Map genome assembly prefix to Ensembl species name
    _GENOME_TO_SPECIES = {
        'hg': 'human', 'mm': 'mouse', 'rn': 'rat',
        'dm': 'drosophila_melanogaster', 'ce': 'caenorhabditis_elegans',
        'sacCer': 'saccharomyces_cerevisiae',
    }
    prefix = next((k for k in _GENOME_TO_SPECIES if genome.startswith(k)), None)

scripts/export_all.py

"""
Export all ChIP-Atlas enrichment results (Step 4).

Exports:
1. analysis_object.pkl - Complete results for downstream use
2. enrichment_results_all.csv - All experiments analyzed
3. enrichment_results_significant.csv - Significant results (q < 0.05)
4. enrichment_results_top20.csv - Top 20 by significance (q-value)
5. failed_genes.txt - Genes that couldn't be mapped
6. summary_report.md - Human-readable summary
"""

import os
import pickle
from datetime import datetime


def export_all(results, output_dir="chipatlas_results"):
    """
    Export all ChIP-Atlas enrichment results with pickle object.

    Parameters:
    -----------
    results : dict
        Results from run_enrichment_workflow()
    output_dir : str
        Output directory

    Returns:
    --------
    None
        Prints export messages and verification

    Verification:
    -------------
    Prints "=== Export Complete ===" when done
    """

    os.makedirs(output_dir, exist_ok=True)

    print("\n" + "="*70)
    print("EXPORTING RESULTS")
    print("="*70 + "\n")

    # 1. Save analysis object as pickle (CRITICAL for downstream skills)
    print("1. Saving analysis objects for downstream use...")

    analysis_object = {
        'enrichment_results': results['enrichment_results'],
        'input_genes': results['input_genes'],
        'input_regions': results['input_regions'],
        'failed_genes': results.get('failed_genes', []),
        'metadata': results['metadata'],
        'parameters': results['parameters'],
        'n_genes': len(results['input_genes']),
        'n_regions': len(results['input_regions']),
        'n_experiments': len(results['enrichment_results']),
        'n_significant': len(results['enrichment_results'][
            results['enrichment_results']['q_value'] < 0.05
        ]),
        'timestamp': datetime.now().isoformat()
    }

    pickle_path = os.path.join(output_dir, "analysis_object.pkl")
    with open(pickle_path, 'wb') as f:
        pickle.dump(analysis_object, f)

    print(f"   Saved: {pickle_path}")
    print(f"   (Load with: import pickle; obj = pickle.load(open('analysis_object.pkl', 'rb')))")

    # 2. Export all results to CSV
    print("\n2. Exporting enrichment results...")

    all_results_path = os.path.join(output_dir, "enrichment_results_all.csv")
    results['enrichment_results'].to_csv(all_results_path, index=False)
    print(f"   Saved: {all_results_path} ({len(results['enrichment_results'])} experiments)")

    # 3. Export significant results (q < 0.05, BH-corrected)
    significant = results['enrichment_results'][
        results['enrichment_results']['q_value'] < 0.05

Companion files