#===============================================================================
# Gene Set Enrichment Analysis (GSEA) Workshop - R Code
# Lab 2c, Pathway and Network Analysis 2026
#
# This file contains only the R code for running fgsea.
#===============================================================================

#===============================================================================
# OPTIONAL EXERCISE: Running GSEA in R with fgsea (30 min)
#===============================================================================
#
# This optional section demonstrates how to run GSEA programmatically in R
# using the fgsea package. This approach is:
#   - Much faster than the desktop application (10-100x)
#   - Reproducible and scriptable
#   - Easy to integrate into existing R pipelines
#   - Produces equivalent results to GSEA desktop
#
# Prerequisites:
#   - Completed Exercise 1 (have the ranked gene list)
#   - R packages: fgsea, ggplot2, ggrepel 
#
#===============================================================================

#-------------------------------------------------------------------------------
# STEP O.1: Install and Load Required Packages
#-------------------------------------------------------------------------------

# Install Bioconductor packages (run once)
# install.packages("BiocManager")
# BiocManager::install("fgsea")

# Install CRAN packages (run once)
# install.packages(c("ggplot2", "ggrepel"))

# Load libraries
library(fgsea)
library(ggplot2)
library(ggrepel) #For placing labels neatly on plots 

#-------------------------------------------------------------------------------
# STEP O.2: Load the Pathway Database (GMT file) and ranked gene list
#-------------------------------------------------------------------------------

# fgsea can read GMT files directly using gmtPathways()
# Use the same GMT file as the desktop GSEA

pathways <- gmtPathways("Human_GOBP_AllPathways_noPFOCR_no_GO_iea_May_01_2026_symbol.gmt")

# Examine the structure
length(pathways)  # Number of gene sets
names(pathways)[1:5]  # First 5 pathway names
head(pathways[[1]])  # Genes in first pathway

# Check gene set sizes
pathway_sizes <- sapply(pathways, length)
summary(pathway_sizes)
hist(pathway_sizes, breaks = 50, main = "Gene Set Size Distribution",
     xlab = "Number of Genes", col = "steelblue")

# Load ranked gene list
ranked_gene_list <- read.table('ranked_basal_vs_classical.rnk', sep='\t')
# fgsea expects a named vector 
gene_ranks <- ranked_gene_list[,2,drop=TRUE]
names(gene_ranks) <- ranked_gene_list[,1]

#-------------------------------------------------------------------------------
# STEP O.4: Run fgsea
#-------------------------------------------------------------------------------

# Run fgsea with multilevel algorithm (recommended)
# This is the same algorithm as GSEA desktop but much faster

set.seed(42)  # For reproducibility
fgsea_results <- fgseaMultilevel(
  pathways = pathways,
  stats = gene_ranks,
  minSize = 15,       # Minimum gene set size
  maxSize = 500,      # Maximum gene set size
  eps = 0,            # Boundary for p-value estimation (0 = compute exactly)
  nPermSimple = 10000 # Number of permutations for simple estimation
)

# View the results
head(fgsea_results[order(pval), ], 10)

# Results columns:
# - pathway: Gene set name
# - pval: Nominal p-value
# - padj: Adjusted p-value (BH-FDR)
# - log2err: Standard error of the p-value estimate
# - ES: Enrichment Score
# - NES: Normalized Enrichment Score
# - size: Number of genes in the gene set (after filtering)
# - leadingEdge: List of leading edge genes


#-------------------------------------------------------------------------------
# STEP O.5: Explore the Results
#-------------------------------------------------------------------------------

# Count significant gene sets
cat("Gene sets tested:", nrow(fgsea_results), "\n")
cat("Significant at FDR < 0.25:", sum(fgsea_results$padj < 0.25, na.rm = TRUE), "\n")
cat("Significant at FDR < 0.05:", sum(fgsea_results$padj < 0.05, na.rm = TRUE), "\n")

# Top upregulated pathways (positive NES)
top_up <- fgsea_results[NES > 0][order(padj), ][1:10, ]
print(top_up[, .(pathway, NES, padj, size)])

# Top downregulated pathways (negative NES)
top_down <- fgsea_results[NES < 0][order(padj), ][1:10, ]
print(top_down[, .(pathway, NES, padj, size)])


#-------------------------------------------------------------------------------
# STEP O.6: Visualize Results - Enrichment Plot
#-------------------------------------------------------------------------------

# Plot enrichment for a specific pathway
# This is equivalent to the GSEA desktop enrichment plot

# Choose a significant pathway
top_pathway <- fgsea_results[order(padj), ]$pathway[1]
cat("Plotting:", top_pathway, "\n")

# Create enrichment plot
plotEnrichment(pathways[[top_pathway]], gene_ranks) +
  labs(title = top_pathway) +
  theme_minimal()

# Save the plot
ggsave("fgsea_enrichment_plot.png", width = 8, height = 5, dpi = 300)


