Computational frameworks for metabolic modeling from snRNA-seq in GBM (De Temmerman et al., 2026)

Paper

• De Temmerman M, Vandemoortele B, Vermeirssen V. “Comparison of different computational frameworks for metabolic modeling from single-cell transcriptomics data in glioblastoma.” bioRxiv (2026). DOI: https://doi.org/10.64898/2026.02.28.706749

Background

The problem

Metabolic reprogramming is a hallmark of glioblastoma (GBM), but how individual cell populations — malignant states, immune cells, stromal populations — each contribute to the metabolic heterogeneity within a tumor remains poorly defined. Bulk metabolomics cannot resolve these contributions, and single-cell metabolomics is still technically limited by low throughput, limited metabolite coverage, and the inability to amplify metabolites. Computational inference from single-cell transcriptomics offers a powerful alternative, but no systematic comparison of different inference methods applied to the same GBM dataset existed.

Study design

The authors applied three complementary computational methods to snRNA-seq data from 6 GBM patient samples (2 per Wang multiomic subtype: nmf1/PN-like, nmf2/MES-like, nmf3/CL-like), all mapped against the Human1 genome-scale metabolic model (GEM) as a unified annotation reference:

Core finding

Across all three methods, tumor-associated macrophages (TAMs) emerge as the metabolically dominant cell population in the GBM TME. TAMs in MES-like tumors show coordinated lipid metabolism driven by five recurrent transcription factors, consistent nucleotide biosynthesis flux, and glutamate-to-glutamine conversion that potentially supports malignant cells.

Methodology details

1. Metabolic pathway activity scoring

Pathway-level analysis adapted from Locasale (Xiao et al., 2019). Metabolic pathway gene sets are derived from Human1 GEM. Expression of enzyme-coding genes is aggregated per cell type, with an inverse-frequency weighting scheme so multi-pathway genes don’t dominate. Statistical significance is assessed by 5000-fold permutation testing against random gene sets of matched size.

Key results:

• Protein degradation, drug metabolism, protein assembly, and even-chain fatty acid desaturation are consistently top-ranking pathways across patients
• TAMs and monocytes drive protein degradation/drug metabolism; oligodendrocytes drive fatty acid desaturation
• MES-like samples show TAM-enriched fatty acid activation and metabolism
• Most malignant cell states show limited pathway activity, except MES-like cells in protein modification and ROS detoxification

2. SCENIC-based metabolic GRN inference

SCENIC (Aibar et al., 2017) reconstructs regulons (TFs + predicted target genes) and quantifies their activity per cell via AUCell. The authors ran SCENIC 10$\times$ per sample (GRNBoost2 + cisTarget + AUCell), keeping consensus edges present in $\ge 80\%$ of runs. Regulon target genes were then overlapped with Human1 metabolic pathway gene sets to build tripartite networks (TF $\rightarrow$ target gene $\rightarrow$ pathway). A connection specificity index (CSI) adapted from Bass et al. was used to identify regulon modules with similar metabolic pathway regulation profiles.

Key results:

• Most regulons (54.3%) are sample-specific, reinforcing the need for per-patient analysis
• A conserved core of 30 regulons spans both immune and malignant compartments across nearly all samples
• Five TAM-specific TFs (CEBPD, FOSL2, HIF1A, IRF5, PRDM1) converge on a shared lipid-regulatory program in MES-like samples — regulating sphingolipid, glycerophospholipid, and arachidonic acid metabolism
• E2F family members drive nucleotide biosynthesis in malignant cells

3. scFEA single-cell flux prediction

scFEA (Alghamdi et al., 2021) uses a graph neural network to predict cell-wise metabolic flux distributions across 168 reaction modules derived from a reduced GEM. The tool integrates single-cell expression with stoichiometric constraints from flux balance analysis. TAM-specific modules were identified when Cohen’s $d$ > 0.8 and $>50\%$ of TAMs showed positive flux.

Key results:

