bio-chipseq-super-enhancers

What this skill does

## Version Compatibility

Reference examples tested with: ROSE (stjude/ROSE, 2018+), ROSE2 (linlabbcm/rose2, 2021+), LILY (BoevaLab/LILY, 2020+), HOMER 4.11+, samtools 1.19+, bedtools 2.31+, GenomicRanges 1.54+.

ROSE is unmaintained Python 2; ROSE2 is the Python 3 port with the same algorithm. LILY (Boeva 2017) is a refactored implementation with input-control background subtraction for low-quality H3K27ac data.

# Super-Enhancer Calling

**"Identify super-enhancers driving cell identity / cancer biology"** -> Stitch nearby active enhancer peaks (H3K27ac, MED1, or BRD4) within a stitching window, exclude proximal-promoter signal, rank by total signal, find the hockey-stick inflection point where signal sharply increases, and classify all stitched regions above the inflection as super-enhancers.

- CLI (ROSE / ROSE2): `python ROSE_main.py -g HG38 -i peaks.gff -r h3k27ac.bam -c input.bam -s 12500 -t 2500`
- CLI (HOMER): `findPeaks tag_dir/ -style super -i input_tag_dir/`
- CLI (LILY): variant with input-control background subtraction
- R (custom hockey-stick): rank enhancers by signal, find tangent-line inflection

The SE concept (Whyte 2013) is a thresholding heuristic on a continuous signal distribution (Pott & Lieb 2015 *Nat Genet*), not a categorical biological category. Functional CRISPR-tiling at SE loci (Hnisz 2017; Dukler 2017) shows only 1-3 constituent elements per SE are essential; the "SE" label is a useful operational definition for BET-inhibitor responsiveness and cell-identity gene regulation, not an absolute biological property.

## Marker Choice: H3K27ac vs MED1 vs BRD4

| Marker | Captures | When to prefer |
|--------|----------|----------------|
| **H3K27ac** | Active regulatory elements broadly | Most widely available; standard for SE definition since Whyte 2013 |
| **MED1** | Mediator complex accumulation (the defining biology) | Direct readout of SE; less common antibody; lower signal-to-noise |
| **BRD4** | BET cofactor accumulation | Most predictive of BET-inhibitor responsiveness; clinical relevance |
| **H3K27ac + MED1 intersection** | High-confidence SE | Gold standard if both available |
| **dELS from ENCODE cCREs** | Cell-type-agnostic distal enhancer registry | Cross-reference; not SE-specific by itself |

**Operational rule:** H3K27ac for discovery; MED1 or BRD4 ChIP for functional / therapeutic claims. SE called on H3K27ac alone may not respond to BET inhibitors; SE called on BRD4 will.

## Algorithmic Taxonomy

| Tool | Method | Strength | Fails when |
|------|--------|----------|------------|
| **ROSE** (Whyte 2013) | Stitch within 12.5 kb, exclude ±2.5 kb of TSS, rank by signal, hockey-stick inflection | Original; widely cited; canonical reference | Python 2 only; unmaintained; ROSE_main.py crashes on Python 3 |
| **ROSE2** (Lin 2016; linlabbcm) | Same algorithm, Python 3 port | Maintained; identical output to ROSE | None vs ROSE |
| **LILY** (Boeva 2017) | ROSE-like with input-control background subtraction | Works on lower-quality H3K27ac data; subtracts input | Adds complexity; less validated; specific to neuroblastoma/glioma in original paper |
| **HOMER `-style super`** | Native ROSE-like in HOMER framework; stitching without TSS exclusion | Integrated with HOMER workflow | Different stitching defaults; not directly comparable to ROSE counts |
| **Custom hockey-stick (R)** | Generic rank-by-signal + tangent inflection | Flexible; works on any signal definition | Reinvents algorithm; verify against ROSE on known dataset |

**Most papers use ROSE/ROSE2 with default stitching (12.5 kb) and TSS exclusion (2.5 kb).** This is the de facto standard for cross-paper comparison. HOMER's `-style super` produces different counts and is not directly comparable.

## Decision Tree: SE Calling Workflow

