peak-annotation

Annotate ChIP-seq peaks to genomic features and genes using ChIPseeker. Assign peaks to promoters, exons, introns, and intergenic regions. Find nearest genes and calculate distance to TSS. Generate an

SKILL.md

GitHub

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.37M human and 967k mouse elements (Moore JE et al 2026 Nature). 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/CA-CTCF/CA-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:

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):

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)

# 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.37M human + 967k mouse elements. Registry V4 uses an 8-class scheme (the older V3 "CTCF-only" and "DNase-H3K4me3" were renamed CA-CTCF and CA-H3K4me3):

Class Definition Marker pattern
PLS (Promoter-Like Signature) ≤ 200 bp of annotated TSS; high DNase + high H3K4me3 DNase + H3K4me3
pELS (Proximal Enhancer-Like Signature) ≤ 2 kb of TSS; enhancer-like (DNase + H3K27ac, low H3K4me3) DNase + H3K27ac
dELS (Distal Enhancer-Like Signature) > 2 kb of TSS; enhancer-like DNase + H3K27ac
CA-H3K4me3 Chromatin-accessible + H3K4me3, not TSS-proximal DNase + H3K4me3
CA-CTCF Chromatin-accessible + CTCF (potential boundary) DNase + CTCF
CA-TF Chromatin-accessible + TF binding DNase + TF
CA Chromatin-accessible only DNase
TF TF-bound, not highly accessible TF
# Download ENCODE cCRE BED from SCREEN (GRCh38, expanded Registry-V4, uncompressed)
wget https://downloads.wenglab.org/Registry-V4/GRCh38-cCREs.bed

# Intersect peaks with cCRE; -wa preserves peak coords, -wb adds cCRE class
bedtools intersect -a peaks.narrowPeak -b GRCh38-cCREs.bed -wa -wb \
    > peaks_ccre.tsv

Cross-referencing peaks against cCREs:

  • Indicates whether peaks overlap canonical regulatory elements
  • Provides the cCRE class (PLS / pELS / dELS / CA-CTCF / CA-H3K4me3)
  • Cell-type-specific activity profiles available via SCREEN web app

GREAT / rGREAT (Regulatory Domain Gene-Set Enrichment)

GREAT (McLean 2010) addresses two problems with standard gene-set enrichment on peaks:

  1. Peak-to-gene assignment via regulatory domains (not nearest TSS)
  2. Statistical correction for region-locus length bias

Regulatory domain rules (default):

  • Basal domain: -5 kb / +1 kb of TSS
  • Extension: up to 1 Mb in each direction, OR until reaching neighbor's basal domain
  • Each peak is assigned to ALL genes whose regulatory domain it overlaps (not just nearest)
library(rGREAT)

# Submit peaks for regulatory-domain gene-set enrichment
res <- great(gr = peaks, gene_sets = 'GO:BP', tss_source = 'TxDb.Hsapiens.UCSC.hg38.knownGene',
              biomart_dataset = 'hsapiens_gene_ensembl')

# Top enriched gene sets
table_results <- getEnrichmentTable(res)
head(table_results)

# Visualization (local great() returns a GreatObject -> plotRegionGeneAssociations)
plotVolcano(res)
plotRegionGeneAssociations(res)

GREAT is most appropriate for distal regulatory elements (enhancer ChIP, ATAC). For promoter-focused marks (H3K4me3), ChIP-Enrich is more standard.

ChIP-Enrich (Locus-Length-Adjusted Gene-Set Enrichment)

Welch 2014: standard gene-set enrichment on peak-associated genes systematically over-counts long genes. ChIP-Enrich models locus length as a covariate.

library(chipenrich)

res <- chipenrich(peaks = 'peaks.bed', genome = 'hg38',
                   genesets = 'GOBP', locusdef = 'nearest_tss',
                   out_name = 'chipenrich_out', n_cores = 4)
# Locus definitions: nearest_tss, nearest_gene, exon, intron, 1kb, 5kb, 10kb
# method= accepts chipenrich (default) or fet; broadenrich() and polyenrich() are separate functions

For broad marks (H3K27me3, H3K9me3): use the separate broadenrich(peaks = 'peaks.bed', genome = 'hg38', genesets = 'GOBP', locusdef = 'nearest_tss') function, which accounts for region width.

ENCODE-rE2G (Modern Enhancer-Gene Linking)

ENCODE-rE2G (2024) replaces ABC for cell types with ENCODE data. Cell-type-specific logistic-regression weights map distal enhancer peaks to target genes with higher accuracy than nearest-TSS or basal+extension.

See atac-seq/enhancer-gene-linking for full workflow; the same model applies to ChIP-seq enhancer marks (H3K27ac, H3K4me1, H3K4me2).

Per-Tool Failure Modes

ChIPseeker -- TxDb / annoDb genome mismatch

Trigger: Using hg19 TxDb on hg38-aligned BAMs / peaks.

Mechanism: Silent; ChIPseeker doesn't verify genome assembly.

Symptom: Annotated gene symbols look reasonable but distance-to-TSS is wrong; promoter / intron classifications drift.

Fix: Match TxDb to BAM alignment genome explicitly; verify with seqlevels(peaks) == seqlevels(txdb).

ChIPseeker -- Default overlap='TSS' decouples gene from feature

Trigger: Default annotation call on peaks in gene bodies.

Mechanism: overlap='TSS' assigns nearest gene by TSS; feature classification is independent of that gene.

Symptom: Annotation reports nearest_gene=X, feature=intron where the intron belongs to a different gene.

Fix: Pass overlap='all' for host-gene-consistent annotation; or accept TSS-only convention and clarify in methods.

