Single-Cell RNA-seq Metabolic Flux Analysis

Overview

This vignette demonstrates metabolic flux analysis for single-cell RNA-seq data using METAFLUX’s community modeling approach.

Biological Context

Tumor Microenvironment Metabolism

Single-cell metabolic analysis reveals:

Metabolic heterogeneity within cell populations
Metabolic cooperation between cell types
Competition for shared nutrients
Metabolic crosstalk via secreted metabolites

Community Modeling Concept

┌─────────────────────────────────────────────────────────────────┐
│                   Tumor Microenvironment                        │
│  ┌──────────┐     ┌──────────┐     ┌──────────┐                │
│  │  Tumor   │ ←→  │  T cell  │ ←→  │  Macro-  │                │
│  │  Cells   │     │          │     │  phage   │                │
│  └────┬─────┘     └────┬─────┘     └────┬─────┘                │
│       │                │                │                       │
│       └────────────────┼────────────────┘                       │
│                        │                                        │
│                   ┌────┴────┐                                   │
│                   │ Shared  │                                   │
│                   │ Medium  │                                   │
│                   └─────────┘                                   │
└─────────────────────────────────────────────────────────────────┘

Data Preparation

Input Requirements

Component	Format	Description
Seurat object	v4/v5	Normalized counts in RNA assay
Cell annotations	Metadata column	Cell type labels
Medium profile	Data frame	Available nutrients

Load Example Data

library(METAFLUX)
library(Seurat)

# Load example single-cell data
data("sc_test_example")
data("human_blood")

# Inspect Seurat object
sc_test_example
table(sc_test_example$Cell_type)

Step 1: Bootstrap Aggregation

Single-cell data is sparse. Bootstrap aggregation reduces noise while preserving biological variation.

# Calculate bootstrap-averaged expression
avg_expr <- calculate_avg_exp(
  myseurat = sc_test_example,
  myident = "Cell_type",        # Column with cell type labels
  n_bootstrap = 100,            # Number of bootstrap iterations
  seed = 42                     # For reproducibility
)

# Check output dimensions
dim(avg_expr)
# Rows: genes
# Columns: cell_type_1, cell_type_2, ..., cell_type_1, ... (repeated n_bootstrap times)

Bootstrap Algorithm

For each bootstrap iteration i:
  1. Resample cells within each cell type (with replacement)
  2. Compute average expression per cell type
  3. Store result

Output: Matrix of [genes × (n_celltypes × n_bootstrap)]

Step 2: Calculate MRAS

# Calculate reaction activity scores
mras <- calculate_reaction_score(avg_expr)

# Verify dimensions
dim(mras)  # 13082 × (n_celltypes × n_bootstrap)

Step 3: Define Cell Type Fractions

Cell type proportions weight the community model:

# Calculate fractions from data
cell_counts <- table(sc_test_example$Cell_type)
fractions <- as.numeric(cell_counts / sum(cell_counts))
names(fractions) <- names(cell_counts)

# Verify: must sum to 1
sum(fractions)  # Should be 1

print(fractions)

Step 4: Community Flux Calculation

# Number of cell types
num_cell <- length(unique(sc_test_example$Cell_type))

# Compute community metabolic fluxes
flux <- compute_sc_flux(
  num_cell = num_cell,
  fraction = fractions,
  fluxscore = mras,
  medium = human_blood,
  parallel = TRUE,              # Enable parallel computing
  n_cores = 4,                  # Number of CPU cores
  use_cpp = TRUE                # Use C++ acceleration
)

# Check output
dim(flux)
head(rownames(flux), 20)

Output Structure

The flux matrix contains:

celltype 1 HMR_0001          # Cell type 1 reactions
celltype 1 HMR_0002
...
celltype 2 HMR_0001          # Cell type 2 reactions
celltype 2 HMR_0002
...
external_medium HMR_XXXX     # External exchange reactions

Step 5: Analyze Results

Extract Cell Type-Specific Fluxes

# Get reaction names
reaction_names <- rownames(flux)

# Extract fluxes for specific cell type
get_celltype_flux <- function(flux_matrix, celltype_num) {
  pattern <- paste0("^celltype ", celltype_num, " ")
  idx <- grep(pattern, rownames(flux_matrix))
  result <- flux_matrix[idx, , drop = FALSE]
  rownames(result) <- gsub(pattern, "", rownames(result))
  return(result)
}

# Extract for each cell type
tumor_flux <- get_celltype_flux(flux, 1)
tcell_flux <- get_celltype_flux(flux, 2)

