pathway_errorbar(
  # ... existing parameters ...
  color_theme = "nature",              # Professional color schemes
  legend_position = "top",             # Flexible legend positioning
  legend_title = "Sample Groups",      # Custom legend titles
  pvalue_format = "smart",             # Smart p-value formatting
  pvalue_stars = TRUE,                 # Significance indicators
  pvalue_colors = TRUE,                # Color-coded significance
  pathway_class_text_color = "auto",   # Auto theme-matched colors
  accessibility_mode = TRUE            # Enhanced accessibility
)

# install.packages("devtools")
devtools::install_github("cafferychen777/ggpicrust2")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

pkgs <- c("phyloseq", "ALDEx2", "SummarizedExperiment", "Biobase", "devtools",
          "ComplexHeatmap", "BiocGenerics", "BiocManager", "metagenomeSeq",
          "Maaslin2", "edgeR", "lefser", "limma", "KEGGREST", "DESeq2")

for (pkg in pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE))
    BiocManager::install(pkg)
}

# If you want to analyze the abundance of KEGG pathways instead of KO within the pathway, please set `ko_to_kegg` to TRUE.
# KEGG pathways typically have more descriptive explanations.

library(readr)
library(ggpicrust2)
library(tibble)
library(tidyverse)
library(ggprism)
library(patchwork)

# Load necessary data: abundance data and metadata
abundance_file <- "path/to/your/abundance_file.tsv"
metadata <- read_delim(
    "path/to/your/metadata.txt",
    delim = "\t",
    escape_double = FALSE,
    trim_ws = TRUE
)

# Run ggpicrust2 with input file path
results_file_input <- ggpicrust2(file = abundance_file,
                                 metadata = metadata,
                                 group = "your_group_column", # For example dataset, group = "Environment"
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 ko_to_kegg = TRUE,
                                 order = "pathway_class",
                                 p_values_bar = TRUE,
                                 x_lab = "pathway_name")

# Run ggpicrust2 with imported data.frame
abundance_data <- read_delim(abundance_file, delim = "\t", col_names = TRUE, trim_ws = TRUE)

# Run ggpicrust2 with input data
results_data_input <- ggpicrust2(data = abundance_data,
                                 metadata = metadata,
                                 group = "your_group_column", # For example dataset, group = "Environment"
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 ko_to_kegg = TRUE,
                                 order = "pathway_class",
                                 p_values_bar = TRUE,
                                 x_lab = "pathway_name")

# Access the plot and results dataframe for the first DA method
example_plot <- results_file_input[[1]]$plot
example_results <- results_file_input[[1]]$results

# Use the example data in ggpicrust2 package
data(ko_abundance)
data(metadata)
results_file_input <- ggpicrust2(data = ko_abundance,
                                 metadata = metadata,
                                 group = "Environment",
                                 pathway = "KO",
                                 daa_method = "LinDA",
                                 ko_to_kegg = TRUE,
                                 order = "pathway_class",
                                 p_values_bar = TRUE,
                                 x_lab = "pathway_name")

# Analyze the EC or MetaCyc pathway
data(metacyc_abundance)
results_file_input <- ggpicrust2(data = metacyc_abundance,
                                 metadata = metadata,
                                 group = "Environment",
                                 pathway = "MetaCyc",
                                 daa_method = "LinDA",
                                 ko_to_kegg = FALSE,
                                 order = "group",
                                 p_values_bar = TRUE,
                                 x_lab = "description")
results_file_input[[1]]$plot
results_file_input[[1]]$results

library(readr)
library(ggpicrust2)
library(tibble)
library(tidyverse)
library(ggprism)
library(patchwork)

# If you want to analyze KEGG pathway abundance instead of KO within the pathway, turn ko_to_kegg to TRUE.
# KEGG pathways typically have more explainable descriptions.

# Load metadata as a tibble
# data(metadata)
metadata <- read_delim("path/to/your/metadata.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)

# Load KEGG pathway abundance
# data(kegg_abundance)
kegg_abundance <- ko2kegg_abundance("path/to/your/pred_metagenome_unstrat.tsv")

# Perform pathway differential abundance analysis (DAA) using ALDEx2 method
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = kegg_abundance, metadata = metadata, group = "Environment", daa_method = "ALDEx2", select = NULL, reference = NULL)