ChIPseeker -- Custom TxDb has no annoDb

Trigger: Building TxDb from GTF and passing annoDb='org.Hs.eg.db'.

Mechanism: Custom TxDb lacks the gene_id-to-symbol mapping that org.Hs.eg.db provides; ChIPseeker silently returns NA for symbols.

Fix: Map symbols separately from the original GTF after annotation; strip Ensembl version suffixes before joining.

HOMER -- Hard-coded promoter window

Trigger: Needing a 2 kb or 5 kb promoter window with HOMER.

Mechanism: HOMER's promoter classification is hard-coded to -1 kb / +100 bp; not configurable.

Fix: Post-hoc reclassify using Distance to TSS column:

awk -F'\t' 'NR>1 { dist = ($10 < 0) ? -$10 : $10; \
    feat = (dist <= 2000) ? "promoter_custom" : $8; \
    print $2, $3, $4, $16, $10, feat }' OFS='\t' annotated.txt

GREAT -- Default regulatory domain inappropriate for some species / cell types

Trigger: Using default basal+extension on insect or compact-genome data.

Mechanism: 1 Mb maximum extension assumes vertebrate-scale enhancer-target distances; not appropriate for organisms with shorter regulatory ranges.

Fix: Adjust extension parameter; for non-default species, configure regulatory domain explicitly.

GREAT -- Hyper-ChIPable peaks inflate enrichment

Trigger: Including unfiltered peaks at rRNA / housekeeping / mtDNA in GREAT analysis.

Mechanism: Hyper-ChIPable artifacts are enriched at highly-transcribed loci; GREAT assigns them to associated genes, inflating GO terms for "translation" and "ribosomal" categories.

Symptom: Top enriched GO terms always include "ribosomal", "translation", "mitochondrion" regardless of biology.

Fix: Blacklist filter + custom hyper-ChIPable filter (top-1% input signal) before GREAT.

ENCODE cCRE -- Cell-type-agnostic vs specific

Trigger: Using the master cCRE BED (cell-type-agnostic) to claim cell-type-specific regulatory activity.

Mechanism: Master cCRE BED is the union across all cell types. Specific activity profile per cell type is a separate dataset.

Fix: Use SCREEN web app or per-cell-type activity profiles for cell-type-specific claims.

Reconciliation: When Methods Disagree

Pattern Likely cause Action
ChIPseeker nearest-TSS gene ≠ HOMER nearest gene Different TSS reference; HOMER uses RefSeq Verify both use same TxDb / RefSeq + UCSC knownGene
GREAT enrichment ≠ ChIP-Enrich enrichment GREAT uses regulatory domain; ChIP-Enrich uses locus length adjustment Both are valid; use GREAT for distal regulatory, ChIP-Enrich for promoter-focused
Peak overlaps cCRE but classified differently than expected Cell-type-specific activity profile not used Check SCREEN per-cell-type profile
Enhancer peak's nearest gene differs from ENCODE-rE2G target ENCODE-rE2G uses cell-type chromatin context Use ENCODE-rE2G for cell-type-specific enhancer-gene claims

Common Errors

Error / symptom Cause Solution
seqlevels mismatch in ChIPseeker chr vs no-chr naming seqlevelsStyle(peaks) <- 'UCSC'
Gene symbols all NA in ChIPseeker Custom TxDb without annoDb Map symbols from original GTF
HOMER reports "no annotation" Genome not installed perl configureHomer.pl -install hg38
rGREAT timeout Large peak set + slow biomart Use pre-computed gene sets; lower peak count
chipenrich slow Default locusdef computed on-the-fly Use built-in locusdef shortcuts (nearest_tss, 1kb)
pyranges feature-overlap result missing strand pyranges 0.x conversion drops strand by default Pass strandedness='same' to overlap operations

References

  • Yu G et al 2015 Bioinformatics 31:2382 (ChIPseeker)
  • Heinz S et al 2010 Mol Cell 38:576 (HOMER annotatePeaks)
  • McLean CY et al 2010 Nat Biotechnol 28:495 (GREAT)
  • Gu Z 2023 Bioinformatics 39:btac745 (rGREAT)
  • Welch RP et al 2014 Nucleic Acids Res 42:e105 (ChIP-Enrich)
  • Cavalcante RG, Lee C, Welch RP, ... Sartor MA 2014 Bioinformatics 30:i393-i400 (Broad-Enrich)
  • ENCODE Project Consortium 2020 Nature 583:699 (cCRE registry v1)
  • Moore JE et al 2026 Nature (expanded cCRE registry, 2.37M human + 967k mouse elements)
  • Fulco CP et al 2019 Nat Genet 51:1664 (ABC model precursor to ENCODE-rE2G)
  • Kundaje lab / ENCODE 2024 (ENCODE-rE2G)
  • SCREEN: screen.encodeproject.org

Related Skills

  • chip-seq/peak-calling - Generate peaks for annotation
  • chip-seq/chipseq-qc - Filter hyper-ChIPable peaks before GREAT / chipenrich
  • chip-seq/super-enhancers - Annotate SE-associated genes
  • chip-seq/differential-binding - Annotate differential peaks
  • atac-seq/enhancer-gene-linking - ENCODE-rE2G workflow for cell-type-specific E-G linking
  • pathway-analysis/go-enrichment - Standard GO enrichment on peak-associated genes
  • pathway-analysis/reactome-pathways - Reactome pathway enrichment
  • genome-intervals/gtf-gff-handling - Parse and convert GTF/GFF
  • genome-intervals/proximity-operations - bedtools closest and window operations