GGpicrust2 Pathway_Errorbar Chart Discrepancy - Urgent Assistance Required

[MizaR108]

Dear Caffery Yang,

I'm writing to you today to express my sincere gratitude for the incredible wealth of knowledge and insights I've gained from the code you've generously shared. Your contributions have been a tremendous asset to my learning journey, and I am truly grateful for your dedication to the field.

I'm currently working on a ggpicrust2 analysis using Metacyc data, and I've encountered a rather perplexing issue with the pathway_errorbar chart. The bar chart's representation appears to be misaligned with the actual abundance values.

To illustrate, I manually calculated the group-wise mean abundance for specific pathways and observed clear visual differences. However, the pathway_errorbar chart often displays contradictory results or significantly underestimates the magnitude of these differences. This discrepancy is concerning as it could potentially lead to erroneous conclusions.

It's worth noting that my abundance data is represented as proportions for each sample. Could this be the source of the problem?

To help you better understand the issue, I've attached the following:

Comparative images: These images illustrate the discrepancies between the manually calculated abundance values and the pathway_errorbar chart.
Code: This is the code used to generate the pathway_errorbar chart.
Data: This is the Metacyc abundance data used in the analysis.
I would be immensely grateful for your assistance in resolving this issue. Your expertise and guidance would be invaluable in helping me obtain accurate and meaningful results from my ggpicrust2 analysis.

Thank you again for your time and consideration.

Sincerely,

Removing rowname

MC_abundance_tb <- as_tibble(MC_abundance_df)
column_to_rownames(MC_abundance_tb, var = "pathway")

Removing description column

MC_abundance_tb <- MC_abundance_tb %>% select(-last_col())

daa analysis using Maaslin2

metacyc_daa_results_df <- pathway_daa(abundance = MC_abundance_tb %>% column_to_rownames("pathway"), metadata = metadata, group = "Hostphylum", daa_method = "Maaslin2")

Finding significant difference (2x) between two groups (diatom, dinoflagellate)

average_expression_MC <- apply(MC_abundance_tb[-1], 1, function(x) {
diatom_mean <- mean(x[diatom_indexes], na.rm = TRUE)
dinoflagellate_mean <- mean(x[dinoflagellate_indexes], na.rm = TRUE)
c(diatom_mean = diatom_mean, dinoflagellate_mean = dinoflagellate_mean)
})
average_expression_matrix_MC <- matrix(average_expression_MC, nrow = 2, ncol = length(MC_abundance_tb$pathway), dimnames = list(c("diatom_mean", "dinoflagellate_mean"), MC_abundance_tb$pathway))
significant_functions_MC <- colnames(average_expression_matrix_MC)[which(average_expression_matrix_MC["diatom_mean", ] > 2 * average_expression_matrix_MC["dinoflagellate_mean", ] | average_expression_matrix_MC["dinoflagellate_mean", ] > 2 * average_expression_matrix_MC["diatom_mean", ])]

change for fitting with pathway ID in abundance data

significant_functions_MC_mod <- str_replace_all(significant_functions_MC, "-", ".")

filtering in daa results with significant difference results

filtered_df_MC <- metacyc_daa_results_df[metacyc_daa_results_df$feature %in% significant_functions_MC_mod,]

combine with annotation information

MC_daa_annotated_df <- left_join(filtered_df_MC, MC_abundance_df[, c("pathway", "Description")], by = c("feature" = "pathway"))

p_adjust value top30 only

MC_daa_annotated_df_top30 <- MC_daa_annotated_df %>% arrange(p_adjust) %>% slice_head(n = 30)

Making pathway_errorbar plot

p_MC_errorbar <- pathway_errorbar(abundance = MC_abundance_df_exep_desc %>% column_to_rownames("pathway"),
daa_results_df = MC_daa_annotated_df_top30,
Group = metadata$Hostphylum,
ko_to_kegg = FALSE,
p_values_threshold = 0.05,
order = "group",
select = NULL,
p_value_bar = TRUE,
colors = NULL,
x_lab = "Description")
print(p_MC_errorbar)

metadata.csv
MetaCyc_pathway_abundance.csv
MC abundance 박스플롯.xlsx