Differential Peaks / DMR

Compare two groups for differential peaks or methylation.

Overview

Problem. Where do binding, accessibility or methylation differ?

Learning goals

• DE thinking extends to genomic regions
• Learn to read BED intervals

Figures

Tutorial

Compare two groups of experiments to identify differential peak regions (DPR) or differentially methylated regions (DMR) using the ChIP-Atlas Diff Analysis API.

When to Use This Skill

Use ChIP-Atlas diff analysis when you need to: - Find differential peaks between two conditions (treated vs control, tissue A vs B) - Identify differentially methylated regions between sample groups (Bisulfite-seq) - Compare chromatin accessibility between cell types using ATAC-seq/DNase-seq data - Leverage edgeR-based statistical framework on ChIP-Atlas public experiment data - Avoid raw data downloads — works directly with experiment accession IDs

Key Concept: Submits two groups of experiment IDs to ChIP-Atlas. Server performs edgeR differential analysis (for DPR) or metilene (for DMR), returning BED files with genomic coordinates, logFC, p-values, q-values (FDR), and per-experiment normalized counts.

Installation

pip install pandas requests numpy plotnine plotnine-prism

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

Inputs

Outputs

Clarification Questions

1. Input Files (ASK THIS FIRST):

   • Do you have experiment IDs (SRX/GSM) for two groups to compare?
   • Expected: At least 2 IDs per group (3+ recommended for statistical power)
   • Formats: Plain text (one ID per line), CSV with ID column, or Python list
   • Or use example data? tp53_k562_dmso_vs_daunorubicin (recommended: TP53 K-562 DMSO vs Daunorubicin, n=8 vs n=8) or tp53_molm13_vs_k562 (MOLM-13 vs K-562, n=4 vs n=4, demonstrates sex-chr QC)

2. Analysis parameters:

   • Species/genome? Human hg38 (default), hg19, mouse mm10/mm9, rat rn6, fly, worm, yeast
   • Analysis type? "diffbind" (default, for ChIP/ATAC/DNase-seq) or "dmr" (for Bisulfite-seq)
   • Group labels? Names for each group (e.g., "Treated" vs "Control"). (Auto-set for example data.)

Standard Workflow

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

Step 1 - Load data:

# Option 1: Example data
from scripts.load_example_data import load_example_data
data = load_example_data("tp53_k562_dmso_vs_daunorubicin")
experiments_a = data['experiments_a']
experiments_b = data['experiments_b']

# Option 2: Your own experiment IDs
# from scripts.load_user_data import load_user_data
# experiments_a = load_user_data("group_a.txt")
# experiments_b = load_user_data("group_b.txt")

Step 2 - Run diff analysis:

from scripts.run_diff_workflow import run_diff_workflow

results = run_diff_workflow(
    experiments_a=experiments_a,
    experiments_b=experiments_b,
    genome=data.get('genome', 'hg38'),
    analysis_type=data.get('analysis_type', 'diffbind'),
    description_a=data.get('description_a', 'group_A'),
    description_b=data.get('description_b', 'group_B'),
    output_dir="diff_analysis_results"
)

Step 3 - Generate visualizations:

from scripts.generate_all_plots import generate_all_plots
generate_all_plots(results, output_dir="diff_analysis_results")

Step 4 - Export results:

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

Common Issues

Interpretation Guidelines

logFC: Log2 fold change — positive = enriched in Group A, negative = enriched in Group B Q-value (FDR): Benjamini-Hochberg corrected p-value — regions with FDR < 0.05 are significant Counts: TMM-normalized read counts per experiment (comma-separated in raw BED) Region Size: 100-500bp (TF sites), 500-2000bp (enhancers), >2000bp (broad marks)

Nearest gene annotations: Top regions in the summary report and top50 CSV are annotated with overlapping or nearby genes via UCSC REST API. Shows "GENE" for direct overlap or "GENE (+2.1kb)" for nearest gene within 5kb. Annotation is optional and skipped gracefully if the UCSC API is unavailable.

Reporting top hits: When summarizing top differential regions, present them in strict Q-value ranking order without skipping unannotated regions. Regions lacking gene annotations are still statistically significant. If you curate a subset (e.g., only gene-annotated hits), explicitly label the table as "selected highlights" — never imply it is a strict significance ranking. Accuracy: When quoting specific statistics (logFC, Q-values, coordinates) in your summary, copy values directly from summary_report.md or the CSV files — do NOT re-derive, round, or paraphrase numerical values from memory, as this introduces transcription errors.