# Filter results for ALDEx2_Welch's t test method
# Please check the unique(daa_results_df$method) and choose one
daa_sub_method_results_df <- daa_results_df[daa_results_df$method == "ALDEx2_Wilcoxon rank test", ]

# Annotate pathway results using KO to KEGG conversion
daa_annotated_sub_method_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_sub_method_results_df, ko_to_kegg = TRUE)

# Generate pathway error bar plot
# Please change Group to metadata$your_group_column if you are not using example dataset
p <- pathway_errorbar(abundance = kegg_abundance, daa_results_df = daa_annotated_sub_method_results_df, Group = metadata$Environment, p_values_threshold = 0.05, order = "pathway_class", select = NULL, ko_to_kegg = TRUE, p_value_bar = TRUE, colors = NULL, x_lab = "pathway_name")

# If you want to analyze EC, MetaCyc, and KO without conversions, turn ko_to_kegg to FALSE.

# Load metadata as a tibble
# data(metadata)
metadata <- read_delim("path/to/your/metadata.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)

# Load KO abundance as a data.frame
# data(ko_abundance)
ko_abundance <- read.delim("path/to/your/pred_metagenome_unstrat.tsv")

# Perform pathway DAA using ALDEx2 method
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = ko_abundance %>% column_to_rownames("#NAME"), metadata = metadata, group = "Environment", daa_method = "ALDEx2", select = NULL, reference = NULL)

# Filter results for ALDEx2_Kruskal-Wallace test method
daa_sub_method_results_df <- daa_results_df[daa_results_df$method == "ALDEx2_Wilcoxon rank test", ]

# Annotate pathway results without KO to KEGG conversion
daa_annotated_sub_method_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_sub_method_results_df, ko_to_kegg = FALSE)

# Generate pathway error bar plot
# Please change column_to_rownames() to the feature column
# Please change Group to metadata$your_group_column if you are not using example dataset
p <- pathway_errorbar(abundance = ko_abundance %>% column_to_rownames("#NAME"), daa_results_df = daa_annotated_sub_method_results_df, Group = metadata$Environment, p_values_threshold = 0.05, order = "group",
select = daa_annotated_sub_method_results_df %>% arrange(p_adjust) %>% slice(1:20) %>% dplyr::select(feature) %>% pull(),
ko_to_kegg = FALSE,
p_value_bar = TRUE,
colors = NULL,
x_lab = "description")

# Workflow for MetaCyc Pathway and EC

# Load MetaCyc pathway abundance and metadata
data("metacyc_abundance")
data("metadata")

# Perform pathway DAA using LinDA method
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = "LinDA")

# Annotate MetaCyc pathway results without KO to KEGG conversion
metacyc_daa_annotated_results_df <- pathway_annotation(pathway = "MetaCyc", daa_results_df = metacyc_daa_results_df, ko_to_kegg = FALSE)

# Generate pathway error bar plot
# Please change column_to_rownames() to the feature column
# Please change Group to metadata$your_group_column if you are not using example dataset
pathway_errorbar(abundance = metacyc_abundance %>% column_to_rownames("pathway"), daa_results_df = metacyc_daa_annotated_results_df, Group = metadata$Environment, ko_to_kegg = FALSE, p_values_threshold = 0.05, order = "group", select = NULL, p_value_bar = TRUE, colors = NULL, x_lab = "description")

# Generate pathway heatmap with clustering
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
feature_with_p_0.05 <- metacyc_daa_results_df %>% filter(p_adjust < 0.05)
pathway_heatmap(
  abundance = metacyc_abundance %>% 
    filter(pathway %in% feature_with_p_0.05$feature) %>% 
    column_to_rownames("pathway"), 
  metadata = metadata, 
  group = "Environment",
  cluster_rows = TRUE,
  clustering_method = "ward.D2"
)

# Generate pathway PCA plot
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
pathway_pca(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment")

# Run pathway DAA for multiple methods
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
methods <- c("ALDEx2", "DESeq2", "edgeR")
daa_results_list <- lapply(methods, function(method) {
  pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = method)
})

