Visualization Guide for METAFLUX

library(METAFLUX)
library(ggplot2)

# Load example data
data("bulk_test_example")
data("human_blood")
data("human_gem")

# Quick preview
cat("Expression matrix:", nrow(bulk_test_example), "genes x", ncol(bulk_test_example), "samples\n")
#> Expression matrix: 58581 genes x 5 samples

# Calculate Metabolic Reaction Activity Scores
mras <- calculate_reaction_score(bulk_test_example)
cat("MRAS matrix:", nrow(mras), "reactions x", ncol(mras), "samples\n")
#> MRAS matrix: 13082 reactions x 5 samples

# Compute metabolic fluxes
flux <- compute_flux(mras = mras, medium = human_blood)
#>   |                                                                              |                                                                      |   0%  |                                                                              |==============                                                        |  20%
#>   |                                                                              |============================                                          |  40%
#>   |                                                                              |==========================================                            |  60%
#>   |                                                                              |========================================================              |  80%
#>   |                                                                              |======================================================================| 100%
cat("Flux matrix:", nrow(flux), "reactions x", ncol(flux), "samples\n")
#> Flux matrix: 13082 reactions x 5 samples

# Prepare data
flux_values <- as.vector(as.matrix(flux))
flux_df <- data.frame(flux = flux_values[abs(flux_values) < 0.5])

# Plot
ggplot(flux_df, aes(x = flux)) +
  geom_histogram(bins = 100, fill = "#3498db", color = "white", alpha = 0.8) +
  labs(
    title = "Distribution of Metabolic Fluxes",
    subtitle = "METAFLUX analysis of bulk RNA-seq data",
    x = "Flux Value",
    y = "Count"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray50")
  )

# Define key reactions
key_rxns <- c(
  "Glucose" = "HMR_9034",
  "Lactate" = "HMR_9135",
  "Glutamine" = "HMR_9063",
  "Pyruvate" = "HMR_9133"
)

# Extract flux values
key_flux <- flux[key_rxns, , drop = FALSE]
rownames(key_flux) <- names(key_rxns)

# Convert to long format
key_df <- data.frame(
  Metabolite = rep(rownames(key_flux), ncol(key_flux)),
  Sample = rep(colnames(key_flux), each = nrow(key_flux)),
  Flux = as.vector(as.matrix(key_flux))
)

# Plot
ggplot(key_df, aes(x = Metabolite, y = Flux, fill = Sample)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(
    title = "Key Metabolite Exchange",
    subtitle = "Negative = uptake, Positive = secretion",
    x = "",
    y = "Flux"
  ) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray50"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 11)
  )

# Calculate variance
flux_var <- apply(flux, 1, var)
top_30 <- names(sort(flux_var, decreasing = TRUE)[1:30])

# Prepare matrix
flux_top <- as.matrix(flux[top_30, ])

# Scale for visualization
flux_scaled <- t(scale(t(flux_top)))

# Convert to data frame
heatmap_df <- data.frame(
  Reaction = rep(rownames(flux_scaled), ncol(flux_scaled)),
  Sample = rep(colnames(flux_scaled), each = nrow(flux_scaled)),
  Flux = as.vector(flux_scaled)
)

# Order reactions by mean flux
rxn_order <- rownames(flux_scaled)[order(rowMeans(flux_scaled))]
heatmap_df$Reaction <- factor(heatmap_df$Reaction, levels = rxn_order)