#-------------------------------------------------------------------------------
# STEP O.7: Visualize Results - Summary Table Plot
#-------------------------------------------------------------------------------

# Create a summary plot of top pathways (both directions)

# Select top 10 from each direction
top_pathways <- rbind(
  fgsea_results[NES > 0][order(padj), ][1:10, ],
  fgsea_results[NES < 0][order(padj), ][1:10, ]
)

# Clean up pathway names for plotting (remove database prefixes)
top_pathways$pathway_short <- gsub("^.*%", "", top_pathways$pathway)
top_pathways$pathway_short <- gsub("_", " ", top_pathways$pathway_short)
top_pathways$pathway_short <- substr(top_pathways$pathway_short, 1, 50)  # Truncate

# Create bar plot
ggplot(top_pathways, aes(x = reorder(pathway_short, NES), y = NES, fill = NES > 0)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = c("TRUE" = "firebrick", "FALSE" = "steelblue"),
                    labels = c("TRUE" = "Upregulated", "FALSE" = "Downregulated"),
                    name = "Direction") +
  labs(
    title = "Top Enriched Pathways (fgsea)",
    subtitle = "Basal vs Classical",
    x = "",
    y = "Normalized Enrichment Score (NES)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

ggsave("fgsea_top_pathways.png", width = 10, height = 8, dpi = 300)


#-------------------------------------------------------------------------------
# STEP O.8: Extract Leading Edge Genes
#-------------------------------------------------------------------------------

# The leadingEdge column contains a list of genes driving each enrichment
# Let's extract them for the top pathway

top_pathway <- fgsea_results[order(padj), ]$pathway[1]
leading_edge <- fgsea_results[pathway == top_pathway, ]$leadingEdge[[1]]

cat("Leading edge genes for:", top_pathway, "\n")
cat("Number of genes:", length(leading_edge), "\n")
print(leading_edge)

# Create a data frame of leading edge genes with their scores
leading_edge_df <- data.frame(
  gene = leading_edge,
  rank_score = gene_ranks[leading_edge]
)
leading_edge_df <- leading_edge_df[order(-leading_edge_df$rank_score), ]
head(leading_edge_df, 20)


#-------------------------------------------------------------------------------
# STEP O.9: Export Results for EnrichmentMap
#-------------------------------------------------------------------------------

# fgsea results can be exported in a format compatible with EnrichmentMap
# We need to create files similar to GSEA desktop output

# Export all results as a table
write.table(fgsea_results[, c('pathway', 'pval', 'padj', 'ES', 'NES', 'size')],
       file = "fgsea_results.txt", sep = "\t")

# Export in GSEA-like format for EnrichmentMap
# Positive results
pos_results <- fgsea_results[NES > 0][order(pval), ]
pos_results$FWER <- p.adjust(pos_results$pval, method = "bonferroni")
pos_results$RANK_AT_MAX <- NA  # fgsea doesn't provide this directly

pos_export <- data.frame(
  NAME = pos_results$pathway,
  GS_SIZE = pos_results$size,
  ES = pos_results$ES,
  NES = pos_results$NES,
  NOM_p_val = pos_results$pval,
  FDR_q_val = pos_results$padj,
  FWER_p_val = pos_results$FWER
)
write.table(pos_export, "fgsea_report_pos.tsv", sep = "\t", 
            row.names = FALSE, quote = FALSE)

# Negative results
neg_results <- fgsea_results[NES < 0][order(pval), ]
neg_results$FWER <- p.adjust(neg_results$pval, method = "bonferroni")

neg_export <- data.frame(
  NAME = neg_results$pathway,
  GS_SIZE = neg_results$size,
  ES = neg_results$ES,
  NES = neg_results$NES,
  NOM_p_val = neg_results$pval,
  FDR_q_val = neg_results$padj,
  FWER_p_val = neg_results$FWER
)
write.table(neg_export, "fgsea_report_neg.tsv", sep = "\t",
            row.names = FALSE, quote = FALSE)

cat("\nFiles exported for EnrichmentMap:\n")
cat("  - fgsea_report_pos.tsv\n")
cat("  - fgsea_report_neg.tsv\n")

#===============================================================================
# OPTIONAL EXERCISE: Questions
#===============================================================================
#
# 1. How do the fgsea results compare to GSEA desktop results?
#    - Are the top pathways similar?
#    - Are the NES values comparable?
#    - Why do we expect fgsea results to differ from GSEA desktop results?
#    - Which approach is better for reproducibility? 
#
# 2. Run fgsea with different parameters:
#    - Try maxSize = 200 instead of 500. How do results change?
#    - Try minSize = 10. What happens?
#
# 3. Extract leading edge genes for the top 5 pathways.
#    - Are there genes that appear in multiple leading edges?
#    - What might this tell you biologically?
#
# 4. Time comparison:
#    - How long did fgsea take to run?
#    - How does this compare to GSEA desktop?
#
#===============================================================================
