Functional Enrichment Analysis

Over-Representation Analysis (ORA) vs. Gene Set Enrichment Analysis (GSEA)

Over-Representation Analysis (ORA)

ORA is a statistical method used to determine whether a predefined set of genes (e.g., genes associated with a specific biological pathway or function) is over-represented in a list of genes of interest (e.g., differentially expressed genes) compared to a reference list (e.g., all expressed genes).

Gene Set Enrichment Analysis (GSEA)

GSEA is a method used to determine whether predefined sets of genes show statistically significant, concordant differences between two biological states (e.g., phenotypes) across an entire ranked list of genes, rather than focusing only on a subset of significant genes.

Differences Between ORA and GSEA

Gene sets and pathway

A gene set is an unordered collection of genes that are functionally related. A pathway can be interpreted as a gene set by ignoring functional relationships among genes.

Gene Ontology (GO)

Gene Ontology defines concepts/classes used to describe gene function, and relationships between these concepts. It classifies functions along three aspects:

• MF: Molecular Function
  • molecular activities of gene products
• CC: Cellular Component
  • where gene products are active
• BP: Biological Process
  • pathways and larger processes made up of the activities of multiple gene products

GO terms are organized in a directed acyclic graph, where edges between terms represent parent-child relationship.

Kyoto Encyclopedia of Genes and Genomes (KEGG)

KEGG is a collection of manually drawn pathway maps representing molecular interaction and reaction networks. These pathways cover a wide range of biochemical processes that can be divided into 7 broad categories:

1. Metabolism
2. Genetic information processing
3. Environmental information processing
4. Cellular processes
5. Organismal systems
6. Human diseases
7. Drug development.

Other gene sets
GO and KEGG are the most frequently used for functional analysis. They are typically the first choice because of their long-standing curation and availability for a wide range of species. Other gene sets include but are not limited to Disease Ontology (DO), Disease Gene Network (DisGeNET), wikiPathways, Molecular Signatures Database MSigDb)

