bio-chipseq-peak-annotation

What this skill does

## Version Compatibility

Reference examples tested with: ChIPseeker 1.38+, GenomicFeatures 1.54+, rtracklayer 1.62+, HOMER 4.11+, rGREAT 2.4+, chipenrich 2.26+, pyranges 0.0.129+, pandas 2.2+.

ENCODE cCRE registry expanded to 2.35M human and 927k mouse elements (Nature 2025; Vu Ernst expansion). SCREEN web app at screen.encodeproject.org provides browser access; ENCODE provides bed files for batch annotation.

# Peak Annotation

**"What genes and regulatory elements do my peaks correspond to?"** -> Assign each peak to a genomic feature (promoter, exon, intron, intergenic), its target gene (via nearest-TSS or host-gene), and where applicable an ENCODE cCRE class (PLS/pELS/dELS/CTCF-only/DNase-H3K4me3).

- R (gene-feature): `ChIPseeker::annotatePeak(peaks, TxDb=txdb)`
- CLI (gene-feature): `annotatePeaks.pl peaks.bed hg38 -gtf annotation.gtf`
- Python (custom): pyranges + pandas
- R (cCRE classification): intersect peaks with ENCODE cCRE BED from SCREEN
- R (gene-set enrichment): `rGREAT::great()` or `chipenrich::chipenrich()`

The single biggest source of misinterpretation is the **nearest-TSS vs host-gene** distinction (see below). For enhancer-driven biology, ENCODE-rE2G or ABC (in atac-seq/enhancer-gene-linking) is more accurate than nearest-TSS.

## Choosing an Annotation Approach

| Context | Recommended | Why |
|---------|-------------|-----|
| Standard genome, pre-built annotations available | ChIPseeker with TxDb package | Simplest; automatic gene symbol mapping via annoDb |
| Custom or project-specific GTF | ChIPseeker + makeTxDbFromGFF, HOMER -gtf, or pyranges | All three handle custom annotations |
| HOMER already in pipeline | HOMER annotatePeaks.pl | Reuses tag directory; combined with motif workflow |
| Fine-grained control | pyranges (Python) | Full control over priority rules, distance calculation |
| Enhancer peaks (distal regulatory) | GREAT / rGREAT | Regulatory domain assignment (basal + extension), not just nearest |
| Cell-type-specific enhancer-gene linking | ENCODE-rE2G | Modern (2024); ABC-trained logistic regression with chromatin context |
| Gene-set enrichment with locus-length adjustment | chipenrich / Broad-Enrich | Corrects for systematic gene-length bias in peak assignment |
| Compare against ENCODE cCRE atlas | SCREEN cCRE BED intersect | Cross-reference standard regulatory registry |
| Promoter-coverage decomposition | bedtools intersect with TSS windows | Quick stats per peak set |

**Critical:** Use the same annotation source as the alignment (UCSC knownGene TxDb with GENCODE GTF alignment causes mismatches). When a specific GTF is provided, use it directly via `makeTxDbFromGFF` rather than a mismatched pre-built TxDb package.

## Nearest-TSS vs Host-Gene Convention

Peak annotation involves two decisions that should be coupled but often aren't:
1. Which gene to assign (target gene)
2. What feature the peak overlaps (promoter / exon / intron / intergenic)

Default tools decouple these, producing internally inconsistent annotations.

| Convention | Gene from | Feature from | Tools |
|------------|-----------|---------------|-------|
| Nearest-TSS (default) | Gene with closest TSS | Physical overlap at peak center | ChIPseeker `overlap='TSS'` (default), HOMER |
| Host-gene priority | Gene whose body contains the peak | Same gene's features | ChIPseeker `overlap='all'` |

**Example failure:** Peak inside gene A's intron, near gene B's TSS. Default tools report `nearest_gene=B, feature=intron` — but the intron belongs to gene A, not gene B. The annotation is internally inconsistent.

### Choosing per Biology

| Context | Convention | Rationale |
|---------|-----------|-----------|
| Distal TF binding (enhancers) | Nearest-TSS, but prefer ENCODE-rE2G / ABC | Enhancers can regulate gene A despite sitting in gene B's intron |
| Histone marks in gene bodies (H3K36me3, H3K27me3) | Host-gene | Mark reflects host transcriptional state |
| Promoter-associated marks (H3K4me3, H3K27ac at promoters) | Either | Most peaks at promoters where conventions agree |
| Custom annotation against project GTF | Host-gene | Internal consistency |
| Reproducing published HOMER results | Nearest-TSS | Matches HOMER default |

