---
title: "Visualization Guide for METAFLUX"
author: "Zaoqu Liu"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Visualization Guide for METAFLUX}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)
```

## Overview

This vignette demonstrates visualization of METAFLUX results with real examples.

## Load Package and Data

```{r load-data}
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")
```

## Step 1: Calculate MRAS

```{r calc-mras}
# Calculate Metabolic Reaction Activity Scores
mras <- calculate_reaction_score(bulk_test_example)
cat("MRAS matrix:", nrow(mras), "reactions x", ncol(mras), "samples\n")
```

## Step 2: Compute Fluxes

```{r compute-flux}
# Compute metabolic fluxes
flux <- compute_flux(mras = mras, medium = human_blood)
cat("Flux matrix:", nrow(flux), "reactions x", ncol(flux), "samples\n")
```

## Visualization 1: Flux Distribution

```{r flux-distribution, fig.cap="Distribution of metabolic fluxes across all 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")
  )
```

## Visualization 2: Key Metabolite Exchange

```{r key-metabolites, fig.cap="Flux values for key metabolic reactions"}
# 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)
  )
```

## Visualization 3: Top Variable Reactions

```{r top-variable, fig.cap="Heatmap of top 30 most variable metabolic reactions"}
# 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)
  )
```

## Visualization 4: Pathway Activity

```{r pathway-activity, fig.cap="Mean absolute flux by metabolic pathway (top 15)"}
# 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)
  )
```

## Visualization 5: Sample Correlation

```{r sample-correlation, fig.cap="Correlation matrix of metabolic flux profiles"}
# 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)
  )
```

## Visualization 6: Exchange Reaction Analysis

```{r exchange-reactions, fig.cap="Distribution of uptake vs secretion reactions"}
# 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")
print(type_counts)

# 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")
  )
```

## Summary Statistics

```{r summary-stats}
cat("========================================\n")
cat("METAFLUX Analysis Summary\n")
cat("========================================\n\n")

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

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

cat("\nKey Metabolites (mean flux):\n")
for (met in names(key_rxns)) {
  val <- mean(flux[key_rxns[met], ])
  cat(sprintf("  - %s: %.4f\n", met, val))
}
```

## Session Information

```{r session-info}
sessionInfo()
```