# Compare results across different methods
comparison_results <- compare_daa_results(daa_results_list = daa_results_list, method_names = c("ALDEx2_Welch's t test", "ALDEx2_Wilcoxon rank test", "DESeq2", "edgeR"))

# Sample usage of the ko2kegg_abundance function
devtools::install_github('cafferychen777/ggpicrust2')

library(ggpicrust2)

# Assume that the KO abundance table is stored in a file named "ko_abundance.tsv"
ko_abundance_file <- "ko_abundance.tsv"

# Convert KO abundance to KEGG pathway abundance
kegg_abundance <- ko2kegg_abundance(file = ko_abundance_file)

# Alternatively, if the KO abundance data is already loaded as a data frame named "ko_abundance"
data("ko_abundance")
kegg_abundance <- ko2kegg_abundance(data = ko_abundance)

# The resulting kegg_abundance data frame can now be used for further analysis and visualization.

# The abundance table is recommended to be a data.frame rather than a tibble.
# The abundance table should have feature names or pathway names as row names, and sample names as column names.
# You can use the output of ko2kegg_abundance
ko_abundance_file <- "path/to/your/pred_metagenome_unstrat.tsv"
kegg_abundance <- ko2kegg_abundance(ko_abundance_file) # Or use data(kegg_abundance)

metadata <- read_delim("path/to/your/metadata.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)

# The default DAA method is "ALDEx2"
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = kegg_abundance, metadata = metadata, group = "Environment", daa_method = "linDA", select = NULL, p.adjust = "BH", reference = NULL)

# If you have more than 3 group levels and want to use the LinDA, limma voom, or Maaslin2 methods, you should provide a reference.
metadata <- read_delim("path/to/your/metadata.txt", delim = "\t", escape_double = FALSE, trim_ws = TRUE)

# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = kegg_abundance, metadata = metadata, group = "Group", daa_method = "LinDA", select = NULL, p.adjust = "BH", reference = "Harvard BRI")

# Other example
data("metacyc_abundance")
data("metadata")
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = "LinDA", select = NULL, p.adjust = "BH", reference = NULL)

library(ggpicrust2)
library(tidyverse)
data("metacyc_abundance")
data("metadata")

# Run pathway_daa function for multiple methods
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
methods <- c("ALDEx2", "DESeq2", "edgeR")
daa_results_list <- lapply(methods, function(method) {
  pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = method)
})

method_names <- c("ALDEx2","DESeq2", "edgeR")
# Compare results across different methods
comparison_results <- compare_daa_results(daa_results_list = daa_results_list, method_names = method_names)

# Make sure to check if the features in `daa_results_df` correspond to the selected pathway

# Annotate KEGG Pathway
data("kegg_abundance")
data("metadata")
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = kegg_abundance, metadata = metadata, group = "Environment", daa_method = "LinDA")

# Generic KO to KEGG pathway annotation (not specific to any organism)
daa_annotated_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df, ko_to_kegg = TRUE)

# Species-specific KEGG pathway annotation (e.g., for human)
human_annotated_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df, ko_to_kegg = TRUE, organism = "hsa")

# Species-specific KEGG pathway annotation (e.g., for E. coli)
ecoli_annotated_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df, ko_to_kegg = TRUE, organism = "eco")

# Annotate KO
data("ko_abundance")
data("metadata")
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(abundance = ko_abundance %>% column_to_rownames("#NAME"), metadata = metadata, group = "Environment", daa_method = "LinDA")
daa_annotated_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df, ko_to_kegg = FALSE)

# Annotate KEGG
# daa_annotated_results_df <- pathway_annotation(pathway = "EC", daa_results_df = daa_results_df, ko_to_kegg = FALSE)

# Annotate MetaCyc Pathway
data("metacyc_abundance")
data("metadata")
# Please change column_to_rownames() to the feature column if you are not using example dataset
# Please change group to "your_group_column" if you are not using example dataset
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = "LinDA")
metacyc_daa_annotated_results_df <- pathway_annotation(pathway = "MetaCyc", daa_results_df = metacyc_daa_results_df, ko_to_kegg = FALSE)

