scMetaTraj: A Metabolism‐Centered Framework for Trajectory Inference in Single‐Cell Data

Metabolic reprogramming is often continuous rather than discrete, yet most single-cell workflows either focus on transcriptome-defined lineages or report metabolism only as a secondary annotation. We developed scMetaTraj, a metabolism-centered R package that scores curated metabolic modules, embeds cells in a metabolic state space, and infers metabolic pseudotime (mPT) from graph distances within that space. scMetaTraj is intended to complement rather than replace classical trajectory tools. Methods such as Monocle 3 and Slingshot infer trajectories in transcriptomic space, whereas scMetabolism, compass, and scFEA quantify metabolism without reconstructing ordered metabolic trajectories. scMetaTraj, instead, performs trajectory inference directly on pathway- and module-level metabolic scores. A metabolism-first perspective is useful when cells with similar transcriptomic identities occupy distinct functional metabolic states or when metabolic transitions occur gradually across activation, stress adaptation, and remodeling. In this setting, distances in the inferred manifold reflect similarity in pathway composition rather than full-transcriptome similarity, enabling the detection of transition zones and intermediate metabolic configurations that may be obscured by transcriptomic clustering alone. An illustrative application is shown in Figure 1, which reconstructs a metabolic gradient within the macrophage lineage, summarizes metabolic cluster positions in a two-axis metabolic state space, compares ranked module activity patterns across clusters, and visualizes the distribution of metabolic subclusters along mPT. Metabolism-centered representation of state variation and subcluster structure. (A) UMAP view of the metabolic gradient within the macrophage lineage, colored by glycolytic axis score. (B) Two-axis metabolic state space summarizing metabolic clusters by glycolytic and mitochondrial axes, with an accompanying bar plot of lipid metabolic intensity across clusters. (C) Module cards showing ranked activity patterns of representative metabolic programs across metabolic clusters. (D) Ridgeline distributions of metabolic subclusters along metabolic pseudotime (mPT). In the PBMC reference dataset, scMetaTraj identified 10 metabolic subclusters distributed over a continuous manifold. Overlay of transcriptome-derived labels showed that the metabolic embedding did not simply replicate canonical PBMC identities; instead, T-cell-like, monocyte-like, dendritic-like, and platelet-like cells were spread across multiple metabolic neighborhoods (Supporting Information S1: Figure S2; Supporting Information S2: Table S2), supporting the view that metabolic heterogeneity is a continuous cellular property. scMetaTraj is implemented as an R package compatible with standard Seurat-based preprocessing workflows. It accepts Seurat objects or gene-by-cell expression matrices and exposes modular functions for metabolic scoring, embedding, clustering, pseudotime inference, cluster profiling, trend analysis, and visualization. For each cell and metabolic module, scMetaTraj computes the arithmetic mean of log-normalized expression values across detected genes assigned to that module. Modules with fewer than three detected genes are excluded. This design mitigates dropout by aggregating information across multiple genes, reducing sensitivity to single-gene sparsity, and avoiding unstable scores from poorly covered modules. The current module collection covers glycolysis, lactate metabolism, tricarboxylic acid cycle, oxidative phosphorylation, fatty acid oxidation, lipid biosynthesis, amino acid metabolism, nucleotide metabolism, pentose phosphate pathway, redox homeostasis, mitochondrial function, and mitochondrial translation. Higher-order axis scores are computed as secondary averages across related modules. Principal component analysis (PCA) of the cell-by-module score matrix defines the analysis space for graph construction, clustering, and pseudotime inference. Uniform Manifold Approximation and Projection (UMAP) is used only for visualization. A k-nearest-neighbor graph is built in PCA space, metabolic subclusters are identified by Louvain community detection, and mPT is calculated as the rescaled Dijkstra shortest-path distance from a selected root cell on the weighted graph. To assess root sensitivity, we compared four root-selection strategies on the PBMC reference dataset. Orientation-invariant pairwise correlations between resulting pseudotime orderings ranged from 0.28 to 1.00 in absolute value, indicating that root choice can materially influence detailed pseudotime ordering rather than implying complete root invariance (Supporting Information S1: Figure S1). The reference analysis used the 10x PBMC 3k dataset (2700 cells after preprocessing). To address the reviewers' request for external validation, we added an independent GEO dataset (GSE226488) containing resting PBMCs and PBMCs stimulated with lipopolysaccharide (LPS) for 4 h. We analyzed GSM7077865 (resting PBMC replicate 1) and GSM7077866 (LPS-stimulated PBMC replicate 1), imported the filtered feature-barcode matrices, performed standard Seurat preprocessing, annotated major PBMC populations by cluster-level marker scoring, and extracted the monocyte-lineage subset for targeted metabolic trajectory analysis. We further compared scMetaTraj metabolic pseudotime with Monocle 3 pseudotime in the GSE226488 monocyte subset. Because pseudotime orientation is arbitrary, the absolute Spearman’s rank correlation was interpreted as the primary concordance metric. The added GEO dataset contributed 14,771 PBMCs in total, from which a targeted monocyte-lineage subset was retained for focused analysis after annotation. scMetaTraj detected a significant condition-associated redistribution along metabolic pseudotime, with median mPT increasing from 0.530 in resting monocytes to 0.661 after 4-h LPS stimulation (Wilcoxon P = 6e-04; Supporting Information S1: Figures S3 and S4). At the single-module level, glycolysis, oxidative phosphorylation, tricarboxylic acid cycle, and lactate metabolism showed directionally consistent but statistically modest mean shifts, suggesting that mPT captured a multivariate redistribution across the metabolic manifold even when no single module dominated the condition contrast by itself. Representative module trends along mPT are shown in Figure 2A and Supporting Information S1: Figure S5. Dynamic metabolic remodeling along metabolic pseudotime. (A) Representative module trends along metabolic pseudotime, with dashed vertical lines indicating inferred transition regions for glycolysis, lactate metabolism, oxidative phosphorylation, pentose phosphate pathway, and the tricarboxylic acid cycle. (B) Cluster-stratified module trends across metabolic pseudotime, showing within-cluster remodeling patterns for the same representative modules. Monocle 3 pseudotime showed only a weak rank correlation with scMetaTraj mPT (Spearman’s rho = −0.14), indicating limited concordance between metabolism-centered and transcriptome-centered orderings in this dataset. We, therefore, present this analysis as an initial reference comparison rather than a systematic benchmark (Supporting Information S1: Figure S6). We additionally benchmarked runtime and memory changes across 1,000, 2,000, and 4000 cells. Scoring and embedding remained fast, whereas graph-based clustering and pseudotime inference scaled more noticeably with cell number, as expected for nearest-neighbor graph construction and shortest-path calculations. We also added PBMC cell-type overlays, PBMC root-sensitivity plots, GSE226488 validation figures, a Monocle 3 comparison plot, and structured supporting information tables (Supporting Information S1: Figure S7). scMetaTraj infers metabolic state from transcript-derived module scores rather than direct measurements of flux, metabolite abundance, or enzyme activity. The resulting mPT should therefore be interpreted as an ordering of transcriptionally supported metabolic programs rather than a direct measurement of biochemical flux. The new GEO validation analysis strengthens generalizability beyond the original PBMC3k reference, but the present revision still includes only one external perturbation-state validation dataset. Additional validation in tumor microenvironments, developmental series, and multimodal metabolic assays, together with a more systematic cross-tool benchmark, will further improve biological interpretation. Future work should also compare transcript-derived metabolic trajectories with orthogonal evidence from metabolomics, proteomics, or flux-estimation frameworks. We thank members of the research community on RedNote and GitHub for their valuable feedback, suggestions, and discussions during the development and testing of scMetaTraj, which greatly helped improve the usability and robustness of the package. During the preparation of this manuscript, ChatGPT (version 5.2, OpenAI) was used for language editing and improvement of grammar and readability. The authors reviewed and edited the output as necessary and take full responsibility for the content of the manuscript. This work was supported by the Prevention and Control of Emerging and Major Infectious Diseases—National Science and Technology Major Project (Nos. 2025ZD01902500; 2025ZD01902501-4). The authors declare no conflicts of interest. The scMetaTraj package source code, metabolic gene sets, revision scripts, and revision outputs are available in the project repository. The new validation dataset analyzed in this revision was obtained from GEO under accession GSE226488. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.