# Plot heatmap
ggplot(heatmap_df, aes(x = Sample, y = Reaction, fill = Flux)) +
  geom_tile() +
  scale_fill_gradient2(
    low = "#2166ac", mid = "white", high = "#b2182b",
    midpoint = 0, name = "Z-score"
  ) +
  labs(
    title = "Top 30 Variable Metabolic Reactions",
    x = "Sample",
    y = "Reaction"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text.y = element_text(size = 6),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Get pathway information
pathway_map <- human_gem$SUBSYSTEM
names(pathway_map) <- human_gem$ID

# Calculate mean flux per pathway
flux_matrix <- as.matrix(flux)
pathways <- unique(pathway_map)

pathway_flux <- sapply(pathways, function(pw) {
  rxns <- names(pathway_map)[pathway_map == pw]
  rxns <- intersect(rxns, rownames(flux_matrix))
  if (length(rxns) > 0) {
    mean(abs(flux_matrix[rxns, ]), na.rm = TRUE)
  } else {
    NA
  }
})

# Remove NA and get top 15
pathway_flux <- sort(pathway_flux[!is.na(pathway_flux)], decreasing = TRUE)
top_pathways <- head(pathway_flux, 15)

# Create data frame
pathway_df <- data.frame(
  Pathway = factor(names(top_pathways), levels = rev(names(top_pathways))),
  Activity = top_pathways
)

# Plot
ggplot(pathway_df, aes(x = Activity, y = Pathway)) +
  geom_bar(stat = "identity", fill = "#2c3e50", width = 0.7) +
  labs(
    title = "Top 15 Active Metabolic Pathways",
    subtitle = "Based on mean absolute flux",
    x = "Mean Absolute Flux",
    y = ""
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray50"),
    axis.text.y = element_text(size = 9)
  )

# Calculate correlation
cor_matrix <- cor(as.matrix(flux), use = "pairwise.complete.obs")

# Convert to long format
cor_df <- data.frame(
  Sample1 = rep(rownames(cor_matrix), ncol(cor_matrix)),
  Sample2 = rep(colnames(cor_matrix), each = nrow(cor_matrix)),
  Correlation = as.vector(cor_matrix)
)

# Plot
ggplot(cor_df, aes(x = Sample1, y = Sample2, fill = Correlation)) +
  geom_tile() +
  geom_text(aes(label = round(Correlation, 2)), size = 4) +
  scale_fill_gradient2(
    low = "#2166ac", mid = "white", high = "#b2182b",
    midpoint = 0.5, limits = c(0, 1)
  ) +
  labs(
    title = "Sample Correlation Matrix",
    subtitle = "Based on metabolic flux profiles",
    x = "", y = ""
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray50"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Identify exchange reactions
exchange_idx <- grep("Exchange", human_gem$SUBSYSTEM)
exchange_rxns <- human_gem$ID[exchange_idx]
exchange_rxns <- intersect(exchange_rxns, rownames(flux))

# Get mean flux
exchange_flux <- rowMeans(flux[exchange_rxns, ])

# Classify
exchange_df <- data.frame(
  Reaction = names(exchange_flux),
  Flux = exchange_flux,
  Type = ifelse(exchange_flux < -0.001, "Uptake",
                ifelse(exchange_flux > 0.001, "Secretion", "Inactive"))
)

# Summary counts
type_counts <- table(exchange_df$Type)
cat("\nExchange reaction summary:\n")
#> 
#> Exchange reaction summary:
print(type_counts)
#> 
#>  Inactive Secretion    Uptake 
#>      1642         3         3

# Plot distribution
ggplot(exchange_df[exchange_df$Type != "Inactive", ], 
       aes(x = Flux, fill = Type)) +
  geom_histogram(bins = 50, alpha = 0.7, position = "identity") +
  scale_fill_manual(values = c("Uptake" = "#e74c3c", "Secretion" = "#27ae60")) +
  labs(
    title = "Exchange Reaction Distribution",
    subtitle = sprintf("Uptake: %d, Secretion: %d, Inactive: %d",
                      type_counts["Uptake"], type_counts["Secretion"], 
                      type_counts["Inactive"]),
    x = "Mean Flux",
    y = "Count"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray50")
  )

cat("========================================\n")
#> ========================================
cat("METAFLUX Analysis Summary\n")
#> METAFLUX Analysis Summary
cat("========================================\n\n")
#> ========================================

cat("Input Data:\n")
#> Input Data:
cat(sprintf("  - Genes: %d\n", nrow(bulk_test_example)))
#>   - Genes: 58581
cat(sprintf("  - Samples: %d\n", ncol(bulk_test_example)))
#>   - Samples: 5

cat("\nOutput:\n")
#> 
#> Output:
cat(sprintf("  - Reactions analyzed: %d\n", nrow(flux)))
#>   - Reactions analyzed: 13082
cat(sprintf("  - Flux range: [%.4f, %.4f]\n", min(flux), max(flux)))
#>   - Flux range: [-0.0068, 0.0075]

cat("\nKey Metabolites (mean flux):\n")
#> 
#> Key Metabolites (mean flux):
for (met in names(key_rxns)) {
  val <- mean(flux[key_rxns[met], ])
  cat(sprintf("  - %s: %.4f\n", met, val))
}
#>   - Glucose: -0.0003
#>   - Lactate: 0.0003
#>   - Glutamine: -0.0004
#>   - Pyruvate: 0.0004

sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3  METAFLUX_2.2.0 rmarkdown_2.31
#> 
#> loaded via a namespace (and not attached):
#>   [1] deldir_2.0-4           pbapply_1.7-4          gridExtra_2.3         
#>   [4] rlang_1.2.0            magrittr_2.0.5         RcppAnnoy_0.0.23      
#>   [7] otel_0.2.0             spatstat.geom_3.8-1    matrixStats_1.5.0     
#>  [10] ggridges_0.5.7         compiler_4.6.0         png_0.1-9             
#>  [13] vctrs_0.7.3            reshape2_1.4.5         stringr_1.6.0         
#>  [16] pkgconfig_2.0.3        fastmap_1.2.0          labeling_0.4.3        
#>  [19] promises_1.5.0         purrr_1.2.2            xfun_0.59             
#>  [22] cachem_1.1.0           jsonlite_2.0.0         goftest_1.2-3         
#>  [25] later_1.4.8            spatstat.utils_3.2-3   irlba_2.3.7           
#>  [28] parallel_4.6.0         cluster_2.1.8.2        R6_2.6.1              
#>  [31] ica_1.0-3              spatstat.data_3.1-9    bslib_0.11.0          
#>  [34] stringi_1.8.7          RColorBrewer_1.1-3     reticulate_1.46.0     
#>  [37] spatstat.univar_3.2-0  parallelly_1.47.0      lmtest_0.9-40         
#>  [40] jquerylib_0.1.4        scattermore_1.2        iterators_1.0.14      
#>  [43] Rcpp_1.1.1-1.1         knitr_1.51             tensor_1.5.1          
#>  [46] future.apply_1.20.2    zoo_1.8-15             sctransform_0.4.3     
#>  [49] httpuv_1.6.17          Matrix_1.7-5           splines_4.6.0         
#>  [52] igraph_2.3.2           tidyselect_1.2.1       abind_1.4-8           
#>  [55] yaml_2.3.12            doParallel_1.0.17      spatstat.random_3.5-0 
#>  [58] codetools_0.2-20       miniUI_0.1.2           spatstat.explore_3.8-1
#>  [61] listenv_1.0.0          lattice_0.22-9         tibble_3.3.1          
#>  [64] plyr_1.8.9             withr_3.0.3            shiny_1.14.0          
#>  [67] S7_0.2.2               ROCR_1.0-12            evaluate_1.0.5        
#>  [70] Rtsne_0.17             future_1.70.0          fastDummies_1.7.6     
#>  [73] survival_3.8-6         polyclip_1.10-7        fitdistrplus_1.2-6    
#>  [76] osqp_1.0.0             pillar_1.11.1          Seurat_5.5.0          
#>  [79] KernSmooth_2.23-26     foreach_1.5.2          plotly_4.12.0         
#>  [82] generics_0.1.4         RcppHNSW_0.7.0         sp_2.2-1              
#>  [85] scales_1.4.0           globals_0.19.1         xtable_1.8-8          
#>  [88] glue_1.8.1             lazyeval_0.2.3         maketools_1.3.2       
#>  [91] tools_4.6.0            sys_3.4.3              data.table_1.18.4     
#>  [94] RSpectra_0.16-2        RANN_2.6.2             buildtools_1.0.0      
#>  [97] dotCall64_1.2          cowplot_1.2.0          grid_4.6.0            
#> [100] tidyr_1.3.2            nlme_3.1-169           patchwork_1.3.2       
#> [103] cli_3.6.6              spatstat.sparse_3.2-0  spam_2.11-4           
#> [106] viridisLite_0.4.3      dplyr_1.2.1            uwot_0.2.4            
#> [109] gtable_0.3.6           sass_0.4.10            digest_0.6.39         
#> [112] progressr_0.19.0       ggrepel_0.9.8          htmlwidgets_1.6.4     
#> [115] SeuratObject_5.4.0     farver_2.1.2           htmltools_0.5.9       
#> [118] lifecycle_1.0.5        httr_1.4.8             mime_0.13             
#> [121] MASS_7.3-65