data("ko_abundance")
data("metadata")
kegg_abundance <- ko2kegg_abundance(data = ko_abundance) # Or use data(kegg_abundance)
# Please change group to "your_group_column" if you are not using example dataset
daa_results_df <- pathway_daa(kegg_abundance, metadata = metadata, group = "Environment", daa_method = "LinDA")
daa_annotated_results_df <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df, ko_to_kegg = TRUE)
# Please change Group to metadata$your_group_column if you are not using example dataset
p <- pathway_errorbar(abundance = kegg_abundance,
           daa_results_df = daa_annotated_results_df,
           Group = metadata$Environment,
           ko_to_kegg = TRUE,
           p_values_threshold = 0.05,
           order = "pathway_class",
           select = NULL,
           p_value_bar = TRUE,
           colors = NULL,
           x_lab = "pathway_name")

# Nature journal style
p_nature <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_annotated_results_df,
  Group = metadata$Environment,
  ko_to_kegg = TRUE,
  color_theme = "nature",              # Professional color theme
  legend_position = "top",             # Legend position control
  legend_title = "Sample Groups",      # Custom legend title
  pvalue_format = "smart",             # Smart p-value formatting
  pvalue_stars = TRUE                  # Significance stars (*, **, ***)
)

# Science journal style with enhanced p-value display
p_science <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_annotated_results_df,
  Group = metadata$Environment,
  ko_to_kegg = TRUE,
  color_theme = "science",             # Science journal theme
  legend_direction = "horizontal",     # Horizontal legend layout
  legend_title_size = 14,              # Larger legend title
  pvalue_format = "combined",          # Show both p-values and stars
  pvalue_colors = TRUE,                # Color-coded significance levels
  pathway_class_text_color = "auto"    # Auto pathway class colors
)

# Comprehensive legend customization
p_custom <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_annotated_results_df,
  Group = metadata$Environment,
  ko_to_kegg = TRUE,
  color_theme = "cell",
  # Legend positioning and styling
  legend_position = "bottom",          # bottom, top, left, right, none
  legend_direction = "horizontal",     # horizontal, vertical
  legend_title = "Treatment Groups",   # Custom title
  legend_title_size = 16,              # Title font size
  legend_text_size = 12,               # Legend text size
  legend_key_size = 1.0,               # Legend key size
  legend_ncol = 3,                     # Number of columns
  # P-value display options
  pvalue_format = "smart",             # numeric, scientific, smart, stars_only, combined
  pvalue_stars = TRUE,                 # Show significance stars
  pvalue_colors = TRUE,                # Color-coded significance
  pvalue_size = 12,                    # P-value text size
  pvalue_angle = 0,                    # P-value text angle
  # Pathway class annotation styling
  pathway_class_text_size = "auto",    # Auto-adjusted text size
  pathway_class_text_color = "auto",   # Theme-matched colors
  pathway_class_text_face = "bold",    # plain, bold, italic
  pathway_class_text_angle = 0         # Text rotation angle
)

# Colorblind-friendly design
p_accessible <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_annotated_results_df,
  Group = metadata$Environment,
  ko_to_kegg = TRUE,
  color_theme = "colorblind_friendly", # Colorblind-safe palette
  accessibility_mode = TRUE,           # Enhanced accessibility
  smart_colors = TRUE,                 # Intelligent color selection
  pvalue_colors = TRUE,                # Color-coded significance
  pvalue_format = "combined"           # Clear p-value display
)

# For EC, MetaCyc, KO analysis without conversions
data("metacyc_abundance")
data("metadata")
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>% column_to_rownames("pathway"), 
                                     metadata = metadata, group = "Environment", daa_method = "LinDA")
metacyc_daa_annotated_results_df <- pathway_annotation(pathway = "MetaCyc", 
                                                      daa_results_df = metacyc_daa_results_df, ko_to_kegg = FALSE)

p_metacyc <- pathway_errorbar(
  abundance = metacyc_abundance %>% column_to_rownames("pathway"),
  daa_results_df = metacyc_daa_annotated_results_df,
  Group = metadata$Environment,
  ko_to_kegg = FALSE,
  p_values_threshold = 0.05,
  order = "group",
  color_theme = "nature",              # Apply professional theme
  legend_title = "Environments",       # Descriptive legend title
  pvalue_format = "smart",             # Smart p-value formatting
  pvalue_stars = TRUE,                 # Include significance indicators
  x_lab = "description"
)

