omics.co.in logo omics.co.in

Chapter 10: Pathway Annotation & Functional Over-Representation: DAVID, Reactome, KEGG & iLINCS

Tested Stacks: GATK 4.6.1 Samtools 1.21 Nextflow 24.10 Seurat v5.1 Python 3.12
Prerequisites: Basic Linux CLI & Sequencing Formats
📚 Chapter 10 of 18 0% Course Completed

1. Over-Representation Analysis (ORA) & Hypergeometric Testing

Identifying a list of 500 differentially expressed genes is only the beginning. Biomedical insight requires contextualizing these genes into biological pathways, molecular functions, and cellular compartments. Over-Representation Analysis evaluates whether a predefined biological pathway contains statistically more DEGs than would be expected by random chance.

The probability of observing \(k\) or more pathway genes in an experimental list of \(n\) DEGs from a genome of \(N\) total genes containing \(M\) pathway members is calculated via the Hypergeometric Distribution / Fisher's Exact Test:

$$P(X \ge k) = \sum_{i=k}^{\min(n, M)} \frac{\binom{M}{i} \binom{N - M}{n - i}}{\binom{N}{n}}$$

2. Curated Biomedical Knowledgebases

DatabaseScope & FocusHierarchical StructurePrimary Clinical Utility
Gene Ontology (GO)Biological Process (BP), Molecular Function (MF), Cellular Component (CC)Directed Acyclic Graph (DAG) with parent-child relationshipsGlobal functional profiling; cellular location
KEGG PathwaysMetabolic maps, signal transduction pathways, human disease cascadesInterconnected graphical metabolic and signaling wiring diagramsDrug targets, oncology signaling pathways (e.g. MAPK, PI3K-Akt)
ReactomeCurated reaction-level biochemical events, complexes, and transitionsHierarchical molecular reaction pathwaysHigh-specificity mechanistic molecular transitions
iLINCS / L1000Integrative transcriptomic signatures of >20,000 small-molecule compoundsPerturbagen-induced gene expression signaturesComputational drug repurposing (reversing disease signatures)

3. Automated Pathway Annotation Script with clusterProfiler & ReactomePA

library(clusterProfiler)
library(org.Hs.eg.db) # Human gene annotation database
library(ReactomePA)

# 1. Extract Upregulated Genes and Convert Ensembl to Entrez IDs
sig_up_genes <- res_df %>%
    filter(Call == "UP") %>%
    pull(Gene_Symbol)

# Convert Gene Identifiers to Entrez IDs (Required by KEGG and Reactome)
entrez_ids <- bitr(
    sig_up_genes,
    fromType = "ENSEMBL",
    toType = "ENTREZID",
    OrgDb = org.Hs.eg.db
)$ENTREZID

# 2. Gene Ontology Biological Process (GO-BP) Over-Representation Analysis
go_bp <- enrichGO(
    gene = entrez_ids,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.05,
    readable = TRUE
)

# 3. Reactome Pathway Enrichment
reactome_res <- enrichPathway(
    gene = entrez_ids,
    pvalueCutoff = 0.05,
    readable = TRUE
)

# 4. Visualization of Enriched Pathways
cat("Top 5 Enriched Reactome Pathways:\n")
head(as.data.frame(reactome_res)[, c("Description", "GeneRatio", "p.adjust")], 5)

# Dotplot of Enriched Biological Processes
dotplot(go_bp, showCategory = 15, title = "Top 15 Enriched GO Biological Processes")
⚠️

Common Pitfalls & Troubleshooting Checklist

  • Java Heap Exhaustion (GATK / Picard): If processes terminate with java.lang.OutOfMemoryError, set your environment variable explicitly: export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp".
  • Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with samtools index -@ 4 sample.bam.
  • Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority: conda config --set channel_priority strict.
  • Chromosome Naming Contamination: Always verify whether reference FASTA headers use chr1 or 1; divergent naming causes silent alignment drops.
Was this guide helpful to your research?
Help us maintain 100% accurate, peer-reviewed bioinformatics standards.
📝 Propose Edit on GitHub ↗