Compare Metabolic Phenotypes

# Compare glucose uptake between cell types
glucose_rxn <- "HMR_9034"

tumor_glucose <- mean(tumor_flux[glucose_rxn, ])
tcell_glucose <- mean(tcell_flux[glucose_rxn, ])

cat(sprintf("Tumor glucose uptake: %.4f\n", tumor_glucose))
cat(sprintf("T cell glucose uptake: %.4f\n", tcell_glucose))

External Medium Exchange

# Extract external exchange reactions
external_idx <- grep("^external_medium", rownames(flux))
external_flux <- flux[external_idx, , drop = FALSE]

# Net nutrient consumption/production
nutrient_balance <- rowMeans(external_flux)
head(sort(nutrient_balance), 10)  # Top uptakes
tail(sort(nutrient_balance), 10)  # Top secretions

Advanced Analysis

Metabolic Competition

Identify metabolites competed for by different cell types:

# Compare lactate exchange across cell types
lactate_rxn <- "HMR_9135"

lactate_by_celltype <- sapply(1:num_cell, function(ct) {
  ct_flux <- get_celltype_flux(flux, ct)
  mean(ct_flux[lactate_rxn, ])
})

names(lactate_by_celltype) <- names(cell_counts)
print(lactate_by_celltype)

Metabolic Cooperation

Identify metabolites exchanged between cell types:

# One cell type produces, another consumes
producer <- which(lactate_by_celltype > 0)
consumer <- which(lactate_by_celltype < 0)

if (length(producer) > 0 && length(consumer) > 0) {
  cat(sprintf("%s produces lactate\n", names(lactate_by_celltype)[producer]))
  cat(sprintf("%s consumes lactate\n", names(lactate_by_celltype)[consumer]))
}

Performance Tips

Parallel Computing

# Check available cores
parallel::detectCores()

# Recommended: use physical cores - 1
n_cores <- max(1, parallel::detectCores(logical = FALSE) - 1)

# Limit to 8 to prevent memory issues
n_cores <- min(n_cores, 8)

Memory Management

# For large datasets, reduce bootstrap iterations
avg_expr <- calculate_avg_exp(
  myseurat = sc_test_example,
  myident = "Cell_type",
  n_bootstrap = 50,    # Reduced from 100
  seed = 42
)

# Or subsample cells first
sc_subset <- subset(sc_test_example, downsample = 5000)

Biological Interpretation

Warburg Effect

# High glycolysis, low oxidative phosphorylation
glycolysis_rxns <- c("HMR_4363", "HMR_4365", "HMR_4367")  # Example
oxphos_rxns <- c("HMR_6916", "HMR_6921")  # Example

tumor_flux <- get_celltype_flux(flux, 1)

glycolysis_mean <- mean(abs(tumor_flux[glycolysis_rxns, ]))
oxphos_mean <- mean(abs(tumor_flux[oxphos_rxns, ]))

warburg_ratio <- glycolysis_mean / oxphos_mean
cat(sprintf("Warburg ratio: %.2f\n", warburg_ratio))

Immunometabolism

# T cell activation markers
glutamine_rxn <- "HMR_9063"  # Glutamine uptake
fa_synthesis <- grep("fatty acid synthesis", rownames(tumor_flux), ignore.case = TRUE)

tcell_flux <- get_celltype_flux(flux, 2)
tcell_glutamine <- mean(abs(tcell_flux[glutamine_rxn, ]))
tcell_fa <- mean(abs(tcell_flux[fa_synthesis, ]))

cat(sprintf("T cell glutamine uptake: %.4f\n", tcell_glutamine))
cat(sprintf("T cell fatty acid synthesis: %.4f\n", tcell_fa))

Output Files

Save results for downstream analysis:

# Save flux matrix
write.csv(flux, "metabolic_flux_singlecell.csv")

# Save cell type-specific fluxes
for (ct in 1:num_cell) {
  ct_flux <- get_celltype_flux(flux, ct)
  filename <- sprintf("flux_celltype_%d.csv", ct)
  write.csv(ct_flux, filename)
}

Next Steps

Visualization: See Visualization Guide
Pathway analysis: Aggregate fluxes by metabolic pathway
Statistical testing: Compare conditions
Integration: Combine with other omics data

Author

Zaoqu Liu - Package maintainer - GitHub: Zaoqu-Liu - ORCID: 0000-0002-0452-742X

Session Information

sessionInfo()