# Create example functional pathway abundance data
abundance_example <- matrix(rnorm(30), nrow = 3, ncol = 10)
colnames(abundance_example) <- paste0("Sample", 1:10)
rownames(abundance_example) <- c("PathwayA", "PathwayB", "PathwayC")

# Create example metadata
# Please change your sample id's column name to sample_name
metadata_example <- data.frame(
  sample_name = colnames(abundance_example),
  group = factor(rep(c("Control", "Treatment"), each = 5)),
  batch = factor(rep(c("Batch1", "Batch2"), times = 5))
)

# Create a basic heatmap
pathway_heatmap(abundance_example, metadata_example, "group")

# Heatmap with row clustering (pathway clustering)
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  cluster_rows = TRUE,
  clustering_method = "ward.D2",
  clustering_distance = "correlation",
  dendro_line_size = 0.8
)

# Heatmap with column clustering (sample clustering)  
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  cluster_cols = TRUE,
  clustering_method = "complete",
  clustering_distance = "euclidean"
)

# Heatmap with both row and column clustering
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  clustering_method = "average",
  clustering_distance = "manhattan"
)

# Faceted heatmap by batch
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  facet_by = "batch",
  colors = c("lightblue", "lightcoral", "lightgreen", "lightyellow")
)

# Custom colorbar settings
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  colorbar_title = "Expression Level",
  colorbar_position = "bottom",
  colorbar_width = 8,
  colorbar_height = 0.8,
  colorbar_breaks = c(-2, -1, 0, 1, 2)
)

# Left-positioned colorbar with custom colors
pathway_heatmap(
  abundance = abundance_example,
  metadata = metadata_example,
  group = "group",
  low_color = "#053061",     # Dark blue
  mid_color = "#f7f7f7",     # Light gray  
  high_color = "#67001f",    # Dark red
  colorbar_position = "left",
  colorbar_title = "Z-Score"
)

library(tidyverse)
library(ggh4x)
library(ggpicrust2)

# Load the data
data("metacyc_abundance")
data("metadata")

# Perform differential abundance analysis
metacyc_daa_results_df <- pathway_daa(
  abundance = metacyc_abundance %>% column_to_rownames("pathway"),
  metadata = metadata,
  group = "Environment",
  daa_method = "LinDA"
)

# Annotate the results
annotated_metacyc_daa_results_df <- pathway_annotation(
  pathway = "MetaCyc",
  daa_results_df = metacyc_daa_results_df,
  ko_to_kegg = FALSE
)

# Filter features with p < 0.05
feature_with_p_0.05 <- metacyc_daa_results_df %>%
  filter(p_adjust < 0.05)

# Create an advanced heatmap with clustering
pathway_heatmap(
  abundance = metacyc_abundance %>%
    right_join(
      annotated_metacyc_daa_results_df %>% select(all_of(c("feature","description"))),
      by = c("pathway" = "feature")
    ) %>%
    filter(pathway %in% feature_with_p_0.05$feature) %>%
    select(-"pathway") %>%
    column_to_rownames("description"),
  metadata = metadata,
  group = "Environment",
  cluster_rows = TRUE,
  clustering_method = "ward.D2",
  clustering_distance = "correlation",
  low_color = "#2166ac",
  mid_color = "#f7f7f7", 
  high_color = "#b2182b",
  colorbar_title = "Standardized Abundance"
)

# Create example functional pathway abundance data
abundance_example <- matrix(rnorm(30), nrow = 3, ncol = 10)
colnames(kegg_abundance_example) <- paste0("Sample", 1:10)
rownames(kegg_abundance_example) <- c("PathwayA", "PathwayB", "PathwayC")

# Create example metadata
metadata_example <- data.frame(sample_name = colnames(kegg_abundance_example),
                                group = factor(rep(c("Control", "Treatment"), each = 5)))
# Perform PCA and create visualizations
pathway_pca(abundance = abundance_example, metadata = metadata_example, "group")

# Create example functional pathway abundance data
data("metacyc_abundance")
data("metadata")

pathway_pca(abundance = metacyc_abundance %>% column_to_rownames("pathway"), metadata = metadata, group = "Environment")

