Pathway & Functional Analysis

Step 5: Biological interpretation of DE results

Overview

Pathway analysis helps interpret differential expression results in a biological context by identifying enriched biological processes, molecular functions, and pathways.

Gene Ontology Enrichment

Over-Representation Analysis (ORA)

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(tidyverse)

# Load DE results
de_results <- read_csv("differential_expression_results.csv")

# Get significant upregulated genes
sig_up <- de_results |>
  filter(padj < 0.05, log2FoldChange > 1) |>
  pull(gene_id)

# Convert to Entrez IDs
entrez_up <- bitr(sig_up, 
                  fromType = "ENSEMBL",
                  toType = "ENTREZID",
                  OrgDb = org.Hs.eg.db)$ENTREZID

# GO enrichment analysis
go_bp <- enrichGO(
  gene = entrez_up,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.2,
  readable = TRUE
)

Visualize GO Results

# Dot plot
dotplot(go_bp, showCategory = 20) +
  ggtitle("GO Biological Process Enrichment")

# Bar plot
barplot(go_bp, showCategory = 15) +
  ggtitle("Top Enriched GO Terms")

# Enrichment map
emapplot(pairwise_termsim(go_bp), showCategory = 30)

# Gene-concept network
cnetplot(go_bp, categorySize = "pvalue", foldChange = NULL)

Gene Set Enrichment Analysis (GSEA)

GSEA uses the full ranked gene list rather than a cutoff-based approach.

# Prepare ranked gene list
gene_list <- de_results |>
  filter(!is.na(padj)) |>
  arrange(desc(log2FoldChange))

# Convert to Entrez and get log2FC
gene_list_entrez <- bitr(gene_list$gene_id,
                         fromType = "ENSEMBL",
                         toType = "ENTREZID",
                         OrgDb = org.Hs.eg.db)

ranked_genes <- gene_list |>
  inner_join(gene_list_entrez, by = c("gene_id" = "ENSEMBL")) |>
  arrange(desc(log2FoldChange)) |>
  distinct(ENTREZID, .keep_all = TRUE)

# Create named vector
gene_ranks <- ranked_genes$log2FoldChange
names(gene_ranks) <- ranked_genes$ENTREZID
gene_ranks <- sort(gene_ranks, decreasing = TRUE)

# Run GSEA with GO
gsea_go <- gseGO(
  geneList = gene_ranks,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  minGSSize = 10,
  maxGSSize = 500,
  pvalueCutoff = 0.05,
  verbose = FALSE
)

GSEA Visualization

# Ridge plot
ridgeplot(gsea_go, showCategory = 15) +
  ggtitle("GSEA Ridge Plot")

# GSEA plot for specific pathway
gseaplot2(gsea_go, geneSetID = 1:3, pvalue_table = TRUE)

KEGG Pathway Analysis

# KEGG enrichment
kegg_enrich <- enrichKEGG(
  gene = entrez_up,
  organism = "hsa",
  pvalueCutoff = 0.05
)

# Visualize
dotplot(kegg_enrich, showCategory = 15) +
  ggtitle("KEGG Pathway Enrichment")

# KEGG GSEA
gsea_kegg <- gseKEGG(
  geneList = gene_ranks,
  organism = "hsa",
  minGSSize = 10,
  pvalueCutoff = 0.05
)

MSigDB Gene Sets

library(msigdbr)

# Get Hallmark gene sets
hallmark <- msigdbr(species = "Homo sapiens", category = "H")

hallmark_list <- hallmark |>
  split(x = .$entrez_gene, f = .$gs_name)

# GSEA with Hallmark
gsea_hallmark <- GSEA(
  geneList = gene_ranks,
  TERM2GENE = hallmark |> select(gs_name, entrez_gene),
  pvalueCutoff = 0.05
)

dotplot(gsea_hallmark, showCategory = 20) +
  ggtitle("Hallmark Gene Set Enrichment")

Comparing Multiple Conditions

# For multiple comparisons, use compareCluster
gene_clusters <- list(
  Upregulated = entrez_up,
  Downregulated = entrez_down
)

compare_go <- compareCluster(
  geneCluster = gene_clusters,
  fun = "enrichGO",
  OrgDb = org.Hs.eg.db,
  ont = "BP"
)

dotplot(compare_go, showCategory = 10) +
  ggtitle("GO Enrichment Comparison")

Export Results

# Export enrichment results
go_results <- as.data.frame(go_bp)
write_csv(go_results, "go_enrichment_results.csv")

gsea_results <- as.data.frame(gsea_go)
write_csv(gsea_results, "gsea_results.csv")

kegg_results <- as.data.frame(kegg_enrich)
write_csv(kegg_results, "kegg_enrichment_results.csv")

Summary

This completes the bulk RNA-seq analysis pipeline. You now have:

1. ✅ Quality-controlled and trimmed reads
2. ✅ Aligned reads to reference genome
3. ✅ Gene count matrix
4. ✅ Differential expression results
5. ✅ Pathway and functional enrichment analysis

TipNext Steps

Consider additional analyses: - Gene regulatory network analysis - Transcription factor enrichment - Integration with other omics data - Validation experiments (qPCR)