Gene annotations: Nearest gene annotations indicate genomic proximity, not functional regulation. Describe associations factually (e.g., "a differential region near PCNA") without asserting direct regulatory relationships. Do not use individual gene annotations as validation of the analysis — gene proximity does not prove functional regulation.

Experimental design: edgeR treats all experiments within a group as biological replicates. If experiments within a group differ in important ways (e.g., different genotypes, tissues, or treatments), state this clearly — the differential analysis captures the average effect between groups, and within-group heterogeneity is treated as noise. When summarizing results to the user, design caveats MUST be stated prominently — not buried at the end. If within-group experiments are not true biological replicates, this is the single most important limitation to communicate.

MA plot interpretation: The MA plot shows mean expression (average log2 counts, x-axis) vs log fold change (y-axis). Normal pattern: funnel/fan shape narrowing at high counts (more power to detect small changes). Red flags: asymmetric significant hits clustering at low counts in one direction (potential normalization artifact), horizontal bands of significant hits (batch effects), or all significant hits at very low counts (noise-driven). CRITICAL: Do NOT claim MA plot asymmetry in your summary unless the automated QC check flagged it. The quantitative threshold is >1 log2 unit difference in median mean count between directions. Visual impressions of asymmetry below this threshold are not meaningful — small differences (e.g., 0.2 log2 units) are negligible and do not support claims of noise-driven signal. If no MA asymmetry QC warning appears, describe the MA plot as showing a normal pattern.

Biological interpretation of top hits: When discussing top differential regions near known genes, briefly note the gene's known function for context but do NOT claim the differential peak proves direct regulation. Genomic proximity does not establish causation. For exploratory analyses, suggest formal pathway/GO enrichment analysis as a next step to identify systematic patterns. Distinguish between confirmatory hits (regions near genes expected to differ) and exploratory hits (unexpected associations requiring validation). When the target factor is well-characterized (e.g., TP53, CTCF, ESR1, NF-kB), note whether canonical target genes appear in the top differential regions — the absence of expected targets is itself a notable finding that may indicate the differential signal is dominated by artifacts or indirect effects. After excluding artifact regions, comment on the effect sizes (logFC) of remaining biologically interpretable hits — modest effect sizes (|logFC| < 1.5) combined with annotations near pseudogenes or lncRNAs suggest limited genuine differential binding.

Citations: When writing your final summary, cite the data source (e.g., GSE accession) and key tools: ChIP-Atlas (Zou et al. 2024), edgeR (Robinson et al. 2010) for DPR, or metilene (Jühling et al. 2016) for DMR. Full references are in the References section below.

Caveats: Results depend on available public data quality. Always check QC warnings in the summary report. Validate key regions with orthogonal methods. See references/output_format.md for detailed output format.

Suggested Next Steps

1. Visualize in IGV using the .igv.xml session file or .igv.bed file
2. Annotate regions with nearby genes using bedtools or GREAT
3. Motif analysis on differential peaks (HOMER, MEME-ChIP)
4. Integrate with expression data to link peaks to gene regulation

Related Skills

• chip-atlas-peak-enrichment - Enrichment of ChIP-seq peaks near gene lists
• bulk-rnaseq-counts-to-de-deseq2 - Differential expression to identify genes for downstream analysis

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
• Robinson et al. (2010). edgeR: a Bioconductor package for differential expression analysis. Bioinformatics 26(1):139-140.
• Jühling et al. (2016). metilene: fast and sensitive calling of differentially methylated regions. Genome Research 26(2):256-262.
• ChIP-Atlas: https://chip-atlas.org
• API documentation: See references/chipatlas_diff_api_format.md
• Statistical methods: See references/diff_analysis_methods.md

Code preview

scripts/__init__.py

# chip-atlas-diff-analysis scripts package

scripts/annotate_genes.py

"""
Annotate differential regions with nearest/overlapping genes via UCSC REST API.

Uses the UCSC Genome Browser REST API to find genes overlapping or near
each genomic region. Requires internet (same as ChIP-Atlas API dependency).

Graceful fallback: if UCSC API is unavailable, returns empty annotations
without raising errors.
"""

