Pathway-Level Prediction and Gene Anchoring in PDAC Identify a TST-Linked Fatty-Acid/Glutathione Axis

1. Introduction

Pancreatic ductal adenocarcinoma (PDAC) has been recognized as a cancer with high mortality and limited durable benefit from available treatments [1,2]. The serum marker CA 19-9 has been used widely for monitoring, yet diagnostic and prognostic precision-particularly in early disease-has remained limited [3]. Broader biomarker efforts, including proteomic and metabolomic panels and liquid-biopsy assays, have shown incremental advances but have not consistently supported treatment selection in practice [4,5,6].

At the biological level, PDAC has been reported to rely on reprogrammed metabolism, with recurring roles for glutamine and lipid pathways that sustain growth under nutrient stress [7,8]. The tumor microenvironment has also been described as fibrotic and immunosuppressive, a state linked to weak antitumor immunity and modest clinical responses [9,10,11]. Single-cell studies have further identified antigen-presenting cancer-associated fibroblasts that shape local immunity [12].

Bulk RNA-seq permits these dimensions to be examined in parallel: immune-cell composition can be approximated by digital cytometry, and pathway activities can be summarized by gene-set scoring, yielding complementary views of tumor-intrinsic programs and microenvironmental context [13-16]. Findings from The Winter Lab have further shown that PDAC cells maintain survival under nutrient and redox stress, and that redox-support pathways influence treatment response; related mechanistic reports have highlighted links between glutamine use and antioxidant defenses in sustaining growth under hostile conditions [8,17,18]. Yet, the broader clinical significance of these dependencies across large, heterogeneous patient cohorts has not been characterized.

A practical question has remained open: which immunometabolic pathways in PDAC can be predicted consistently from routine transcriptomes, and how can those predictions be tied to specific genes in a way that is interpretable and transferable across cohorts? To address this, a compact, pathway-level map was produced that quantifies which pathways are recoverable from bulk expression, measures the added value of immune-composition context, and anchors each pathway to specific gene signals for transparent interpretation.

2. Materials and Methods

2.1. Data Acquisition and Quality Control

2.1.1. Primary Dataset

Bulk RNA-Seq data for 178 primary PDAC samples were downloaded from the AACR Project GENIE database [19] and quantified with RSEM [20]. Six outlier samples were removed based on interquartile-range criteria applied to total read counts and the fraction of zero counts per sample. Genes with missing or duplicate identifiers were collapsed or dropped, and expression values were transformed using log₂(x+1). To reduce noise, genes with zero counts in more than 50% of samples or with variance in the bottom 20^th percentile were filtered out, resulting in the main (GENIE) dataset for model training and analysis.

2.1.2. Validation Dataset

Bulk RNA-Seq count data for 357 primary PDAC samples were downloaded from the Gene Expression Omnibus under accession GEO:GSE71729 [21]. To ensure equitable comparison of results, the GSE71729 cohort was processed through the same QC pipeline described above for GENIE (outlier removal by IQR on library size and zero-count fraction; log₂(x+1) transformation; removal of low-variance and sparsely expressed genes), resulting in the external validation (GSE71729) dataset.

2.1.3. Cross-Cohort Harmonization

For further like-to-like comparison of models’ results, both primary (GENIE) and external validation (GSE71729) datasets were intersected (6,432 common genes). Raw counts were then normalized and variance-stabilized in one step using DESeq2’s size-factor normalization and variance-stabilizing transformation [22]. Finally, residual batch effects between cohorts were corrected with ComBat [23,24]. Post-QC cohort characteristics (sample/gene counts, library sizes, stage, and metastasis distributions) are summarized in Table 1 (Results).

2.2. Tumor Microenvironment Estimation and Demographic Stratification

Immune cell fractions were inferred using CIBERSORTx [14] with the LM22 signature matrix [16] (LM22 models immune lineages and does not include fibroblasts/CAFs). Output scores were normalized across the 22 LM22 cell types to sum to 1.0 per sample. Eight cell types were prespecified for downstream stratification: naïve B cells, memory B cells, CD8⁺ T cells, CD4⁺ memory resting T cells, regulatory T cells (Tregs), resting NK cells, M1 macrophages, and M2 macrophages. Fractions were summarized as medians with interquartile ranges (IQRs) at the cohort level and across AJCC stage and sex; “predominance” refers to the cell type with the highest cohort-wide median estimated fraction (ties broken by the larger 75th percentile). As a sensitivity check, exploratory non-parametric comparisons across AJCC stage (Kruskal–Wallis) and by sex (Mann–Whitney U) were performed with Benjamini–Hochberg false discovery rate control (two-sided, q < 0.05) [25,26,27].

2.3. Metabolic and Immune Pathway Activity

