Bio Chipseq Peak Annotation

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

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.

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

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

# 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.