• TAMs form distinct clusters in flux-based UMAP space across all samples
• Purine and pyrimidine synthesis modules show the strongest and most consistent TAM-specific flux (AICAR $\rightarrow$ IMP; UMP $\rightarrow$ $\beta$-alanine; dTMP $\rightarrow$ succinyl-CoA)
• Glutamate $\rightarrow$ glutamine conversion is a recurrent TAM-specific flux across all samples (the “metabolic tug-of-war” with tumor cells)
• Fatty acid $\rightarrow$ acetyl-CoA flux is TAM-specific only in nmf3 (CL-like) samples — absent in nmf2 MES-like samples despite strong transcriptional signals there, revealing method-specific sensitivity differences

Cross-method comparison

The three methods span a hierarchy of increasing granularity:

Pathway activity (which pathways are upregulated?) → GRN inference (which TFs drive those pathways?) → Flux prediction (which reactions carry flux?)

Where all three converge — especially on TAM metabolic prominence — the findings are likely robust. Where they diverge (e.g., fatty acid metabolism in MES-like TAMs visible transcriptionally but not in flux), it reflects methodological differences in how each approach captures metabolic state from transcript data.

Pseudo data: input → method → output

Method 1: Metabolic pathway activity scoring

Input:

+ Human1 GEM pathway-gene mapping (e.g., Glycolysis → {Gene_A, Gene_B}; Lipid metabolism → {Gene_C, Gene_D})

Process: For each cell type, calculate weighted mean expression of pathway genes (weight = 1/number_of_pathways_per_gene). Compare against 5000 permuted random gene sets of same size to get empirical $p$-value.

Output:

Advantages:

• Simple, interpretable, and broadly applicable
• Directly maps to biological pathway knowledge
• Works well for identifying “which cell types are metabolically dominant in which pathways”
• Computationally lightweight

Limitations

Some caveats to note for this first method:

1. Cell-type-level resolution only (no single-cell granularity)
2. Ignores network structure and reaction connectivity
3. Cannot identify regulatory mechanisms behind pathway changes

Method 2: SCENIC GRN inference + metabolic overlap

Input:

+ TF binding motif databases (cisTarget) + Human1 pathway-gene mapping

Process:

1. GRNBoost2 infers co-expression modules (TF → candidate target genes)
2. cisTarget prunes to targets with matching TF binding motifs
3. AUCell scores regulon activity per cell
4. Overlap regulon target genes with Human1 metabolic pathway gene sets → tripartite network (TF $\rightarrow$ target gene $\rightarrow$ pathway)
5. Compute CSI between regulon pairs to find modules with similar metabolic regulation profiles

Output:

Regulon activity scores (per cell):

Regulon-pathway overlap matrix:

Advantages:

• Reveals the regulatory logic behind metabolic states (which TFs drive which pathways)
• Single-cell resolution for regulon activity
• Can identify conserved vs. patient-specific regulatory programs
• CSI analysis groups TFs with similar metabolic-regulatory strategies

Limitations

Some caveats to note for the GRN approach:

1. Regulon composition is sample-specific (hard to directly compare regulon content across patients)
2. Only examines transcriptional regulation — misses post-transcriptional, allosteric, and metabolite-level control
3. Does not account for pathway connectivity or stoichiometric constraints

Method 3: scFEA flux prediction

Input:

+ scFEA metabolic module definitions (168 modules) + stoichiometric constraint matrix (cmMat_c70_m168.csv)

Process: A graph neural network maps gene expression → predicted flux for each of the 168 reaction modules per cell, subject to metabolite balance constraints (flux balance analysis). TAM-specificity is assessed via Cohen’s $d$ (TAM vs. non-TAM) and proportion of TAMs with positive flux.

Output:

Advantages:

• Incorporates stoichiometric/network constraints — not just expression levels
• Sub-pathway, reaction-level granularity (168 modules)
• Single-cell flux profiles enable metabolic UMAP clustering
• Captures metabolic dependencies invisible to transcript-level analysis alone

Limitations

Some caveats to note for flux prediction:

1. Depends on a reduced/linearized GEM — cannot capture all metabolic reactions (e.g., fatty acid oxidation in MES-TAMs missed)
2. Computational cost is higher than pathway scoring
3. Module definitions may not align with standard pathway annotations (requires additional mapping)
4. Flux predictions are inferred, not measured — validation against metabolomics is essential

