bio-chipseq-super-enhancers

## Version Compatibility

Reference examples tested with: GenomicRanges 1.54+, bedtools 2.31+, ggplot2 3.5+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Super-Enhancer Calling

**"Identify super-enhancers from H3K27ac ChIP-seq"** → Stitch nearby enhancer peaks and rank by signal to find large regulatory domains controlling cell identity genes.
- CLI: `ROSE_main.py -g hg38 -i peaks.gff -r chip.bam -c input.bam`

Identify super-enhancers (SEs) - large clusters of enhancers that control cell identity genes.

## Background

Super-enhancers are:
- Large clusters of enhancer regions
- Marked by H3K27ac, Med1, BRD4
- Control cell identity genes
- Often altered in disease/cancer

## ROSE (Rank Ordering of Super-Enhancers)

### Installation

```bash
git clone https://github.com/stjude/ROSE.git
cd ROSE
# Requires samtools, R, bedtools
```

### Input Requirements

1. **BAM file** - H3K27ac ChIP-seq aligned reads
2. **Peak file** - Called peaks (BED or GFF)
3. **Genome annotation** - TSS annotations

### Run ROSE

**Goal:** Identify super-enhancers by stitching nearby enhancer peaks and ranking by H3K27ac signal.

**Approach:** Run ROSE_main.py with a GFF peak file, ChIP-seq BAM, and optional input control to stitch enhancers within 12.5 kb, rank by signal, and identify the inflection point separating super-enhancers from typical enhancers.

```bash
# Basic usage
python ROSE_main.py \
    -g HG38 \
    -i peaks.gff \
    -r h3k27ac.bam \
    -o output_dir \
    -s 12500 \
    -t 2500

# With control/input
python ROSE_main.py \
    -g HG38 \
    -i peaks.gff \
    -r h3k27ac.bam \
    -c input.bam \
    -o output_dir
```

### Key Parameters

| Parameter | Description | Default |
|-----------|-------------|---------|
| `-s` | Stitching distance | 12500 bp |
| `-t` | TSS exclusion | 2500 bp |
| `-c` | Control BAM | None |

### Output Files

```
output_dir/
├── *_AllEnhancers.table.txt        # All enhancer regions
├── *_SuperEnhancers.table.txt      # Super-enhancers only
├── *_Enhancers_withSuper.bed       # BED with SE annotation
└── *_Plot_points.png               # Hockey stick plot
```

## Prepare Input Files

### Convert BED to GFF

```bash
# ROSE requires GFF format for peaks
awk 'BEGIN{OFS="\t"} {print $1,"peaks","enhancer",$2,$3,".",$6,".","ID="NR}' \
    peaks.bed > peaks.gff
```

### Filter Peaks for Enhancers

```bash
# Remove promoter peaks (within 2.5kb of TSS)
bedtools intersect -a peaks.bed -b promoters.bed -v > enhancer_peaks.bed
```

## Alternative: HOMER Super-Enhancers

```bash
# Call super-enhancers with HOMER
findPeaks tag_dir/ -style super -o auto

# Or from existing peaks
findPeaks tag_dir/ -style super -i input_tag_dir/ \
    -typical typical_enhancers.txt \
    -superSlope -1000 \
    > super_enhancers.txt
```

## Alternative: SEanalysis

```bash
# R-based analysis
Rscript << 'EOF'
library(SEanalysis)

# Load H3K27ac signal at enhancers
signal <- read.table('enhancer_signal.txt', header=TRUE)

# Rank and identify super-enhancers
se_result <- identifySE(signal$signal, method='ROSE')

# Get super-enhancer IDs
super_enhancers <- signal$id[se_result$is_super]
write.table(super_enhancers, 'super_enhancers.txt', quote=FALSE, row.names=FALSE)
EOF
```

## Custom Hockey Stick Analysis (R)

**Goal:** Classify enhancers as super-enhancers vs typical using a custom hockey stick plot and inflection-point detection.

**Approach:** Rank enhancers by normalized signal, compute the slope at each point, find where the tangent exceeds 1 (inflection point), and classify all enhancers above the inflection as super-enhancers.