library(ComplexHeatmap)
set.seed(123)
# First metagenome
metagenome1 <- abs(matrix(rnorm(1000), nrow = 100, ncol = 10))
rownames(metagenome1) <- paste0("KO", 1:100)
colnames(metagenome1) <- paste0("sample", 1:10)
# Second metagenome
metagenome2 <- abs(matrix(rnorm(1000), nrow = 100, ncol = 10))
rownames(metagenome2) <- paste0("KO", 1:100)
colnames(metagenome2) <- paste0("sample", 1:10)
# Put the metagenomes into a list
metagenomes <- list(metagenome1, metagenome2)
# Define names
names <- c("metagenome1", "metagenome2")
# Call the function
results <- compare_metagenome_results(metagenomes, names)
# Print the correlation matrix
print(results$correlation$cor_matrix)
# Print the p-value matrix
print(results$correlation$p_matrix)

library(ggpicrust2)
library(tidyverse)

# Load example data
data("ko_abundance")
data("metadata")

# Standard KEGG GSEA analysis
gsea_results_kegg <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",        # KEGG pathways (default)
  method = "fgsea",             # Fast GSEA implementation
  rank_method = "signal2noise", # Ranking method
  nperm = 1000                  # Permutations
)

# View results
head(gsea_results_kegg)

# IMPORTANT: MetaCyc GSEA requires EC (enzyme) abundance data, not pathway abundance
# The metacyc_abundance dataset contains pathway-level data, which should use pathway_daa()

# For MetaCyc GSEA, you need EC abundance from PICRUSt2:
# ec_abundance <- read_delim("EC_metagenome_out/pred_metagenome_unstrat.tsv", delim = "\t")

# Example with simulated EC abundance data for demonstration:
ec_abundance_demo <- matrix(rnorm(100 * 10), nrow = 100, ncol = 10)
rownames(ec_abundance_demo) <- paste0("EC:", sprintf("1.1.1.%d", 1:100))
colnames(ec_abundance_demo) <- paste0("Sample", 1:10)

# MetaCyc GSEA analysis (50+ metabolic pathways)
gsea_results_metacyc <- pathway_gsea(
  abundance = ec_abundance_demo,
  metadata = metadata[1:10,],   # Match sample count
  group = "Environment",
  pathway_type = "MetaCyc",      # MetaCyc pathways
  method = "fgsea",
  rank_method = "log2_ratio",
  nperm = 1000
)

# Note: If you have metacyc_abundance (pathway-level data), use pathway_daa instead:
# data("metacyc_abundance")
# daa_results <- pathway_daa(
#   abundance = metacyc_abundance %>% column_to_rownames("pathway"),
#   metadata = metadata,
#   group = "Environment",
#   daa_method = "LinDA"
# )

# GO analysis with category selection
gsea_results_go_bp <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  pathway_type = "GO",          # Gene Ontology
  go_category = "BP",           # Biological Process
  method = "fgsea",
  rank_method = "t_test",
  nperm = 1000
)

# Analyze different GO categories
for(category in c("BP", "MF", "CC")) {
  results <- pathway_gsea(
    abundance = ko_abundance %>% column_to_rownames("#NAME"),
    metadata = metadata,
    group = "Environment", 
    pathway_type = "GO",
    go_category = category,     # BP, MF, or CC
    method = "fgsea"
  )
  print(paste("GO", category, "- Significant pathways:", sum(results$p.adjust < 0.05)))
}

# Advanced GSEA with multiple ranking methods
ranking_methods <- c("signal2noise", "t_test", "log2_ratio", "diff_abundance")

gsea_comparison <- lapply(ranking_methods, function(method) {
  pathway_gsea(
    abundance = ko_abundance %>% column_to_rownames("#NAME"),
    metadata = metadata,
    group = "Environment",
    rank_method = method,        # Different ranking approaches
    min_size = 5,               # Flexible gene set sizes
    max_size = 300,
    seed = 42                   # Reproducible results
  )
})
names(gsea_comparison) <- ranking_methods

library(ggpicrust2)
library(tidyverse)

# Load example data
data("ko_abundance")
data("metadata")

# Run GSEA analysis for different pathway types
gsea_kegg <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG"
)

gsea_go <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment", 
  pathway_type = "GO",
  go_category = "BP"
)