| Scenario | Recommended pipeline |
|----------|----------------------|
| Standard SE discovery, H3K27ac available | ROSE2 with default `-s 12500 -t 2500`; input control for subtraction |
| Predict BET-inhibitor response | BRD4 ChIP -> ROSE2 (or H3K27ac SE intersected with BRD4 peaks) |
| Compare SE between conditions (drug treatment) | ROSE2 per condition + spike-in normalization (HDACi/BETi/EZH2i need ChIP-Rx) |
| Build core regulatory circuitry | ROSE2 + Saint-Andre 2016 algorithm: identify TFs encoded by SE that bind own SE + cross-bind other SE-encoded TFs |
| Low-quality H3K27ac (low FRiP) | LILY with input subtraction |
| Compare with ENCODE dELS atlas | ROSE2 + intersect with ENCODE cCRE dELS BED |
| Differential SE between conditions | ROSE2 per condition + signal-quantitative differential (DiffBind on SE regions) |

## ROSE / ROSE2 Workflow

**Goal:** Identify super-enhancers by stitching nearby active enhancer peaks within a stitching distance and ranking by total signal.

**Approach:** Convert peaks to GFF, exclude promoter-proximal peaks via `-t` (TSS exclusion window), stitch enhancers within `-s` (default 12.5 kb), rank by total H3K27ac (or MED1/BRD4) signal, find the hockey-stick inflection point, classify regions above as super-enhancers.

```bash
# Install ROSE2 (Python 3 port; unmaintained ROSE Py2 not recommended)
git clone https://github.com/linlabbcm/rose2.git
pip install ./rose2

# Convert peaks BED to GFF (ROSE requires GFF input)
awk 'BEGIN{OFS="\t"} {print $1,"peaks","enhancer",$2,$3,".",$6,".","ID="NR}' \
    peaks.narrowPeak > peaks.gff

# Filter promoter peaks before SE calling (within 2.5 kb of TSS)
# ROSE handles this via -t flag; preferable to pre-filter for clarity
bedtools intersect -a peaks.narrowPeak -b promoters_2kb.bed -v > enhancer_peaks.bed

# Run ROSE2 with input control
rose2 -g HG38 -i peaks.gff \
    -r h3k27ac.bam -c input.bam \
    -o rose_output/ \
    -s 12500 \
    -t 2500
```

ROSE2 outputs:
- `*_AllEnhancers.table.txt` — all stitched enhancer regions ranked by signal
- `*_SuperEnhancers.table.txt` — SE only (above hockey-stick inflection)
- `*_Enhancers_withSuper.bed` — BED with SE / TE classification
- `*_Plot_points.png` — hockey-stick plot

## Cross-Condition SE Comparison

This is the analysis most often done wrong. SE calling thresholds depend on absolute signal, so any global shift (HDACi, BETi, EZH2i) confounds direct SE-count comparison.

**Wrong approach:** Call ROSE2 on condition A and condition B separately, intersect SE BEDs, report "gained/lost SE."

**Right approach:**
1. Spike-in normalize signal between conditions (see chip-seq/spike-in-normalization)
2. Build a union SE set from both conditions
3. Quantify signal at union SE regions per condition (DiffBind on the union)
4. Apply differential testing with appropriate normalization (background-bin TMM or spike-in)

```r
library(DiffBind)
# Union of SE BED files from condition A and B
union_se <- rtracklayer::import('union_SE.bed')
# Run DiffBind quantification on this region set with spike-in normalization
```

For BET-inhibitor experiments: the biology IS that all SE decrease globally; spike-in is mandatory.

## Core Regulatory Circuitry (Saint-André 2016)

The CRC algorithm identifies master TF networks from SE annotations:

1. List all TFs encoded by SE-associated genes
2. For each such TF, check if its motif appears in its own SE (auto-regulation)
3. Build a graph where TFs encoded by SE-A bind to motifs in SE-B
4. Identify highly-interconnected sub-networks (CRC)

```bash
# Install CRC pipeline (linlabbcm CRC)
git clone https://github.com/linlabbcm/CRC2.git

# Requires: SE BED, TF motif annotations, gene-SE mapping
python CRC2/crc.py -e SE_table.txt -g genes.gtf -b h3k27ac.bam
```

CRC outputs the connected components of the regulatory network. Master TFs typically appear in the largest component with high out-degree.

## ENCODE dELS Cross-Reference

ENCODE distal Enhancer-Like Signatures (dELS) are the cell-type-agnostic regulatory atlas (see chip-seq/peak-annotation). Cross-referencing SE against dELS:

- Validates SE constituents are at canonical regulatory elements
- Identifies SE constituents NOT in t