import pandas as pd
import requests

UCSC_API_URL = "https://api.genome.ucsc.edu/getData/track"
FLANK_BP = 5000


def annotate_nearest_genes(df, genome="hg38", max_regions=50):
    """
    Annotate regions with overlapping/nearest gene symbols via UCSC API.

    Parameters
    ----------
    df : pd.DataFrame
        DataFrame with chrom, chromStart, chromEnd columns.
        Should be pre-sorted by significance (top regions first).
    genome : str
        UCSC genome assembly (hg38, mm10, etc.)
    max_regions : int
        Maximum number of regions to annotate (default: 50).

    Returns
    -------
    pd.Series
        Gene annotations indexed to match df. Each entry is a string
        like "SBF2" (overlapping) or "SBF2 (+2.1kb)" (nearest within flank).
        Empty string for unannotated regions.
    """
    annotations = pd.Series("", index=df.index)

    regions_to_query = df.head(max_regions)
    if len(regions_to_query) == 0:
        return annotations

    print(f"   Annotating top {len(regions_to_query)} regions with nearest genes (UCSC API)...")

    n_annotated = 0
    n_failed = 0

    for idx, row in regions_to_query.iterrows():
        try:
            gene_str = _query_genes_for_region(
                row["chrom"],
                int(row["chromStart"]),
                int(row["chromEnd"]),
                genome=genome,
            )
            if gene_str:
                annotations.at[idx] = gene_str
                n_annotated += 1
        except Exception:
            n_failed += 1
            if n_failed >= 3:
                print(f"   Warning: UCSC API unavailable ({n_failed} failures). "
                      f"Skipping remaining gene annotations.")
                break

    print(f"   Annotated {n_annotated}/{len(regions_to_query)} regions with gene symbols")
    return annotations


def _query_genes_for_region(chrom, start, end, genome="hg38"):
    """
    Query UCSC REST API for genes overlapping or near a region.

    Returns
    -------
    str
        Gene symbol(s) or empty string.
        Direct overlap: "GENE1, GENE2"

scripts/export_all.py

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

Exports:
1. analysis_object.pkl - Complete results for downstream use
2. diff_regions_all.csv - All differential regions
3. diff_regions_significant.csv - Significant regions (FDR < 0.05)
4. diff_regions_top50.csv - Top 50 by significance
5. summary_report.md - Human-readable summary
"""

import os
import pickle
from datetime import datetime

try:
    try:
        from scripts.annotate_genes import annotate_nearest_genes
    except ImportError:
        from annotate_genes import annotate_nearest_genes
    _HAS_GENE_ANNOTATION = True
except ImportError:
    _HAS_GENE_ANNOTATION = False


def export_all(results, output_dir="diff_analysis_results", qvalue_threshold=0.05):
    """
    Export all ChIP-Atlas diff analysis results with pickle object.

    Parameters
    ----------
    results : dict
        Results from run_diff_workflow()
    output_dir : str
        Output directory
    qvalue_threshold : float
        Maximum Q-value (FDR) for "significant" regions (default: 0.05)

    Verification
    ------------
    Prints "=== Export Complete ===" when done
    """
    os.makedirs(output_dir, exist_ok=True)

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

    df = results["diff_regions"]
    params = results["parameters"]

    qc_warnings = results.get("qc_warnings", [])
    unfiltered_df = results.get("diff_regions_unfiltered", df)

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

    n_a = (df["direction"] == "A_enriched").sum() if len(df) > 0 else 0
    n_b = (df["direction"] == "B_enriched").sum() if len(df) > 0 else 0

    analysis_object = {
        "diff_regions": df,
        "diff_regions_unfiltered": unfiltered_df,
        "qc_warnings": qc_warnings,
        "experiments_a": results["experiments_a"],
        "experiments_b": results["experiments_b"],
        "raw_files": results.get("raw_files", {}),
        "log_content": results.get("log_content", ""),
        "parameters": params,
        "n_regions_total": len(df),
        "n_regions_unfiltered": len(unfiltered_df),
        "n_regions_a_enriched": int(n_a),
        "n_regions_b_enriched": int(n_b),
        "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)

Companion files