Extending to spatial context

Why spatial matters for metabolic inference

All three methods in this study operate on dissociated single-cell/nucleus data and therefore lose the spatial organization of the tumor. In GBM, spatial context is critical: metabolic states are shaped by gradients of oxygen, nutrients, and signaling molecules that radiate from blood vessels, necrotic cores, and the infiltrative margin. A TAM adjacent to the necrotic core likely has a very different metabolic program from a TAM at the leading edge, even if both are transcriptionally classified as “TAM” in a dissociated dataset. Without spatial information, these metabolic subtypes are averaged together or, at best, separated only when they differ enough in gene expression to form distinct clusters.

Spatial transcriptomics technologies — Visium, MERFISH, Slide-seq, and 10X Xenium — now provide the means to anchor transcriptional profiles to tissue coordinates. This opens the possibility of spatially resolved metabolic inference: applying the same three computational frameworks (pathway activity, GRN, flux prediction) while preserving the spatial relationships between cells and their microenvironment. The Human1 GEM annotation framework used in this study would transfer directly to spatial data, providing a unified metabolic reference regardless of the measurement platform.

Practical approaches

Spatial pathway activity scoring is the most immediately accessible extension. For spot-based platforms like Visium (each spot covers ~5–10 cells), one can compute weighted pathway scores per spot using the same Human1 gene sets and then overlay these onto the tissue image. This would reveal, for example, whether lipid metabolism is concentrated around vascular niches or whether glycolysis peaks in hypoxic regions surrounding necrosis. For single-molecule platforms (MERFISH, Xenium) that achieve true single-cell resolution, pathway scoring could directly identify spatial clusters of metabolically similar cells — but the limited gene panels (hundreds to a few thousand genes) would require careful selection of metabolic gene targets during panel design.

Spatial GRN inference is more challenging but tractable. Tools like SpatialDM and spatially informed versions of SCENIC (e.g., based on SCENIC+) can incorporate spatial proximity to identify spatially coherent regulon modules. One could ask: do the five TAM-specific TFs (CEBPD, FOSL2, HIF1A, IRF5, PRDM1) that converge on lipid metabolism co-activate in spatially localized TAM niches? If so, this would transform the finding from “TAMs regulate lipid metabolism” to “TAMs in region X regulate lipid metabolism in concert with neighboring MES-like cells.” Ligand–receptor analysis tools (CellChat, COMMOT for spatial data) could then connect the metabolic GRN findings to intercellular communication axes, testing whether the glutamate-glutamine “metabolic tug-of-war” detected by scFEA actually occurs between spatially adjacent TAMs and tumor cells.

Spatial flux modeling represents the frontier. The core challenge is that scFEA and related tools assume cell autonomy — each cell’s flux distribution is predicted independently. In reality, metabolite exchange between neighboring cells is central to tumor metabolism: TAMs export glutamine that tumor cells import; tumor cells secrete lactate that rewires macrophage polarization. A spatially informed flux model would need to impose inter-cellular metabolite balance constraints — effectively treating the tissue as a multi-compartment system where metabolites flow between adjacent spots or cells. This is conceptually related to community-level FBA methods (e.g., SteadyCom, OptCom), which model metabolic interactions between microbial species, but would need substantial adaptation for mammalian tissue geometry. Early steps could include computing pairwise metabolite exchange scores between spatially adjacent cell types and testing whether scFEA-predicted export fluxes in one cell type correlate with import fluxes in neighbors.

What spatial resolution would add to this study’s conclusions

The most impactful spatial extension of this study’s findings would be to test whether the TAM metabolic states identified computationally correspond to spatially organized niches. The MES-like subtype’s lipid-regulatory TAM program might localize to myelin-rich border zones where macrophages phagocytose debris, while nucleotide biosynthesis-active TAMs might cluster near proliferative tumor zones. If these spatial predictions hold, the three-method framework used here — activity → regulation → flux — could be applied as a spatial pipeline, moving from “what metabolic programs exist” to “where do they operate and with whom do they interact.” This would substantially strengthen the translational relevance of these computational approaches and directly inform spatially targeted therapeutic strategies.