Pathway activity was quantified using GSVA [15] on gene sets from Gene Ontology Biological Process (GO-BP) [28], KEGG [29], Reactome [30], and the MSigDB Hallmark collection [31]. Candidate pathways were ranked by between-tumor variance of GSVA scores. The five immune pathways with the highest variance were retained. Ten metabolic pathways for downstream analyses were then selected a priori according to three criteria: (1) robust between-tumor variability, (2) established relevance to PDAC progression or therapy resistance, and (3) evidence of crosstalk with immune programs (assessed via correlations with immune-pathway scores). Anchoring to prior PDAC literature on pathway relevance was used for justification of the selected metabolic set [32,33,34,35,36,37,38,39].

2.4. Machine Learning Framework

2.4.1. Model specification, training, and ensembling

Models were constructed to predict each of the 15 pathways. For each pathway, the response variable was defined as its GSVA score. Predictor sets were composed of: (1) GSVA scores of the remaining 14 pathways (the target excluded), (2) the 22 LM22 immune cell fractions estimated by CIBERSORTx, and (3) clinical covariates (AJCC stage, metastasis state, sex). After stratification by stage (and, where applicable, metastasis and sex), an 80/20 train/test split was applied. Hyperparameters for Elastic Net (EN), Random Forest (RF), and XGBoost (XGB) were tuned by five-fold cross-validation on the training set.

Predictions from the three base learners were also averaged to form an averaged ensemble; stacking was evaluated but did not surpass averaging and is not emphasized. For each pathway and cohort, the winner was defined as the single base learner achieving the highest held-out test R². For consistent comparison within and across cohorts, both the winner’s test R² and the averaged-ensemble test R² were reported, and per-model means across pathways (EN/RF/XGB) were summarized to contextualize pathway-specific winners. All modeling was implemented in Python using scikit-learn [40] and XGBoost [41].

2.4.2. Performance evaluation and benchmarking

Model accuracy was measured by the coefficient of determination (R²) on the held-out test set, with mean ± standard error reported across cross-validation folds.

2.4.3. External validation

Trained models were applied to the GSE71729 cohort for external validation [21].

2.4.4. Controlled feature-set comparison (GSVA-only, TME-only, Integrated)

A controlled feature-set comparison was conducted to quantify the contribution of tumor microenvironment (TME) features to pathway-activity prediction. Three predictor configurations were evaluated per pathway: GSVA-only (the 14 non-target pathway GSVA scores), TME-only (the 22 LM22 immune fractions), and Integrated (GSVA-only ∪ TME-only), with the same clinical covariates (AJCC stage, metastasis state, sex) included in all configurations. The identical split, stratification, preprocessing (z-scoring within the training set and applying the parameters to the test set), learners (Elastic Net, Random Forest, XGBoost), and five-fold cross-validation (described in 2.4.1) were used without modification. For each configuration, models were trained and evaluated on the held-out test set; per-pathway test R² values were recorded for EN, RF, XGB, and the averaged ensemble, and summary statistics were computed across pathways within each configuration. The incremental benefit of integration was quantified as ΔR² = R²_Integrated− max(R²_GSVA-only, R²_TME-only), calculated for both the per-pathway winner (the base learner with the highest held-out test R²) and the averaged ensemble.

2.5. Permutation Importance

In addition to permutation importance, model-internal importances were computed per learner (absolute coefficients for EN, mean decrease in impurity for RF, and gain for XGB), normalized within learner and pathway, and averaged across learners to produce consensus rankings used in downstream summaries. Permutation importance was used to quantify each feature’s contribution to out-of-sample performance [42]. Within every 5-fold cross-validation split, models were fit on the training fold and a baseline R² was recorded on the held-out fold. For each predictor (pathway GSVA score, CIBERSORTx fraction, or clinical covariate), its values in the held-out fold were randomly permuted 10,000 times, the model was re-evaluated, and ΔR² = R²_perm − R²_base was stored (more negative ΔR² indicates greater importance). Mean ΔR² across permutations and folds with 95% CIs are reported, and two-sided empirical p-values are computed as the fraction of permutations with ΔR² ≥ 0, adjusted across features per pathway using Benjamini–Hochberg FDR control (q < 0.05) [25].

Permutation importance was computed in two modes: (i) a single-feature mode, in which each predictor was permuted individually in the held-out fold to obtain its marginal ΔR²; and (ii) a groupwise mode for correlated predictors (|r| ≥ 0.7), in which the entire correlated set was permuted jointly to obtain a correlation-aware, joint ΔR², mitigating bias from multicollinearity.

2.6. Gene-level correlation analysis

For each pathway, Pearson correlations were computed between each gene’s log₂(x+1) expression and that pathway’s GSVA score across all samples. The top 10 genes per pathway were retained by the absolute correlation coefficient (|r|). Correlations were computed genome-wide and were not restricted to model features. This analysis was used to contextualize model-derived feature importance at the gene level; genes were not used as predictors in the models.

