17Functional enrichment analysis of genomic coordinations
Functional analysis using NGS data (e.g., RNA-Seq and ChIP-Seq) can be performed by linking coding and non-coding regions to coding genes via the ChIPseeker(Yu et al. 2015) package, which can annotate genomic regions to their nearest genes, host genes, and flanking genes respectively. In addition, it provides a function, seq2gene(), that simultaneously considers host genes, promoter regions and flanking genes from intergenic regions that may be under control via cis-regulation. This function maps genomic regions to genes in a many-to-many manner and facilitates functional analysis.
17.1 Chapter overview
Aspect
Genomic coordination enrichment
Questions
How can genomic regions such as peaks be mapped to candidate genes before functional enrichment?
Input
Annotated genomic regions, a genome/annotation configuration, and a region-to-gene mapping rule such as nearest, host, promoter, or flanking genes.
Methods
Region annotation and seq2gene() many-to-many coordination, followed by standard ORA, GSEA, comparison, or visualization workflows.
Main functions
annotatePeak(), seq2gene(), and downstream enrichGO() / enricher() / GSEA() after constructing gene-level inputs.
Output
Region annotations, candidate gene mappings, and downstream enrichResult / gseaResult objects.
Main limitations
Region-to-gene assignment is an inference, not proof of regulation. Results depend on genome build, annotation release, promoter/flanking parameters, and many-to-many mapping choices.
After mapping regions to genes, return to the shared data contract and choose an enrichment route. Preserve the genome build, annotation release, and mapping rule with the result.
17.2 Annotating peaks
annotatePeak() places every peak in its genomic context relative to transcript start sites. The example uses the peak table bundled with ChIPseeker and TxDb.Hsapiens.UCSC.hg19.knownGene, so it runs offline against a fixed annotation release rather than a downloaded one:
The printed GRanges object confirms the chromosome style (chr10) and the six metadata columns carried over from the narrowPeak-style table. That chromosome naming must agree with the naming style of the TxDb used downstream: a GRanges holding hg19 coordinates annotated against an hg38 TxDb produces the wrong features rather than an error, so record the genome build with the result.
anno <-annotatePeak(peaks,tssRegion =c(-3000, 3000),TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene)
annotatePeak() prints a progress banner while it works, and that banner carries the build timestamp, so it is suppressed here (output: false) and the resulting summary is shown on its own:
Annotated peaks generated by ChIPseeker
150/150 peaks were annotated
Genomic Annotation Summary:
Feature Frequency
8 Promoter (<=1kb) 11.3333333
9 Promoter (1-2kb) 4.0000000
10 Promoter (2-3kb) 6.0000000
4 5' UTR 0.6666667
3 3' UTR 3.3333333
1 1st Exon 6.0000000
6 Other Exon 7.3333333
2 1st Intron 20.0000000
7 Other Intron 25.3333333
5 Distal Intergenic 16.0000000
The printed summary is the distribution of annotation features relative to tssRegion. as.data.frame() exposes the per-peak call together with the gene responsible for it:
annotation is a feature label — Promoter (<=1kb), 1st Intron, 3' UTR, Distal Intergenic, and so on — and geneId is the Entrez identifier of the gene carrying that feature. distanceToTSS is signed: a negative value means the peak lies upstream of the transcript start site.
17.3 From regions to genes with seq2gene()
annotatePeak() assigns one label per peak. seq2gene() answers a different question — which genes may be regulated by each region — by collecting the genes whose promoter, gene body, or flanking window overlaps the region. The mapping is many-to-many: one region can link several genes, and several regions can converge on the same gene.
Check what that vector actually contains before feeding it anywhere. Its documented return value is a gene vector, and here only 12 of 152 entries are bare Entrez identifiers; the rest are composite transcriptId.geneId labels inherited from the TxDb’s gene naming. Splitting them with a pattern such as sub(".*[._]", "", x) looks reasonable and is wrong — on ENST00000446298.1_3 it returns "3", and those bogus values are simply not found during ID lookup, so enrichment quietly runs on a shrunk gene set.
The robust identifier source is annotatePeak()’s geneId column, which is numeric Entrez by construction. Use seq2gene() for region-to-gene connectivity — how many genes a region could reach, and which regions share a gene — and take the identifier list from the annotation table.
Read those three numbers together: 150 input peaks produce 152 region-to-gene links under the chosen window, and the annotated peaks resolve to 149 distinct candidate genes. The distinct count can be smaller than the peak count because several peaks may map to the same gene — here 1 peak(s) are extra assignments to an already-counted gene (150 peaks → 149 genes). A peak is dropped only when its window reaches no annotated gene at all and its geneId is NA; in this example 0 peaks have an NA geneId, including the 24 Distal Intergenic peaks, which annotatePeak() assigns to their nearest flanking gene. Report both the deduplication and any truly unassigned peaks before interpreting the result.
tssRegion and flankDistance are the decisions that change the gene set most: widening them recruits more genes and pulls the signal toward broadly annotated housekeeping terms. Record both, plus the genome build, with the result.
17.4 Enriching the mapped genes
Once regions have become a character vector of identifiers, nothing downstream is specific to ChIP-Seq. The vector enters the same engine and knowledge connector used anywhere else in the book:
This returns 0 terms, and that is the honest output rather than a broken pipeline: the peak table bundled with ChIPseeker is synthetic, not a transcription-factor experiment, so there is no biological theme in it to recover. Even with the thresholds loosened to pvalueCutoff = 1, the best adjusted p-value reachable from these peaks is roughly 0.42 — nowhere near interpretable. Everything above still ran end to end, the object is a normal enrichResult, and every downstream route from the data contract (simplify(), dotplot(), compareCluster()) applies once you point this chain at your own peak set.
An empty ORA result is a diagnostic, not a dead end. Before re-checking the identifiers, look at the two decisions this chapter owns: the tssRegion / flankDistance window, and how many peaks survived to contribute a geneId. A background that is too wide, or a window that silently drops most regions, both empty the result table without raising an error.
Region-to-gene assignment is an inference. A peak inside a gene body, or in the neighbouring intergenic window, does not establish that the gene is regulated by the assayed factor. Report how many regions produced how many candidate genes before interpreting any term.