Bulk RNA-seq Metabolic Flux Analysis
Overview
This vignette demonstrates the complete workflow for metabolic flux analysis from bulk RNA-seq data using METAFLUX.
Load Package and Data
library(METAFLUX)
library(ggplot2)
# Load example bulk RNA-seq data
data("bulk_test_example")
data("human_blood")
data("human_gem")
# Inspect data
cat("Gene expression matrix:\n")
#> Gene expression matrix:
cat(" Dimensions:", dim(bulk_test_example), "\n")
#> Dimensions: 58581 5
cat(" Samples:", colnames(bulk_test_example), "\n")
#> Samples: Sample1 Sample2 Sample3 Sample4 Sample5
cat(" Example genes:", head(rownames(bulk_test_example), 5), "\n")
#> Example genes: TSPAN6 TNMD DPM1 SCYL3 C1orf112Step 1: Calculate MRAS
Metabolic Reaction Activity Scores (MRAS) integrate gene expression with metabolic network topology.
# Calculate MRAS
mras <- calculate_reaction_score(bulk_test_example)
# Check output
cat("\nMRAS matrix:\n")
#>
#> MRAS matrix:
cat(" Dimensions:", dim(mras), "\n")
#> Dimensions: 13082 5
cat(" Score range:", range(as.matrix(mras)), "\n")
#> Score range: 0 1MRAS Distribution
mras_values <- as.vector(as.matrix(mras))
mras_df <- data.frame(MRAS = mras_values[mras_values > 0])
ggplot(mras_df, aes(x = MRAS)) +
geom_histogram(bins = 50, fill = "#3498db", color = "white", alpha = 0.8) +
labs(
title = "Distribution of Metabolic Reaction Activity Scores",
x = "MRAS Value",
y = "Count"
) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
Distribution of MRAS values
Step 2: Compute Metabolic Fluxes
# Run flux balance analysis
flux <- compute_flux(mras = mras, medium = human_blood)
#> | | | 0% | |============== | 20%
#> | |============================ | 40%
#> | |========================================== | 60%
#> | |======================================================== | 80%
#> | |======================================================================| 100%
# Check results
cat("\nFlux matrix:\n")
#>
#> Flux matrix:
cat(" Dimensions:", dim(flux), "\n")
#> Dimensions: 13082 5
cat(" Flux range:", range(flux), "\n")
#> Flux range: -0.006807341 0.007501395Flux Interpretation
# Classify reactions by flux direction
flux_mean <- rowMeans(flux)
flux_class <- data.frame(
Direction = c("Production (+)", "Consumption (-)", "Inactive (≈0)"),
Count = c(
sum(flux_mean > 0.001),
sum(flux_mean < -0.001),
sum(abs(flux_mean) <= 0.001)
)
)
ggplot(flux_class, aes(x = Direction, y = Count, fill = Direction)) +
geom_bar(stat = "identity", width = 0.6) +
geom_text(aes(label = Count), vjust = -0.5, size = 5) +
scale_fill_manual(values = c("#27ae60", "#e74c3c", "#95a5a6")) +
labs(
title = "Reaction Classification by Flux Direction",
x = "",
y = "Number of Reactions"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "none"
) +
ylim(0, max(flux_class$Count) * 1.1)
Interpretation of flux values
Step 3: Key Metabolic Indicators
Central Carbon Metabolism
# Define key reactions
central_carbon <- c(
"Glucose uptake" = "HMR_9034",
"Lactate secretion" = "HMR_9135",
"Pyruvate transport" = "HMR_9133",
"Glutamine uptake" = "HMR_9063"
)
# Extract and plot
cc_flux <- flux[central_carbon, , drop = FALSE]
rownames(cc_flux) <- names(central_carbon)
cc_df <- data.frame(
Reaction = rep(rownames(cc_flux), ncol(cc_flux)),
Sample = rep(colnames(cc_flux), each = nrow(cc_flux)),
Flux = as.vector(as.matrix(cc_flux))
)
ggplot(cc_df, aes(x = Reaction, y = Flux, fill = Sample)) +
geom_bar(stat = "identity", position = "dodge") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_fill_brewer(palette = "Set2") +
labs(
title = "Central Carbon Metabolism",
subtitle = "Key reaction fluxes across samples",
x = "",
y = "Flux"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 30, hjust = 1)
)
Central carbon metabolism flux
Biomass Production
biomass_flux <- flux["biomass_human", ]
biomass_df <- data.frame(
Sample = names(biomass_flux),
Flux = as.numeric(biomass_flux)
)
ggplot(biomass_df, aes(x = Sample, y = Flux, fill = Sample)) +
geom_bar(stat = "identity", width = 0.6) +
scale_fill_brewer(palette = "Set2") +
labs(
title = "Biomass Production Rate",
subtitle = "Proxy for cellular growth rate",
x = "Sample",
y = "Biomass Flux"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "none"
)
Biomass production flux (growth rate proxy)
Step 4: Pathway Analysis
# Map reactions to pathways
pathway_map <- human_gem$SUBSYSTEM
names(pathway_map) <- human_gem$ID
# Calculate pathway activity
pathways <- unique(pathway_map)
pathway_activity <- sapply(pathways, function(pw) {
rxns <- names(pathway_map)[pathway_map == pw]
rxns <- intersect(rxns, rownames(flux))
if (length(rxns) > 0) mean(abs(flux[rxns, ])) else NA
})
# Top 12 pathways
top_pw <- sort(pathway_activity[!is.na(pathway_activity)], decreasing = TRUE)[1:12]
pw_df <- data.frame(
Pathway = factor(names(top_pw), levels = rev(names(top_pw))),
Activity = top_pw
)
ggplot(pw_df, aes(x = Activity, y = Pathway)) +
geom_bar(stat = "identity", fill = "#2c3e50") +
labs(
title = "Top 12 Active Metabolic Pathways",
x = "Mean Absolute Flux",
y = ""
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.y = element_text(size = 9)
)
Top metabolic pathways by activity
Step 5: Quality Control
Steady-State Verification
# Load stoichiometric matrix
Hgem <- METAFLUX:::Hgem
# Check steady-state constraint: S*v = 0
sv_violations <- sapply(1:ncol(flux), function(i) {
sv <- Hgem$S %*% flux[, i]
max(abs(sv))
})
cat("Steady-state constraint check:\n")
#> Steady-state constraint check:
cat(" Max violations per sample:", round(sv_violations, 6), "\n")
#> Max violations per sample: 0 0 0 0 0
cat(" All < 0.001:", all(sv_violations < 0.001), "\n")
#> All < 0.001: TRUEGene Coverage
gene_num <- METAFLUX:::gene_num
input_genes <- rownames(bulk_test_example)
metabolic_genes <- rownames(gene_num)
coverage <- sum(input_genes %in% metabolic_genes)
total <- length(metabolic_genes)
cov_df <- data.frame(
Category = c("Detected", "Missing"),
Count = c(coverage, total - coverage)
)
ggplot(cov_df, aes(x = "", y = Count, fill = Category)) +
geom_bar(stat = "identity", width = 1) +
coord_polar("y") +
scale_fill_manual(values = c("Detected" = "#27ae60", "Missing" = "#e74c3c")) +
labs(
title = "Metabolic Gene Coverage",
subtitle = sprintf("%.1f%% (%d/%d genes)", coverage/total*100, coverage, total)
) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5)
)
Metabolic gene coverage
Summary
cat("================================================\n")
#> ================================================
cat("METAFLUX Bulk RNA-seq Analysis Complete\n")
#> METAFLUX Bulk RNA-seq Analysis Complete
cat("================================================\n\n")
#> ================================================
cat("Input:\n")
#> Input:
cat(sprintf(" Genes: %d\n", nrow(bulk_test_example)))
#> Genes: 58581
cat(sprintf(" Samples: %d\n", ncol(bulk_test_example)))
#> Samples: 5
cat(sprintf(" Metabolic gene coverage: %.1f%%\n", coverage/total*100))
#> Metabolic gene coverage: 96.2%
cat("\nOutput:\n")
#>
#> Output:
cat(sprintf(" Reactions: %d\n", nrow(flux)))
#> Reactions: 13082
cat(sprintf(" Flux range: [%.4f, %.4f]\n", min(flux), max(flux)))
#> Flux range: [-0.0068, 0.0075]
cat("\nKey findings:\n")
#>
#> Key findings:
cat(sprintf(" Active reactions: %d\n", sum(abs(rowMeans(flux)) > 0.001)))
#> Active reactions: 39
cat(sprintf(" Glucose uptake (mean): %.4f\n", mean(flux["HMR_9034", ])))
#> Glucose uptake (mean): -0.0003
cat(sprintf(" Lactate secretion (mean): %.4f\n", mean(flux["HMR_9135", ])))
#> Lactate secretion (mean): 0.0003Session Information
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