3. Results

3.1. Post-QC Cohort Characteristics

After QC, the GENIE expression matrix comprised 172 high-quality PDAC samples and 6,432 genes. Post-QC library sizes in GENIE averaged 1.38 × 10⁵ reads per sample (SD =1.34 × 10⁴; range, 1.26 × 10⁵–1.64 × 10⁵). Among GENIE patients, metastasis state was M0 (no distant metastasis) in 84 samples (45.7%), M1 (distant metastasis) in 5 samples (2.7%), and MX (unknown) in 95 samples (51.6%) (Table 1). Figure 1 shows the primary dataset (GENIE) before and after quality control. Multiple PCA outliers were identified; all extreme outliers were removed, while a single marginal outlier-located near the main cluster-was retained for downstream analyses.

3.2. Tumor Microenvironment Composition and Demographic Stratification

Digital cytometry (CIBERSORTx/LM22) across 172 PDAC tumors indicated M2 macrophages as the predominant immune population, while CD8⁺ T cells and regulatory T cells (Tregs) constituted the principal lymphoid compartments (Figure 2A). Figure 2B stratifies eight representative immune cell types by AJCC pathologic stage and sex; groupwise medians suggested CD8⁺ levels peak at stage IIB and Tregs increase with stage, with no material sex-specific differences. Exploratory stage- and sex-based comparisons did not remain significant after FDR control (all q >0.05).

3.3. Predicting pathway activity from gene expression

Pathway variability and selection were defined a priori as described in Methods (2.3); The accompanying rationale and figures are provided in Supplementary Text S1 and Supplementary Figures S1–S3.

3.3.1. Within-cohort model performance (primary cohort)

Higher accuracy was observed for metabolic targets than for immune targets.

• Metabolic (n = 10): the winner mean test R² was 0.825 (EN was the per-pathway winner for 10/10 targets); the averaged-ensemble mean was 0.714 (Δ = 0.111). Per-model means (test R²) were: EN 0.825, XGB 0.744, RF 0.572.
• Immune (n = 5): the winner mean test R² was 0.648 (winners distributed as XGB = 2, RF =2, EN = 1); the averaged-ensemble mean was 0.603 (Δ = 0.045). Per-model means were: XGB 0.622, RF 0.598, EN 0.589.

Pathway-level performance distributions are shown in Figure 3A (boxplots for all 15 pathways, sorted by the mean test R² across base learners; metabolic = red; immune = blue). Model-by-model bar plots are provided in Supplementary Figure S4.

3.3.2. External validation

Generalizability was assessed on the external validation cohort GSE71729 (n = 351) using identical preprocessing, feature construction, and the same winner rule; the metabolic-over-immune advantage was preserved.

• Metabolic (n = 10): the winner mean test R² was 0.685 (EN remained the per-pathway winner for 10/10 targets); the averaged-ensemble mean was 0.584 (Δ = 0.101). Per-model means were: EN 0.685, XGB 0.544, RF 0.522.
• Immune (n = 5): the winner mean test R² was 0.562 (winners mixed across models); the averaged-ensemble mean was 0.543 (Δ = 0.019). Per-model means were: XGB 0.562, RF 0.538, EN 0.529.

Validation distributions are shown in Figure 3B; pathway-wise bar plots are provided in Supplementary Figure S5.

3.4. Integration of GSVA and TME features

Model performance was compared across GSVA-only, TME-only, and Integrated predictor sets using the identical pipeline. In Figure 4, each panel (A) metabolic, (B) immune displays three box-and-whisker plots summarizing the distribution of test R² values for the corresponding configuration; model identities are not shown. For metabolic targets, the Integrated configuration shows a modest upward shift in median and upper quartile relative to GSVA-only and clearly exceeds TME-only. For immune targets, the Integrated configuration shows a marked upward shift compared with either single-source configuration, while the GSVA-only and TME-only distributions remain lower. These results indicate that combining GSVA and TME features yields higher predictive accuracy than using either feature set alone.

3.5. Feature Importance and Permutation Testing

Feature contributions were assessed by permutation importance (as described in Methods 2.5), with higher scores indicating larger performance loss when a feature was perturbed. In Figure 5, metabolic features (Figure 5A) are shown in red and immune features (Figure 5B) in blue; clinical covariates appear at the right of each panel. Among metabolic terms, broad small-molecule and lipid programs dominated. Monocarboxylic-acid metabolic process and sphingolipid metabolic process exhibited the highest median importances with wide upper tails, followed by dicarboxylic-acid metabolic process, steroid biosynthetic process, and cellular modified amino-acid biosynthetic process. Redox and amino-acid–linked terms-including NADPH regeneration, L-amino-acid metabolic process, serine-family amino-acid biosynthetic process, and arginine catabolic process- showed consistent but slightly lower medians. Clinical covariates ranked near the bottom.