# Annotate results for better pathway names
kegg_annotated <- gsea_pathway_annotation(gsea_kegg, pathway_type = "KEGG")
go_annotated <- gsea_pathway_annotation(gsea_go, pathway_type = "GO")

# 1. Enhanced dot plot with pathway names
dot_plot <- visualize_gsea(
  gsea_results = kegg_annotated,
  plot_type = "dotplot",
  n_pathways = 20,
  sort_by = "p.adjust",
  pathway_label_column = "pathway_name"  # Use descriptive names
)

# 2. Professional bar plot 
bar_plot <- visualize_gsea(
  gsea_results = go_annotated,
  plot_type = "barplot",
  n_pathways = 15,
  sort_by = "NES"
)

# 3. Advanced network analysis
network_plot <- visualize_gsea(
  gsea_results = kegg_annotated,
  plot_type = "network", 
  n_pathways = 20,
  network_params = list(
    similarity_measure = "jaccard",    # Pathway similarity metric
    similarity_cutoff = 0.3,           # Connection threshold  
    layout = "fruchterman",            # Network layout
    node_color_by = "NES",             # Color by enrichment score
    edge_width_by = "similarity"       # Edge width by similarity
  )
)

# 4. Expression heatmap with clustering
heatmap_plot <- visualize_gsea(
  gsea_results = kegg_annotated,
  plot_type = "heatmap",
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  n_pathways = 15,
  heatmap_params = list(
    cluster_rows = TRUE,
    cluster_columns = TRUE,
    show_rownames = TRUE,
    annotation_colors = list(
      Group = c("Aquatic" = "#2166ac", "Terrestrial" = "#762a83")
    )
  )
)

# Compare visualizations across pathway databases
kegg_plot <- visualize_gsea(kegg_annotated, plot_type = "dotplot", n_pathways = 10)
go_plot <- visualize_gsea(go_annotated, plot_type = "dotplot", n_pathways = 10)

# Combine plots for comparison
library(patchwork)
combined_plot <- kegg_plot / go_plot + 
  plot_annotation(title = "KEGG vs GO Pathway Enrichment")

# Network comparison between databases  
kegg_network <- visualize_gsea(kegg_annotated, plot_type = "network", n_pathways = 15)
go_network <- visualize_gsea(go_annotated, plot_type = "network", n_pathways = 15)

# Create enrichment plot for top pathways
top_pathways <- kegg_annotated[order(kegg_annotated$p.adjust)[1:5], ]

enrichment_plots <- lapply(1:5, function(i) {
  visualize_gsea(
    gsea_results = kegg_annotated,
    plot_type = "enrichment_plot",
    n_pathways = 1,
    sort_by = "p.adjust"
  )
})

# Multi-panel enrichment display
enrichment_combined <- wrap_plots(enrichment_plots, ncol = 2)

library(ggpicrust2)
library(tidyverse)

# Load example data
data("ko_abundance")
data("metadata")

# Perform GSEA analysis
gsea_results <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment"
)

# Perform DAA analysis
daa_results <- pathway_daa(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  daa_method = "ALDEx2"
)

# Compare GSEA and DAA results
comparison <- compare_gsea_daa(
  gsea_results = gsea_results,
  daa_results_df = daa_results,
  gsea_pvalue_cutoff = 0.05,
  daa_pvalue_cutoff = 0.05,
  plot_type = "venn"  # Can be "venn", "upset", or "both"
)

# View the comparison plot
comparison$plot

# View the overlapping pathways
head(comparison$overlap)

library(ggpicrust2)
library(tidyverse)

# Load example data and perform GSEA
data("ko_abundance")
data("metadata")

gsea_results <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment"
)

# Annotate GSEA results
annotated_results <- gsea_pathway_annotation(
  gsea_results = gsea_results,
  pathway = "KO"
)

# View the annotated results
head(annotated_results)

pathway_errorbar(abundance = abundance,
                 daa_results_df = daa_results_df,
                 Group = metadata$Environment,
                 ko_to_kegg = TRUE,
                 p_values_threshold = 0.05,
                 order = "pathway_class",
                 select = NULL,
                 p_value_bar = TRUE,
                 colors = NULL,
                 x_lab = "pathway_name")

library(patchwork)

library(ggprism)