R/pathway_errorbar.R

In ggpicrust2: Make 'PICRUSt2' Output Analysis and Visualization Easier

#' The function pathway_errorbar() is used to visualize the results of functional pathway differential abundance analysis as error bar plots.
#'
#' @name pathway_errorbar
#' @param abundance A data frame with row names representing pathways and column names representing samples. Each element represents the relative abundance of the corresponding pathway in the corresponding sample.
#' @param daa_results_df A data frame containing the results of the differential abundance analysis of the pathways, generated by the pathway_daa function. x_lab should be a column name of daa_results_df.
#' @param Group A data frame or a vector that assigns each sample to a group. The groups are used to color the samples in the figure.
#' @param ko_to_kegg A logical parameter indicating whether there was a convertion that convert ko abundance to kegg abundance.
#' @param p_values_threshold A numeric parameter specifying the threshold for statistical significance of differential abundance. Pathways with p-values below this threshold will be considered significant.
#' @param order A parameter controlling the ordering of the rows in the figure. The options are: "p_values" (order by p-values), "name" (order by pathway name), "group" (order by the group with the highest mean relative abundance), or "pathway_class" (order by the pathway category).
#' @param select A vector of pathway names to be included in the figure. This can be used to limit the number of pathways displayed. If NULL, all pathways will be displayed.
#' @param p_value_bar A logical parameter indicating whether to display a bar showing the p-value threshold for significance. If TRUE, the bar will be displayed.
#' @param colors A vector of colors to be used to represent the groups in the figure. Each color corresponds to a group. If NULL, colors will be selected based on the color_theme.
#' @param x_lab A character string to be used as the x-axis label in the figure. The default value is "description" for KOs'descriptions and "pathway_name" for KEGG pathway names.
#' @param log2_fold_change_color A character string specifying the color for log2 fold change bars. Default is "#87ceeb" (light blue). Can also be "auto" to use theme-based colors.
#' @param max_features A numeric parameter specifying the maximum number of features to display before issuing a warning. Default is 30. Set to a higher value to display more features, or Inf to disable the limit entirely.
#' @param color_theme A character string specifying the color theme to use. Options include: "default", "nature", "science", "cell", "nejm", "lancet", "colorblind_friendly", "viridis", "plasma", "minimal", "high_contrast", "pastel", "bold". Default is "default".
#' @param pathway_class_colors A vector of colors for pathway class annotations. If NULL, colors will be selected from the theme.
#' @param smart_colors A logical parameter indicating whether to use intelligent color selection based on data characteristics. Default is FALSE.
#' @param accessibility_mode A logical parameter indicating whether to use accessibility-friendly colors. Default is FALSE.
#' @param legend_position A character string specifying legend position. Options: "top", "bottom", "left", "right", "none". Default is "top".
#' @param legend_direction A character string specifying legend direction. Options: "horizontal", "vertical". Default is "horizontal".
#' @param legend_title A character string for legend title. If NULL, no title is displayed.
#' @param legend_title_size A numeric value specifying legend title font size. Default is 12.
#' @param legend_text_size A numeric value specifying legend text font size. Default is 10.
#' @param legend_key_size A numeric value specifying legend key size in cm. Default is 0.8.
#' @param legend_ncol A numeric value specifying number of columns in legend. If NULL, automatic layout is used.
#' @param legend_nrow A numeric value specifying number of rows in legend. If NULL, automatic layout is used.
#' @param pvalue_format A character string specifying p-value format. Options: "numeric", "scientific", "smart", "stars_only", "combined". Default is "smart".
#' @param pvalue_stars A logical parameter indicating whether to display significance stars. Default is TRUE.
#' @param pvalue_colors A logical parameter indicating whether to use color coding for significance levels. Default is FALSE.
#' @param pvalue_size A numeric value or "auto" for p-value text size. Default is "auto".
#' @param pvalue_angle A numeric value specifying p-value text angle in degrees. Default is 0.
#' @param pvalue_thresholds A numeric vector of significance thresholds. Default is c(0.001, 0.01, 0.05).
#' @param pvalue_star_symbols A character vector of star symbols for significance levels. Default is c("***", "**", "*").
#' @param pathway_class_text_size A numeric value or "auto" for pathway class text size. Default is "auto".
#' @param pathway_class_text_color A character string for pathway class text color. Use "auto" for theme-based color. Default is "black".
#' @param pathway_class_text_face A character string for pathway class text face. Options: "plain", "bold", "italic". Default is "bold".
#' @param pathway_class_text_angle A numeric value specifying pathway class text angle in degrees. Default is 0.
#' @param pathway_class_position A character string specifying pathway class position. Options: "left", "right", "none". Default is "left".
#' @param pathway_names_text_size A numeric value or "auto" for pathway names (y-axis labels) text size. Default is "auto".
#' @importFrom stats sd
#' @return A ggplot2 plot showing the error bar plot of the differential abundance analysis results for the functional pathways.
#' The plot visualizes the differential abundance results of a specific differential abundance analysis method. The corresponding dataframe contains the results used to create the plot.
#' @export
#'
#' @examples
#' \dontrun{
#' # Example 1: Analyzing KEGG pathway abundance
#' metadata <- read_delim(
#'   "path/to/your/metadata.txt",
#'   delim = "\t",
#'   escape_double = FALSE,
#'   trim_ws = TRUE
#' )
#'
#' # data(metadata)
#'
#' kegg_abundance <- ko2kegg_abundance(
#'   "path/to/your/pred_metagenome_unstrat.tsv"
#' )
#'
#' # data(kegg_abundance)
#'
#' # Please change group to "your_group_column" if you are not using example dataset
#' group <- "Environment"
#'
#' daa_results_df <- pathway_daa(
#'   abundance = kegg_abundance,
#'   metadata = metadata,
#'   group = group,
#'   daa_method = "ALDEx2",
#'   select = NULL,
#'   reference = NULL
#' )
#'
#' # Please check the unique(daa_results_df$method) and choose one
#' daa_sub_method_results_df <- daa_results_df[daa_results_df$method
#' == "ALDEx2_Welch's t test", ]
#'
#' daa_annotated_sub_method_results_df <- pathway_annotation(
#'   pathway = "KO",
#'   daa_results_df = daa_sub_method_results_df,
#'   ko_to_kegg = TRUE
#' )
#'
#' # Please change Group to metadata$your_group_column if you are not using example dataset
#' Group <- metadata$Environment
#'
#' p <- pathway_errorbar(
#'   abundance = kegg_abundance,
#'   daa_results_df = daa_annotated_sub_method_results_df,
#'   Group = Group,
#'   p_values_threshold = 0.05,
#'   order = "pathway_class",
#'   select = daa_annotated_sub_method_results_df %>%
#'   arrange(p_adjust) %>%
#'   slice(1:20) %>%
#'   select("feature") %>% pull(),
#'   ko_to_kegg = TRUE,
#'   p_value_bar = TRUE,
#'   colors = NULL,
#'   x_lab = "pathway_name",
#'   log2_fold_change_color = "#FF5733" # Custom color for log2 fold change bars
#' )
#'
#' # Example 2: Analyzing EC, MetaCyc, KO without conversions
#' metadata <- read_delim(
#'   "path/to/your/metadata.txt",
#'   delim = "\t",
#'   escape_double = FALSE,
#'   trim_ws = TRUE
#' )
#' # data(metadata)
#'
#' metacyc_abundance <- read.delim("path/to/your/metacyc_abundance.tsv")
#'
#' # data(metacyc_abundance)
#'
#' group <- "Environment"
#'
#' daa_results_df <- pathway_daa(
#'   abundance = metacyc_abundance %>% column_to_rownames("pathway"),
#'   metadata = metadata,
#'   group = group,
#'   daa_method = "LinDA",
#'   select = NULL,
#'   reference = NULL
#' )
#'
#'
#' daa_annotated_results_df <- pathway_annotation(
#'   pathway = "MetaCyc",
#'   daa_results_df = daa_results_df,
#'   ko_to_kegg = FALSE
#' )
#'
#' Group <- metadata$Environment
#'
#' p <- pathway_errorbar(
#'   abundance = metacyc_abundance %>% column_to_rownames("pathway"),
#'   daa_results_df = daa_annotated_results_df,
#'   Group = Group,
#'   p_values_threshold = 0.05,
#'   order = "group",
#'   select = NULL,
#'   ko_to_kegg = FALSE,
#'   p_value_bar = TRUE,
#'   colors = NULL,
#'   x_lab = "description",
#'   log2_fold_change_color = "#006400" # Dark green for log2 fold change bars
#' )
#' }
utils::globalVariables(c("group", "name", "value", "feature", "negative_log10_p", "group_nonsense", "nonsense", "pathway_class", "p_adjust", "log_2_fold_change", "transform_sample_counts", "column_to_rownames", "txtProgressBar", "setTxtProgressBar", "utils"))
pathway_errorbar <-
  function(abundance,
           daa_results_df,
           Group,
           ko_to_kegg = FALSE,
           p_values_threshold = 0.05,
           order = "group",
           select = NULL,
           p_value_bar = TRUE,
           colors = NULL,
           x_lab = NULL,
           log2_fold_change_color = "#87ceeb",
           max_features = 30,
           color_theme = "default",
           pathway_class_colors = NULL,
           smart_colors = FALSE,
           accessibility_mode = FALSE,
           # New legend control parameters
           legend_position = "top",
           legend_direction = "horizontal",
           legend_title = NULL,
           legend_title_size = 12,
           legend_text_size = 10,
           legend_key_size = 0.8,
           legend_ncol = NULL,
           legend_nrow = NULL,
           # P-value display parameters
           pvalue_format = "smart",
           pvalue_stars = TRUE,
           pvalue_colors = FALSE,
           pvalue_size = "auto",
           pvalue_angle = 0,
           pvalue_thresholds = c(0.001, 0.01, 0.05),
           pvalue_star_symbols = c("***", "**", "*"),
           # Pathway class annotation parameters
           pathway_class_text_size = "auto",
           pathway_class_text_color = "black",
           pathway_class_text_face = "bold",
           pathway_class_text_angle = 0,
           pathway_class_position = "left",
           # Pathway names text size parameter
           pathway_names_text_size = "auto") {
    # Add more complete input validation at the beginning of the function
    if(!is.matrix(abundance) && !is.data.frame(abundance)) {
      stop("'abundance' must be a matrix or data frame")
    }

    if(!is.data.frame(daa_results_df)) {
      stop("'daa_results_df' must be a data frame")
    }

    # Check required columns
    required_cols <- c("feature", "method", "group1", "group2", "p_adjust")
    missing_cols <- setdiff(required_cols, colnames(daa_results_df))
    if(length(missing_cols) > 0) {
      stop("Missing required columns in daa_results_df: ",
           paste(missing_cols, collapse = ", "))
    }

    # Add data validation at the beginning of the function
    if (length(Group) != ncol(abundance)) {
      stop("Length of Group must match number of columns in abundance matrix")
    }

    # Check the number of significant features
    # Calculate number of significant features (for reference)
    # sig_features <- sum(daa_results_df$p_adjust < 0.05)

    # Identify pathways with missing annotation
    missing_pathways <- daa_results_df[is.na(daa_results_df$pathway_name), "feature"]

    # Inform the user about the missing annotations
    if(length(missing_pathways) > 0) {
      message("The following pathways are missing annotations and have been excluded: ",
              paste(missing_pathways, collapse = ", "))
      message("You can use the 'pathway_annotation' function to add annotations for these pathways.")
    }

    # Get the names of all columns in the data frame
    column_names <- colnames(daa_results_df)

    # Check if there are more than two 'group' columns
    group_columns <- grepl("^group", column_names)
    if (sum(group_columns) > 2) {
      p_value_bar <- FALSE
      message(
        "There are more than two 'group' columns in the 'daa_results_df' data frame. As a result, it is not possible to compute the log2 fold values. The 'p_value_bar' has been automatically set to FALSE."
      )
    }

    # Set x_lab if not provided
    if (is.null(x_lab)){
      if (ko_to_kegg == TRUE){
        x_lab <- "pathway_name"
      }else{
        x_lab <- "description"
      }

      if (is.null(daa_results_df$pathway_name) && is.null(daa_results_df$description)){
        message(
          "Please utilize the 'pathway_annotation' function to annotate the 'daa_results_df' data frame."
        )
      }
    }

    if (!(x_lab %in% colnames(daa_results_df))){
      message(
        "The 'x_lab' you defined does not exist as a column in the 'daa_results_df' data frame."
      )
    }

    # Exclude rows with missing pathway annotation (after x_lab is set)
    original_rows <- nrow(daa_results_df)
    daa_results_df <- daa_results_df[!is.na(daa_results_df[,x_lab]),]
    filtered_rows <- nrow(daa_results_df)
    
    # Check if all rows were filtered out due to missing annotations
    if (filtered_rows == 0) {
      warning(
        sprintf("All %d rows were excluded due to missing '%s' annotations. ", original_rows, x_lab),
        "This typically occurs when no statistically significant pathways were found. ",
        "Cannot create error bar plot with empty data. Returning NULL.",
        call. = FALSE
      )
      return(NULL)
    }
    
    if (original_rows > filtered_rows) {
      message(sprintf("Excluded %d rows with missing '%s' annotations.", original_rows - filtered_rows, x_lab))
    }

    if (length(unique(daa_results_df$method)) != 1) {
      stop(
        "The 'method' column in the 'daa_results_df' data frame contains more than one method. Please filter it to contain only one method."
      )
    }

    if (length(unique(daa_results_df$group1)) != 1 || length(unique(daa_results_df$group2)) != 1) {
      stop(
        "The 'group1' or 'group2' column in the 'daa_results_df' data frame contains more than one group. Please filter each to contain only one group."
      )
    }

    # Enhanced color selection using the new color theme system
    n_groups <- nlevels(as.factor(Group))
    
    # Source the color themes and legend utilities if functions are not available
    if (!exists("get_color_theme")) {
      source(file.path(dirname(parent.frame()$ofile), "color_themes.R"))
    }
    if (!exists("format_pvalue_smart")) {
      source(file.path(dirname(parent.frame()$ofile), "legend_annotation_utils.R"))
    }
    
    # Use smart color selection if requested
    if (smart_colors || accessibility_mode) {
      smart_selection <- smart_color_selection(
        n_groups = n_groups,
        has_pathway_class = ko_to_kegg,
        data_type = "abundance",
        accessibility_mode = accessibility_mode
      )
      color_theme <- smart_selection$theme_name
      if (smart_colors) {
        message("Smart color selection: ", smart_selection$reason)
      }
    }
    
    # Get the color theme
    theme_colors <- get_color_theme(color_theme, n_groups)
    
    # Set group colors
    if (is.null(colors)) {
      colors <- theme_colors$group_colors[1:n_groups]
    }
    
    # Set pathway class colors if not provided
    if (is.null(pathway_class_colors)) {
      pathway_class_colors <- theme_colors$pathway_class_colors
    }
    
    # Set log2 fold change color
    if (log2_fold_change_color == "auto") {
      log2_fold_change_color <- theme_colors$fold_change_single
    }

    errorbar_abundance_mat <- as.matrix(abundance)

    daa_results_filtered_df <-
      daa_results_df[daa_results_df$p_adjust < p_values_threshold,]

    if (!is.null(select)) {
      daa_results_filtered_sub_df <-
        daa_results_filtered_df[daa_results_filtered_df$feature %in% select, ]
    } else {
      daa_results_filtered_sub_df <- daa_results_filtered_df
    }

    if (nrow(daa_results_filtered_sub_df) > max_features) {
      message(
        paste0(
          "The number of features with statistical significance exceeds ", max_features, ", which may lead to suboptimal visualization. ",
          "Please use 'select' to reduce the number of features or increase the 'max_features' parameter.\n",
          "Currently, you have these features: ",
          paste(paste0('"', daa_results_filtered_sub_df$feature, '"'), collapse = ", "), ".\n",
          "You can find the statistically significant features with the following command:\n",
          "daa_results_df %>% filter(p_adjust < 0.05) %>% select(c(\"feature\",\"p_adjust\"))"
        )
      )
      warning("The number of features with statistical significance exceeds ", max_features, ", which may lead to suboptimal visualization. Setting max_features=Inf will disable this warning.")
      # Changed from stop() to warning() to allow users to proceed with visualization if desired
    }

    if (nrow(daa_results_filtered_sub_df) == 0){
      stop(
        "Visualization with 'pathway_errorbar' cannot be performed because there are no features with statistical significance. ",
        "For possible solutions, please check the FAQ section of the tutorial."
      )
    }
    # Convert to relative abundance
    relative_abundance_mat <- apply(t(errorbar_abundance_mat), 1, function(x)
      x / sum(x))

    # Subset to only include the features present in daa_results_filtered_sub_df$feature
    sub_relative_abundance_mat <- relative_abundance_mat[rownames(relative_abundance_mat) %in% daa_results_filtered_sub_df$feature, , drop = FALSE]

    # Create a matrix for the error bars
    # Fix mapping error: ensure correct mapping between samples and their experimental groups
    if (length(Group) != ncol(abundance)) {
      stop("Length of Group must match number of columns in abundance matrix")
    }
    
    # Create a mapping from sample names to their corresponding groups
    # Note: We can't directly use the Group vector as its order may not match the column names in sub_relative_abundance_mat
    sample_to_group_map <- stats::setNames(as.character(Group), colnames(abundance))
    
    # Get the correct groups based on the column names in sub_relative_abundance_mat
    correct_groups <- sample_to_group_map[colnames(sub_relative_abundance_mat)]
    
    # Create the error bar matrix using the correctly mapped groups
    error_bar_matrix <- cbind(
      sample = colnames(sub_relative_abundance_mat),
      group = correct_groups,
      t(sub_relative_abundance_mat)
    )
    error_bar_df <- as.data.frame(error_bar_matrix)

    error_bar_df$group <- factor(correct_groups, levels = levels(as.factor(Group)))

      error_bar_pivot_longer_df <- tidyr::pivot_longer(error_bar_df,-c(sample, group))

    error_bar_pivot_longer_tibble <-
      mutate(error_bar_pivot_longer_df, group = as.factor(group))

    error_bar_pivot_longer_tibble$sample <-
      factor(error_bar_pivot_longer_tibble$sample)

    error_bar_pivot_longer_tibble$name <-
      factor(error_bar_pivot_longer_tibble$name)

    error_bar_pivot_longer_tibble$value <-
      as.numeric(error_bar_pivot_longer_tibble$value)

    error_bar_pivot_longer_tibble_summarised <-
      error_bar_pivot_longer_tibble %>% group_by(name, group) %>%
      summarise(mean = mean(value), sd = stats::sd(value))

    error_bar_pivot_longer_tibble_summarised <-
      error_bar_pivot_longer_tibble_summarised %>% mutate(group2 = "nonsense")

    switch(
      order,
      "p_values" = {
        order <- order(daa_results_filtered_sub_df$p_adjust)
      },
      "name" = {
        order <- order(daa_results_filtered_sub_df$feature)
      },
      "group" = {
        # Initialize pro column with default value 1
        daa_results_filtered_sub_df$pro <- 1
        
        for (i in levels(error_bar_pivot_longer_tibble_summarised$name)) {
          # Get subset for current feature
          error_bar_pivot_longer_tibble_summarised_sub <-
            error_bar_pivot_longer_tibble_summarised[error_bar_pivot_longer_tibble_summarised$name == i,]
          
          # Find group with maximum mean abundance
          pro_group <-
            error_bar_pivot_longer_tibble_summarised_sub[error_bar_pivot_longer_tibble_summarised_sub$mean ==
                                                            max(error_bar_pivot_longer_tibble_summarised_sub$mean),]$group
          pro_group <- as.vector(pro_group)
          
          # Find indices of rows matching the current feature, excluding NA values
          idx <- which(daa_results_filtered_sub_df$feature == i & !is.na(daa_results_filtered_sub_df$feature))
          
          # Only assign values if valid indices exist
          if (length(idx) > 0) {
            daa_results_filtered_sub_df$pro[idx] <- pro_group
          }
        }
        
        # Order by group and p-value
        order <-
          order(daa_results_filtered_sub_df$pro,
                daa_results_filtered_sub_df$p_adjust)
      },
      "pathway_class" = {
        if (!"pathway_class" %in% colnames(daa_results_filtered_sub_df)) {
          stop(
            "The 'pathway_class' column is missing in the 'daa_results_filtered_sub_df' data frame. ",
            "Please use the 'pathway_annotation' function to annotate the 'pathway_daa' results."
          )
        }
        order <- order(
          daa_results_filtered_sub_df$pathway_class,
          daa_results_filtered_sub_df$p_adjust
        )
      },
      {
        order <- order
      }
    )

    daa_results_filtered_sub_df <-
      daa_results_filtered_sub_df[order,]

    error_bar_pivot_longer_tibble_summarised_ordered <-
      data.frame(
        name = NULL,
        group = NULL,
        mean = NULL,
        sd = NULL
      )

    for (i in daa_results_filtered_sub_df$feature) {
      error_bar_pivot_longer_tibble_summarised_ordered <-
        rbind(
          error_bar_pivot_longer_tibble_summarised_ordered,
          error_bar_pivot_longer_tibble_summarised[error_bar_pivot_longer_tibble_summarised$name ==
                                                     i,]
        )
    }

    if (ko_to_kegg == FALSE) {
      # Match by feature name using the match function
      matched_indices <- match(
        error_bar_pivot_longer_tibble_summarised_ordered$name,
        daa_results_filtered_sub_df$feature
      )
      error_bar_pivot_longer_tibble_summarised_ordered[, x_lab] <-
        daa_results_filtered_sub_df[matched_indices, x_lab]
    }

    if (ko_to_kegg == TRUE) {
      error_bar_pivot_longer_tibble_summarised_ordered$pathway_class <-
        rep(daa_results_filtered_sub_df$pathway_class,
            each = length(levels(
              factor(error_bar_pivot_longer_tibble_summarised_ordered$group)
            )))
    }

    error_bar_pivot_longer_tibble_summarised_ordered$name <- factor(error_bar_pivot_longer_tibble_summarised_ordered$name, levels = rev(daa_results_filtered_sub_df$feature))

    # Calculate smart text size for pathway names (y-axis labels)
    pathway_names_final_text_size <- if (pathway_names_text_size == "auto") {
      if (exists("calculate_smart_text_size")) {
        calculate_smart_text_size(nrow(daa_results_filtered_sub_df),
                                base_size = 10, min_size = 8, max_size = 14)
      } else {
        10  # Default size (current hardcoded value)
      }
    } else {
      pathway_names_text_size
    }

    # Calculate smart text size for pathway class annotations
    pathway_class_final_text_size <- if (pathway_class_text_size == "auto") {
      if (exists("calculate_smart_text_size")) {
        calculate_smart_text_size(nrow(daa_results_filtered_sub_df),
                                base_size = 10, min_size = 3, max_size = 4)
      } else {
        3.5  # Default size
      }
    } else {
      pathway_class_text_size
    }

    # Set pathway class text color based on theme if "auto"
    pathway_class_final_text_color <- if (pathway_class_text_color == "auto") {
      # Get theme colors for pathway class
      current_theme <- get_color_theme(color_theme)
      current_theme$pathway_class_colors[1]  # Use first theme color
    } else {
      pathway_class_text_color
    }

    bar_errorbar <-
      ggplot2::ggplot(error_bar_pivot_longer_tibble_summarised_ordered,
             ggplot2::aes(mean, name, fill = group)) +
      ggplot2::geom_errorbar(
        ggplot2::aes(xmax = mean + sd, xmin = 0),
        position = ggplot2::position_dodge(width = 0.8),
        width = 0.5,
        linewidth = 0.5,
        color = "black"
      ) +
      ggplot2::geom_bar(stat = "identity",
               position = ggplot2::position_dodge(width = 0.8),
               width = 0.8) +
      GGally::geom_stripped_cols(width = 10) +
      ggplot2::scale_fill_manual(values = colors) +
      ggplot2::scale_color_manual(values = colors) +
      ggprism::theme_prism(base_size = 12) +
      ggplot2::scale_x_continuous(expand = c(0, 0)) +
      ggplot2::scale_y_discrete(labels = rev(daa_results_filtered_sub_df[, x_lab])) +
      ggplot2::labs(x = "Relative Abundance", y = NULL) +
      ggplot2::theme(
        axis.ticks.y = ggplot2::element_blank(),
        axis.line.y = ggplot2::element_blank(),
        axis.line.x = ggplot2::element_line(linewidth = 0.5),
        axis.ticks.x = ggplot2::element_line(linewidth = 0.5),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        axis.text = ggplot2::element_text(size = 10, color = "black"),
        axis.text.x = ggplot2::element_text(margin = ggplot2::margin(r = 0)),
        axis.text.y = ggplot2::element_text(
          size = pathway_names_final_text_size,
          color = "black",
          margin = ggplot2::margin(b = 6)
        ),
        axis.title.x =  ggplot2::element_text(
          size = 10,
          color = "black",
          hjust = 0.5
        ),
        legend.position = legend_position,
        legend.key.size = ggplot2::unit(legend_key_size, "cm"),
        legend.direction = legend_direction,
        legend.justification = "center",
        legend.text = ggplot2::element_text(size = legend_text_size, face = "bold"),
        legend.box.just = "center",
        legend.title = if (!is.null(legend_title)) {
          ggplot2::element_text(size = legend_title_size, face = "bold")
        } else {
          ggplot2::element_blank()
        },
        plot.margin = ggplot2::margin(0, 0.5, 0.5, 0, unit = "cm")
      ) +
      # Add guide specifications for multi-column/row layouts
      {if (!is.null(legend_ncol) || !is.null(legend_nrow)) {
        ggplot2::guides(fill = ggplot2::guide_legend(
          ncol = legend_ncol,
          nrow = legend_nrow,
          byrow = TRUE
        ))
      } else {
        ggplot2::guides()
      }} +
      ggplot2::coord_cartesian(clip = "off") +
      ggplot2::theme(
        axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
        plot.margin = ggplot2::unit(c(1, 1, 1, 1), "cm")
      )

    if (ko_to_kegg == TRUE) {
      pathway_class_group_mat <-
        daa_results_filtered_sub_df$pathway_class %>%
        table() %>%
        data.frame() %>% tibble::column_to_rownames(".")
      pathway_class_group <- data.frame(.= unique(daa_results_filtered_sub_df$pathway_class),Freq = pathway_class_group_mat[unique(daa_results_filtered_sub_df$pathway_class),])
      start <-
        c(1, rev(pathway_class_group$Freq)[1:(length(pathway_class_group$Freq) - 1)]) %>%
        cumsum()
      end <- cumsum(rev(pathway_class_group$Freq))
      ymin <- start - 1 / 2
      ymax <- end + 1 / 2
      nPoints <- length(start)
      pCol <- pathway_class_colors[1:nPoints]
      pFill <- pCol
      for (i in 1:nPoints)  {
        bar_errorbar <- bar_errorbar +
          ggplot2::annotation_custom(
            grob = grid::rectGrob(
              gp = grid::gpar(
                col = pCol[i],
                fill = pFill[i],
                lty = NULL,
                lwd = NULL,
                alpha = 0.2
              )
            ),
            xmin = ggplot2::unit(-2, 'native'),
            xmax = ggplot2::unit(0, 'native'),
            ymin = ggplot2::unit(ymin[i], 'native'),
            ymax = ggplot2::unit(ymax[i], 'native')
          )
      }
    }
    # Add necessary columns
    daa_results_filtered_sub_df$negative_log10_p <- -log10(daa_results_filtered_sub_df$p_adjust)
    daa_results_filtered_sub_df$group_nonsense <- "nonsense"
    
    # Only add log_2_fold_change if it doesn't exist
    if (!"log_2_fold_change" %in% colnames(daa_results_filtered_sub_df)) {
      daa_results_filtered_sub_df$log_2_fold_change <- rep(NA, nrow(daa_results_filtered_sub_df))
    }

    for (i in daa_results_filtered_sub_df$feature){
      # Get mean values for this feature
      feature_means <- error_bar_pivot_longer_tibble_summarised_ordered[error_bar_pivot_longer_tibble_summarised_ordered$name %in% i,]

      # Get group1 and group2 names for this feature from DAA results
      feature_row <- daa_results_filtered_sub_df[daa_results_filtered_sub_df$feature == i, ]
      group1_name <- feature_row$group1[1]
      group2_name <- feature_row$group2[1]

      # Get mean values for each group in the correct order
      mean_group1 <- feature_means[feature_means$group == group1_name, ]$mean
      mean_group2 <- feature_means[feature_means$group == group2_name, ]$mean

      # Calculate log2 fold change as group2/group1 (consistent with calculate_abundance_stats)
      if (length(mean_group1) > 0 && length(mean_group2) > 0) {
        # Add small pseudocount to avoid log(0)
        pseudocount <- 1e-10
        log2_fc <- log2((mean_group2 + pseudocount) / (mean_group1 + pseudocount))
        daa_results_filtered_sub_df[daa_results_filtered_sub_df$feature==i,]$log_2_fold_change <- log2_fc
      } else {
        # Fallback to original calculation if group matching fails
        mean <- feature_means$mean
        daa_results_filtered_sub_df[daa_results_filtered_sub_df$feature==i,]$log_2_fold_change <- log2(mean[1]/mean[2])
      }
    }
    daa_results_filtered_sub_df$feature <- factor(daa_results_filtered_sub_df$feature,levels = rev(daa_results_filtered_sub_df$feature))
    p_values_bar <- daa_results_filtered_sub_df %>%
      ggplot2::ggplot(ggplot2::aes(feature, log_2_fold_change, fill = group_nonsense)) +
      ggplot2::geom_bar(stat = "identity",
               position = ggplot2::position_dodge(width = 0.8),
               width = 0.8) +
      ggplot2::labs(y = "log2 fold change", x = NULL) +
      GGally::geom_stripped_cols() +
      ggplot2::scale_fill_manual(values = log2_fold_change_color) +
      ggplot2::scale_color_manual(values = log2_fold_change_color) +
      ggplot2::geom_hline(ggplot2::aes(yintercept = 0),
                 linetype = 'dashed',
                 color = 'black') +
      ggprism::theme_prism(base_size = 12) +
      ggplot2::scale_y_continuous(expand = c(0, 0)) +
      ggplot2::theme(
        axis.ticks.y = ggplot2::element_blank(),
        axis.line.y = ggplot2::element_blank(),
        axis.line.x = ggplot2::element_line(linewidth = 0.5),
        axis.ticks.x = ggplot2::element_line(linewidth = 0.5),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        axis.text = ggplot2::element_text(size = 10, color = "black"),
        axis.text.y = ggplot2::element_blank(),
        axis.text.x = ggplot2::element_text(
          size = 10,
          color = "black",
          margin = ggplot2::margin(b = 6)
        ),
        axis.title.x =  ggplot2::element_text(
          size = 11,
          color = "black",
          hjust = 0.5
        ),
        legend.position = "non"
      ) +
      ggplot2::coord_flip()

    if (ko_to_kegg == TRUE) {
      pathway_class_y <- (ymax + ymin) / 2 - 0.5
      pathway_class_plot_df <-
        as.data.frame(
          cbind(
            nonsense = "nonsense",
            pathway_class_y = pathway_class_y,
            pathway_class = rev(unique(
              daa_results_filtered_sub_df$pathway_class
            ))
          )
        )
      pathway_class_plot_df$pathway_class_y <-
        as.numeric(pathway_class_plot_df$pathway_class_y)
      pathway_class_annotation <-
        pathway_class_plot_df %>% ggplot2::ggplot(ggplot2::aes(nonsense, pathway_class_y)) + ggplot2::geom_text(
          ggplot2::aes(nonsense, pathway_class_y, label = pathway_class),
          size = pathway_class_final_text_size,
          color = pathway_class_final_text_color,
          fontface = pathway_class_text_face,
          family = "sans",
          angle = pathway_class_text_angle,
          hjust = if (pathway_class_position == "left") 1 else 0
        ) +
        ggplot2::scale_y_discrete(position = "right") +
        ggprism::theme_prism(base_size = 12) +
        ggplot2::theme(
          axis.ticks = ggplot2::element_blank(),
          axis.line = ggplot2::element_blank(),
          panel.grid.major.y = ggplot2::element_blank(),
          panel.grid.major.x = ggplot2::element_blank(),
          panel.background = ggplot2::element_blank(),
          axis.text = ggplot2::element_blank(),
          plot.margin = ggplot2::unit(c(0, 0.2, 0, 0), "cm"),
          axis.title.y =  ggplot2::element_blank(),
          axis.title.x = ggplot2::element_blank(),
          legend.position = "non"
        )
    }

    # Enhanced p-value formatting using the new system
    format_p_value <- function(p) {
      # Use smart formatting if available, otherwise fallback to simple formatting
      if (exists("format_pvalue_smart")) {
        format_pvalue_smart(p, 
                           format = pvalue_format,
                           stars = pvalue_stars,
                           thresholds = pvalue_thresholds,
                           star_symbols = pvalue_star_symbols)
      } else {
        # Fallback to original formatting
        ifelse(p < 0.001, sprintf("%.1e", p), sprintf("%.3f", p))
      }
    }
    



    # Calculate p-value colors if enabled
    pvalue_text_colors <- if (pvalue_colors && exists("get_significance_colors")) {
      get_significance_colors(daa_results_filtered_sub_df$p_adjust,
                            thresholds = pvalue_thresholds,
                            colors = c("#d73027", "#fc8d59", "#fee08b"),
                            default_color = "black")
    } else {
      rep("black", nrow(daa_results_filtered_sub_df))
    }
    
    # Calculate smart text size for p-values
    pvalue_text_size <- if (pvalue_size == "auto") {
      if (exists("calculate_smart_text_size")) {
        calculate_smart_text_size(nrow(daa_results_filtered_sub_df), base_size = 10)
      } else {
        3.5  # Default size
      }
    } else {
      pvalue_size
    }

    daa_results_filtered_sub_df$unique <-
      nrow(daa_results_filtered_sub_df) - seq_len(nrow(daa_results_filtered_sub_df)) + 1
    p_annotation <- daa_results_filtered_sub_df %>%
      ggplot2::ggplot(ggplot2::aes(group_nonsense, p_adjust)) +
      ggplot2::geom_text(
        ggplot2::aes(group_nonsense, unique, 
                    label = format_p_value(p_adjust),
                    color = I(pvalue_text_colors)),
        size = pvalue_text_size,
        fontface = "bold",
        family = "sans",
        angle = pvalue_angle
      ) +
      ggplot2::labs(y = "p-value (adjusted)") +
      ggplot2::scale_y_discrete(position = "right") +
      ggprism::theme_prism(base_size = 12) +
      ggplot2::theme(
        axis.ticks = ggplot2::element_blank(),
        axis.line = ggplot2::element_blank(),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        panel.background = ggplot2::element_blank(),
        axis.text = ggplot2::element_blank(),
        plot.margin = ggplot2::unit(c(0, 0.2, 0, 0), "cm"),
        axis.title.y =  ggplot2::element_text(
          size = 11,
          color = "black",
          vjust = 0
        ),
        axis.title.x = ggplot2::element_blank(),
        legend.position = "non"
      )
    if (p_value_bar == TRUE) {
      if (ko_to_kegg == TRUE) {
        combination_bar_plot <-
          pathway_class_annotation + bar_errorbar + p_values_bar + p_annotation + patchwork::plot_layout(ncol = 4, widths =
                                                                                                c(1, 1.2, 0.5, 0.1))
      }
      else{
        combination_bar_plot <-
          bar_errorbar + p_values_bar + p_annotation + patchwork::plot_layout(ncol = 3, widths = c(2.3, 0.7, 0.3))
      }
    }else{
      if (ko_to_kegg == TRUE) {
        combination_bar_plot <-
          pathway_class_annotation + bar_errorbar + p_annotation + patchwork::plot_layout(ncol = 3, widths =
                                                                                                           c(1, 1.2, 0.1))
      }
      else{
        combination_bar_plot <-
          bar_errorbar + p_annotation + patchwork::plot_layout(ncol = 2, widths = c(2.5,  0.2))
      }
    }
    return(combination_bar_plot)
  }