/histone-aggregation

name: histone-aggregation
description: Build comprehensive histone mark maps by aggregating narrowPeak data across multiple ENCODE experiments, donors, and labs. Use when the user wants to answer "where is this histone mark present in my tissue?" by combining peak calls from multiple studies into a union peak set with confidence annotations. Handles cross-lab batch effects, broad vs narrow marks, and ENCODE blocklist filtering.

Aggregate Histone ChIP-seq Peaks Across Studies

When to Use

• User wants to combine histone ChIP-seq peaks across multiple ENCODE experiments for a tissue or cell type
• User asks "where is H3K27ac in pancreas?" or "build a histone mark map for liver"
• User needs a union peak set from multiple donors, labs, or replicates
• User wants to create a consensus binding map from multiple ChIP-seq datasets
• Example queries: "aggregate all H3K4me3 peaks in brain", "combine histone marks across donors", "build enhancer map from H3K27ac data"

Build a comprehensive map of histone mark binding for a tissue/cell type by merging narrowPeak files from multiple ENCODE experiments into a union peak set.

Scientific Rationale

**The question**: "Does my tissue have this histone mark, and at what genomic locations?"

This is a **detection/cataloging** question, not a differential one. Once a histone mark passes noise thresholds (ENCODE IDR, quality metrics), detection is binary — the mark is either bound or not. If detected in one donor but not another, that region is still a real binding site. Individual variation and technical differences (lab, depth, antibody lot) explain *absence*, not that *presence* is spurious.

**Therefore: we want the UNION of all detections, not a consensus.**

Literature Support

• **ChIP-Atlas** (Oki et al. 2018, EMBO Reports, 597 citations): Integrated >70,000 public ChIP-seq datasets using union of all peak calls
• **ENCODE Phase 3** (Gorkin et al. 2020, Nature, 301 citations): Created unified chromatin state annotations by integrating all peaks across 1,128 ChIP-seq experiments
• **ENCODE Blacklist** (Amemiya et al. 2019, Scientific Reports, 1,372 citations): Defined the comprehensive set of problematic genomic regions to filter from all functional genomics analyses. Essential quality step. [DOI](https://doi.org/10.1038/s41598-019-45839-z)
• **Perna et al. 2024** (BMC Genomics): Found top 25% signalValue peaks most consistent across different processing pipelines — use as per-sample noise filter
• **ChIP-R** (Newell et al. 2020, 26 citations): Rank-product method for combining peaks from multiple replicates without BAMs, works directly on narrowPeak files
• **MSPC** (Jalili et al. 2021, BMC Bioinformatics): Rescues weak-but-real binding sites that IDR discards by exploiting replicates to lower calling thresholds — more sensitive alternative for union-based approaches
• **Hecht et al. 2023** (PLoS Comp Bio): Probability-of-Being-Signal (PBS) approach for cross-dataset comparison with differing read depths

Step 1: Find All Available Experiments

Search for all histone ChIP-seq data for the target mark and tissue:

encode_search_experiments(
    assay_title="Histone ChIP-seq",
    target="H3K4me1",        # or H3K27ac, H3K4me3, H3K27me3, etc.
    organ="pancreas",         # user's tissue of interest
    biosample_type="tissue",  # or "cell line", "primary cell"
    limit=100
)

Present a summary table to the user showing:

• Number of experiments found
• Labs represented
• Number of unique donors/biosamples
• Any audit flags

Use `encode_get_facets` first if unsure what's available:

encode_get_facets(assay_title="Histone ChIP-seq", organ="pancreas")

Step 2: Quality-Gate Each Experiment

For each experiment, check quality before including:

encode_get_experiment(accession="ENCSR...")

Include if:

• Audit status: no ERROR flags (WARNING is acceptable)
• Has IDR thresholded peaks (passed replicate concordance)
• Sequencing depth meets ENCODE standards (10M+ for narrow marks, 20M+ for broad)

Exclude if:

• ERROR audit flags
• Only pseudoreplicated peaks (no IDR = did not pass reproducibility)
• Known antibody issues (check audit details)

Track all included experiments:

encode_track_experiment(accession="ENCSR...")

Step 3: Download IDR Thresholded NarrowPeak Files

For each passing experiment, get the peak files:

encode_list_files(
    experiment_accession="ENCSR...",
    file_format="bed",
    output_type="IDR thresholded peaks",
    assembly="GRCh38"
)

**File selection priority:** 1. **IDR thresholded peaks** (gold standard — passed replicate concordance) 2. **Optimal IDR peaks** (pooled replicates — most complete set) 3. **Replicated peaks** (alternative peak caller output)

Prefer `preferred_default=True` files when available.

Download all selected files:

encode_download_files(
    file_accessions=["ENCFF...", "ENCFF...", ...],
    download_dir="/path/to/data/narrowpeaks",
    organize_by="flat"
)

Step 4: Per-Sample Noise Filtering

**IMPORTANT**: Filter BEFORE merging, not after.

4a. ENCODE Blocklist Filtering (Amemiya et al. 2019)

Remove artifact-prone regions (centromeres, telomeres, rDNA repeats, satellite repeats):

# Download ENCODE blocklist for GRCh38 from:
# https://github.com/Boyle-Lab/Blacklist/blob/master/lists/hg38-blacklist.v2.bed.gz
# For mm10: https://github.com/Boyle-Lab/Blacklist/blob/master/lists/mm10-blacklist.v2.bed.gz
bedtools intersect -a sample.narrowPeak -b hg38-blacklist.v2.bed -v > sample.filtered.narrowPeak

4b. SignalValue Filtering (Perna et al. 2024)

Filter each sample's peaks to retain those above the 25th percentile signalValue (column 7 in narrowPeak). The top 75% of peaks by signalValue are the most reliable across processing pipelines:

# Calculate the 25th percentile of the signalValue DISTRIBUTION for this sample
# (This is a true quantile, not 25% of the range)
TOTAL=$(wc -l < sample.filtered.