As for immune feature patterns, the largest and most consistent importances were observed for processes linked to cytotoxic and antiviral immunity, including positive regulation of natural killer (NK) cell–mediated cytotoxicity, positive regulation of NK cell–mediated immunity, and granzyme-mediated programmed cell death. Pathways reflecting antigen-clearance and inflammatory signaling, including Fc-gamma receptor signaling, defense response to virus, cellular response to type I interferon, and positive regulation of complement activation, also showed elevated median importances. T-cell context features such as CD4⁺ αβ T-cell activation and mononuclear cell differentiation contributed at intermediate levels. Clinical covariates (stage, sex, metastasis) displayed low median importances relative to immune pathway features.

Together, these patterns indicate that predictive signal for metabolic targets relied most on lipid handling and core small-molecule/amino-acid metabolism, whereas for immune targets it was concentrated in NK/cytotoxic and interferon/antiviral programs, with clinical covariates contributing minimally. These feature-level patterns were concordant with gene–pathway correlations (in 3.6; e.g., lipid and small-molecule metabolism for metabolic targets; cytotoxic/interferon programs for immune targets), while exact one-to-one correspondence was not expected due to aggregation and multicollinearity.

3.6. Gene-Level Correlations

Gene–pathway associations were summarized using Pearson correlations between gene expression and pathway GSVA scores (top 10 by |r| per pathway), with the heatmap shown in Figure 6 (red = positive; blue = negative). Among metabolic pathways, the strongest association was observed for SC4MOL with the sterol biosynthetic process. TST then showed equally strong positive correlations with both the fatty-acid biosynthetic process and the glutathione metabolic process. Additional positive correlations were observed for AGMAT with the arginine metabolic process, HOOK1 with the glutamate metabolic process, and TMEM52 with the serine-family amino-acid metabolic process. Negative correlations included a strong association of ZNF25 with the glucose-6-phosphate metabolic process, as well as C14orf181 with the sphingolipid biosynthetic process, FGD1 with the one-carbon metabolic process, and HLX with the glutamine-family amino-acid metabolic process.

For immune pathways, pathway-concordant patterns were reproduced, and all correlations were positive: PRF1 with the CD8-activation signature, CXCL10 with IFNG signaling, MRC1 (CD206) with macrophage activation, GNLY with NK cytotoxicity, and FOXP3 with the Treg signature. These gene-level results were consistent with the feature-importance themes reported above and provided mechanistic anchors for the pathway readouts.

4. Discussion

A pathway-centered analysis of bulk tumor transcriptomes and immune-composition profiles was conducted to clarify which PDAC pathways are recoverable from routine data and how those signals can be anchored for interpretation. Three conclusions were supported. First, the ten preselected metabolic pathways were predicted accurately from expression in the primary cohort, most often by a sparse linear model, and this relative advantage over immune-pathway prediction was reproduced in external validation (Figure 3). Second, the five immune pathways showed mixed difficulty, with improvements observed when immune-composition features were integrated with expression, indicating that context from the tumor microenvironment added information beyond expression alone (Figure 4). Third, the cohort displayed diverse immune composition without cohort-level shifts by stage or sex (Figure 2), reducing concern that broad composition differences explain the modeling contrasts between metabolic and immune targets.

Feature-importance profiles and permutation-based robustness checks supported these findings (Figure 5). The gene-level analysis provided concrete anchors that connected model behavior to recognizable biology (Figure 6). In particular, TST emerged as a cross-pathway signal linking fatty-acid biosynthesis and glutathione metabolism, consistent with coupling between lipid anabolism and antioxidant capacity in PDAC. This observation aligns with prior work indicating that survival under nutrient and redox stress is supported by metabolic programs with antioxidant components; such context helps interpret why redox-linked nodes may track with lipid-anabolic demands in tumors.

Limitations and future direction: Bulk RNA-seq infers pathway activity indirectly and does not measure flux or resolve cell-type–specific effects. Single-cell or spatial assays, together with metabolomic and flux measurements, would be expected to refine immune-pathway models and clarify mechanism. These constraints do not alter the central observations that metabolic pathway activities were recovered with high accuracy from bulk expression, that immune pathway activities were more context dependent, and that integration with immune-composition information improved immune-pathway prediction.

5. Conclusions

A compact, pathway-level map of PDAC was produced in which (i) metabolic pathways were predicted accurately and reproducibly across cohorts, (ii) immune pathways showed variable but informative performance that improved with immune-composition context, and (iii) feature-importance, permutation testing, and external validation supported robustness. Gene-level anchors-including TST as a top signal shared by fatty-acid biosynthesis and glutathione metabolism-provide concrete entry points for follow-up studies and potential translational use.