When a task says "nearest gene," clarify which definition. For most annotation purposes where gene + feature should be consistent, use host-gene; for distal enhancer biology, use a proper enhancer-gene linker (ENCODE-rE2G, ABC).

## Coordinate Systems and TSS

BED uses 0-based half-open `[start, end)`. GTF uses 1-based closed `[start, end]`. Mixing without conversion shifts annotations by one base.

**Peak center (BED):** `(start + end) // 2`

**TSS from GTF (1-based to 0-based):**
- Plus-strand: `tss_0based = start - 1`
- Minus-strand: `tss_0based = end`

**Signed distance** (negative = upstream of TSS):
- Plus-strand: `distance = peak_center - tss`
- Minus-strand: `distance = -(peak_center - tss)`

## ChIPseeker (R)

**Goal:** Assign each ChIP-seq peak to a gene and a feature category using a transcript database.

**Approach:** Load the TxDb (pre-built or custom-built from GTF), pass peaks to `annotatePeak()` with the desired `tssRegion` window and `overlap` convention (host-gene vs nearest-TSS), then export the annotated data frame with gene symbols mapped from `annoDb` or the original GTF.

**Standard genome:**

```r
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

peaks <- readPeakFile('peaks.narrowPeak')
peak_anno <- annotatePeak(peaks,
                           TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene,
                           tssRegion = c(-2000, 2000),
                           annoDb = 'org.Hs.eg.db',
                           overlap = 'all')   # host-gene convention
anno_df <- as.data.frame(peak_anno)
```

**Custom GTF** (use makeTxDbFromGFF; map symbols from original GTF since custom TxDb objects lack annoDb mappings):

```r
library(GenomicFeatures)
library(rtracklayer)

txdb <- makeTxDbFromGFF('genes.gtf.gz', format = 'gtf')
peaks <- readPeakFile('peaks.bed')
peak_anno <- annotatePeak(peaks, TxDb = txdb, tssRegion = c(-2000, 2000),
                           overlap = 'all')

gtf <- import('genes.gtf.gz')
gene_map <- unique(data.frame(
    gene_id = sub('\\..*', '', gtf$gene_id),
    symbol = gtf$gene_name, stringsAsFactors = FALSE))
gene_map <- gene_map[!is.na(gene_map$symbol), ]
anno_df <- as.data.frame(peak_anno)
anno_df$gene_id_base <- sub('\\..*', '', anno_df$geneId)
anno_df$SYMBOL <- gene_map$symbol[match(anno_df$gene_id_base, gene_map$gene_id)]
```

GENCODE gene IDs have version suffixes (`ENSG00000142192.25`); strip before joining.

**Promoter window:** `tssRegion = c(-2000, 2000)` is common; `c(-3000, 3000)` is ChIPseeker default. Match to analysis requirements.

**Feature priority:** Default `Promoter > 5'UTR > 3'UTR > Exon > Intron > Downstream > Intergenic`. A peak in both a promoter (gene A) and an intron (gene B) receives "Promoter (gene A)" by default.

## HOMER annotatePeaks.pl (CLI)

```bash
# Standard genome (HOMER's installed annotation)
annotatePeaks.pl peaks.bed hg38 > annotated.txt

# Custom GTF (overrides HOMER's default)
annotatePeaks.pl peaks.bed hg38 -gtf genes.gtf > annotated.txt

# Without installed genome, GTF only
annotatePeaks.pl peaks.bed none -gtf genes.gtf > annotated.txt

# Generate annotation statistics
annotatePeaks.pl peaks.bed hg38 -gtf genes.gtf -annStats stats.txt > annotated.txt
```

HOMER's 19-column output: columns 8 (Annotation), 10 (Distance to TSS), 16 (Gene Name) are the primary annotation columns.

**HOMER promoter window is fixed at -1kb / +100bp** — not configurable via flags. For custom windows, reclassify using the Distance to TSS column post-hoc.

## ENCODE cCRE Classification

The ENCODE Registry of candidate cis-Regulatory Elements (cCREs) provides 2.35M human + 927k mouse elements classified into 5 categories:

| Class | Definition | Marker pattern |
|-------|------------|-----------------|
| **PLS** (Promo