```r
library(ggplot2)

# Load enhancer signal data
enhancers <- read.table('enhancer_signal.txt', header=TRUE)

# Rank by signal
enhancers <- enhancers[order(enhancers$signal), ]
enhancers$rank <- 1:nrow(enhancers)

# Find inflection point (tangent = 1)
# Normalize ranks and signal to 0-1
enhancers$rank_norm <- enhancers$rank / max(enhancers$rank)
enhancers$signal_norm <- enhancers$signal / max(enhancers$signal)

# Calculate slope at each point
n <- nrow(enhancers)
slopes <- diff(enhancers$signal_norm) / diff(enhancers$rank_norm)
inflection <- which(slopes > 1)[1]

# Classify
enhancers$type <- ifelse(enhancers$rank >= inflection, 'Super-Enhancer', 'Typical')

# Plot
ggplot(enhancers, aes(rank, signal, color = type)) +
    geom_point(size = 0.5) +
    scale_color_manual(values = c('Super-Enhancer' = 'red', 'Typical' = 'grey60')) +
    geom_vline(xintercept = inflection, linetype = 'dashed') +
    labs(x = 'Enhancer Rank', y = 'H3K27ac Signal', title = 'Super-Enhancer Identification') +
    theme_bw()

ggsave('hockey_stick_plot.pdf', width = 8, height = 6)

# Output super-enhancers
super_enhancers <- enhancers[enhancers$type == 'Super-Enhancer', ]
write.table(super_enhancers, 'super_enhancers.txt', sep = '\t', quote = FALSE, row.names = FALSE)
```

## Calculate Enhancer Signal

```bash
# Get H3K27ac signal at peak regions
bedtools multicov -bams h3k27ac.bam -bed enhancer_peaks.bed > enhancer_counts.txt

# Normalize by peak size
awk 'BEGIN{OFS="\t"} {
    size = $3 - $2
    rpm = ($NF / TOTAL_READS) * 1e6
    rpkm = rpm / (size / 1000)
    print $0, rpkm
}' enhancer_counts.txt > enhancer_signal.txt
```

## Downstream Analysis

### Gene Assignment

```bash
# Assign super-enhancers to nearest genes
bedtools closest -a super_enhancers.bed -b genes.bed -d > se_gene_assignment.txt
```

### Compare Conditions

**Goal:** Find super-enhancers gained or lost between two experimental conditions.

**Approach:** Convert super-enhancer tables to GRanges objects and use subsetByOverlaps with invert to identify condition-specific super-enhancers.

```r
# Load SE from two conditions
se1 <- read.table('condition1_SE.txt', header=TRUE)
se2 <- read.table('condition2_SE.txt', header=TRUE)

# Find differential super-enhancers
library(GenomicRanges)
gr1 <- makeGRangesFromDataFrame(se1)
gr2 <- makeGRangesFromDataFrame(se2)

# Gained in condition 2
gained <- subsetByOverlaps(gr2, gr1, invert=TRUE)

# Lost in condition 2
lost <- subsetByOverlaps(gr1, gr2, invert=TRUE)
```

### Enrichment of Disease Variants

```bash
# Check if GWAS SNPs enriched in super-enhancers
bedtools intersect -a gwas_snps.bed -b super_enhancers.bed -wa -wb > snps_in_SE.txt

# Calculate enrichment
total_snps=$(wc -l < gwas_snps.bed)
snps_in_se=$(wc -l < snps_in_SE.txt)
se_coverage=$(awk '{sum += $3-$2} END {print sum}' super_enhancers.bed)
genome_size=3000000000

expected=$(echo "$total_snps * $se_coverage / $genome_size" | bc -l)
enrichment=$(echo "$snps_in_se / $expected" | bc -l)
echo "Enrichment: $enrichment"
```

## Complete Workflow

```bash
#!/bin/bash
set -euo pipefail

H3K27AC_BAM=$1
PEAKS_BED=$2
OUTPUT_DIR=$3

mkdir -p $OUTPUT_DIR

echo "=== Convert peaks to GFF ==="
awk 'BEGIN{OFS="\t"} {print $1,"peaks","enhancer",$2,$3,".",$6,".","ID="NR}' \
    $PEAKS_BED > $OUTPUT_DIR/peaks.gff

echo "=== Run ROSE ==="
python ROSE_main.py \
    -g HG38 \
    -i $OUTPUT_DIR/peaks.gff \
    -r $H3K27AC_BAM \
    -o $OUTPUT_DIR \
    -s 12500 \
    -t 2500

echo "=== Summary ==="
n_typical=$(grep -c "Typical" $OUTPUT_DIR/*_AllEnhancers.table.txt || echo 0)
n_super=$(wc -l < $OUTPUT_DIR/*_SuperEnhancers.table.txt)

echo "Typical enhancers: $n_typical"
echo "Super-enhancers: $n_super"
```

## Related Skills

- chip-seq/peak-calling - Call H3K27ac peaks first
- chip-seq/peak-annotation - Annotate SE to genes
- chip-seq/differential-binding - Compare SE between conditions
- data-visualization/genome-tracks - Visualize SE regions