# Install and load necessary packages
install.packages("pheatmap")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("biomaRt")
BiocManager::install("DESeq2")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("ReactomePA")
BiocManager::install("fgsea")
BiocManager::install("msigdbr")
BiocManager::install("EnhancedVolcano")
install.packages('ggridges')
library(ReactomePA)
library(ggridges)
library(EnhancedVolcano)
library(pheatmap)
library(biomaRt)
library(DESeq2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(fgsea)
library(msigdbr)

# Load count data (replace the file path with your actual path)
count_data <- read.csv("/content/count_data.csv", row.names = 1)

# Create metadata for the full dataset
metadata <- data.frame(
  row.names = colnames(count_data),
  condition = c(rep("Normal", 3), rep("HER2", 5), rep("NonTNBC", 6), rep("TNBC", 5))
)

# Convert non-integer values to integers by rounding
count_data <- round(count_data)

# Create DESeq2 dataset for the full dataset
dds <- DESeqDataSetFromMatrix(countData = count_data, colData = metadata, design = ~ condition)

# Pre-filtering (Optional) Remove rows with low counts to improve computation
dds <- dds[rowSums(counts(dds)) > 1, ]

# Normalize the data and run DESeq2 analysis
dds <- DESeq(dds)

# Extract normalized counts
normalized_counts <- counts(dds, normalized = TRUE)

# Perform differential expression analysis for each condition vs Normal
res_her2_vs_normal <- results(dds, contrast = c("condition", "HER2", "Normal"))
res_non_tnbc_vs_normal <- results(dds, contrast = c("condition", "NonTNBC", "Normal"))
res_tnbc_vs_normal <- results(dds, contrast = c("condition", "TNBC", "Normal"))

# Define significance thresholds
padj_threshold <- 0.05
log2fc_threshold <- 1

# Extract significant results (adjusted p-value < 0.05, log2FoldChange > 1)
extract_significant <- function(res, padj_threshold, log2fc_threshold) {
  subset(res, padj < padj_threshold & abs(log2FoldChange) > log2fc_threshold)
}

sig_HER2_vs_Normal <- extract_significant(res_her2_vs_normal, padj_threshold, log2fc_threshold)
sig_NonTNBC_vs_Normal <- extract_significant(res_non_tnbc_vs_normal, padj_threshold, log2fc_threshold)
sig_TNBC_vs_Normal <- extract_significant(res_tnbc_vs_normal, padj_threshold, log2fc_threshold)

# Top genes for each comparison
top_genes <- function(sig_res) {
  rownames(head(sig_res[order(sig_res$padj), ], 50))
}

top_genes_HER2 <- top_genes(sig_HER2_vs_Normal)
top_genes_NonTNBC <- top_genes(sig_NonTNBC_vs_Normal)
top_genes_TNBC <- top_genes(sig_TNBC_vs_Normal)

# Convert Ensembl IDs to Gene Symbols and Entrez IDs
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
convert_ids <- function(ensembl_ids) {
    ids <- getBM(filters = "ensembl_gene_id",
                 attributes = c("ensembl_gene_id", "hgnc_symbol", "entrezgene_id"),
                 values = ensembl_ids,
                 mart = mart)
    ids <- ids[match(ensembl_ids, ids$ensembl_gene_id), ]  # Ensure the order matches the input
    ids$hgnc_symbol[is.na(ids$hgnc_symbol)] <- ids$ensembl_gene_id[is.na(ids$hgnc_symbol)]
    return(ids)
}

# Convert top genes to gene symbols and Entrez IDs
top_genes_HER2_ids <- convert_ids(top_genes_HER2)
top_genes_NonTNBC_ids <- convert_ids(top_genes_NonTNBC)
top_genes_TNBC_ids <- convert_ids(top_genes_TNBC)

# Add gene symbols to significant results
add_gene_symbols <- function(sig_res) {
  sig_res$gene_symbol <- convert_ids(rownames(sig_res))
  sig_res <- sig_res[, c("gene_symbol", names(sig_res)[-ncol(sig_res)])]
  sig_res[order(sig_res$padj), ]
}

sig_HER2_vs_Normal <- add_gene_symbols(sig_HER2_vs_Normal)
sig_NonTNBC_vs_Normal <- add_gene_symbols(sig_NonTNBC_vs_Normal)
sig_TNBC_vs_Normal <- add_gene_symbols(sig_TNBC_vs_Normal)

# Save significant results to CSV files
write.csv(sig_HER2_vs_Normal, "sig_her2_vs_normal.csv", row.names = FALSE)
write.csv(sig_NonTNBC_vs_Normal, "sig_non_tnbc_vs_normal.csv", row.names = FALSE)
write.csv(sig_TNBC_vs_Normal, "sig_tnbc_vs_normal.csv", row.names = FALSE)

# Extract Entrez IDs for enrichment analysis
extract_entrez_ids <- function(sig_res) {
  na.omit(convert_ids(rownames(sig_res))$entrezgene_id)
}

entrez_ids_her2_vs_normal <- extract_entrez_ids(sig_HER2_vs_Normal)
entrez_ids_non_tnbc_vs_normal <- extract_entrez_ids(sig_NonTNBC_vs_Normal)
entrez_ids_tnbc_vs_normal <- extract_entrez_ids(sig_TNBC_vs_Normal)

# Perform GO enrichment analysis
perform_go_enrichment <- function(entrez_ids, comparison_name) {
  ego <- enrichGO(gene          = entrez_ids,
                  OrgDb         = org.Hs.eg.db,
                  keyType       = "ENTREZID",
                  ont           = "ALL",
                  pAdjustMethod = "BH",
                  pvalueCutoff  = 0.05,
                  qvalueCutoff  = 0.05)
  write.csv(as.data.frame(ego), paste0("go_enrichment_", comparison_name, ".csv"))
  return(ego)
}

ego_her2_vs_normal <- perform_go_enrichment(entrez_ids_her2_vs_normal, "her2_vs_normal")
ego_non_tnbc_vs_normal <- perform_go_enrichment(entrez_ids_non_tnbc_vs_normal, "non_tnbc_vs_normal")
ego_tnbc_vs_normal <- perform_go_enrichment(entrez_ids_tnbc_vs_normal, "tnbc_vs_normal")

# Perform KEGG pathway enrichment analysis
perform_kegg_enrichment <- function(entrez_ids, comparison_name) {
  ekegg <- enrichKEGG(gene         = entrez_ids,
                      organism     = "hsa",
                      pvalueCutoff = 0.05)
  write.csv(as.data.frame(ekegg), paste0("kegg_enrichment_", comparison_name, ".csv"))
  return(ekegg)
}

ekegg_her2_vs_normal <- perform_kegg_enrichment(entrez_ids_her2_vs_normal, "her2_vs_normal")
ekegg_non_tnbc_vs_normal <- perform_kegg_enrichment(entrez_ids_non_tnbc_vs_normal, "non_tnbc_vs_normal")
ekegg_tnbc_vs_normal <- perform_kegg_enrichment(entrez_ids_tnbc_vs_normal, "tnbc_vs_normal")

# Perform Reactome Pathway enrichment analysis
reactome_her2_vs_normal <- enrichPathway(gene          = entrez_ids_her2_vs_normal,
                                         organism      = "human",
                                         pvalueCutoff  = 0.05,
                                         pAdjustMethod = "BH")

reactome_non_tnbc_vs_normal <- enrichPathway(gene          = entrez_ids_non_tnbc_vs_normal,
                                             organism      = "human",
                                             pvalueCutoff  = 0.05,
                                             pAdjustMethod = "BH")

reactome_tnbc_vs_normal <- enrichPathway(gene          = entrez_ids_tnbc_vs_normal,
                                         organism      = "human",
                                         pvalueCutoff  = 0.05,
                                         pAdjustMethod = "BH")

# Save Reactome enrichment results to CSV files
write.csv(as.data.frame(reactome_her2_vs_normal), "reactome_enrichment_her2_vs_normal.csv")
write.csv(as.data.frame(reactome_non_tnbc_vs_normal), "reactome_enrichment_non_tnbc_vs_normal.csv")
write.csv(as.data.frame(reactome_tnbc_vs_normal), "reactome_enrichment_tnbc_vs_normal.csv")

go_enrichment_her2_vs_normal.csv

go_enrichment_tnbc_vs_normal (4).csv

go_enrichment_non_tnbc_vs_normal (5).csv