The training and investigation dataset was obtained from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). Two additional datasets were obtained from the Gene Expression Omnibus (GEO; GSE62452) (https://www.ncbi.nlm.nih.gov/gds/) and International Cancer Genome Consortium (ICGC; https://icgc.org/) databases that served as external validation sets (Figure 1). All three datasets are publicly available and accessible online. The reasons for selecting these datasets for study were that they originated from PC patients, encompassed survival and transcriptional information, and boasted a sufficient sample size. After downloading the data, the following three steps were performed during the initial data processing: Patients with missing or zero overall survival (OS) were excluded, the transcriptome data from TCGA and ICGC databases were unified as log₂ (x+1), and gene IDs were converted to gene symbols.

To assess the activity of different phenotypes, co-expressed genes were initially selected both in TCGA and GEO. To eliminate non-significant genes, DEGs between tumor and non-tumor tissues were then screened using the “limma” package in GEO, with a false discovery rate (FDR) < 0.05 as the cutoff. The intersection gene dataset, both in co-expressed genes and DEGs, was taken for subsequent analysis. The activities of different metabolic pathways and DDR for each sample were scored using the single-sample gene set enrichment of the R package GSVA in TCGA, utilizing the expression data of the intersection gene dataset. The metabolic signature gene set, which includes seven metabolic pathways (lipid, carbohydrate, amino acid, integrated energy, nucleotide, vitamin cofactor, and TCA), was previously reported (Sung and Cheong, 2021). In total, 276 DDR gene sets were obtained from previous studies and used for the knowledge-based curation of DDR pathways (Li et al., 2022). Immune infiltration analysis included the immune cell profiling (ICI) and estimate algorithm with the intersection gene dataset. ICI involved 22 types of immune cells using the “CIBERSORT” package, while the estimate algorithm included four types of scoring for immune, stromal, estimate, and tumor purity using the “ESTIMATE” package (Liu et al., 2022).

WGCNA was used to identify potential hub genes that share similar functions by constructing a co-expression network analysis and mining modules (Su et al., 2021). We applied WGCNA to define gene modules associated with metabolic phenotypic features derived from the intersection gene dataset in TCGA. Specifically, samples with outliers were removed first, and the function power estimate was used to define the optimal soft threshold. The subsequent steps and parameters of WGCNA included constructing correlation networks, setting the maximun number of gene modules to 30, and calculating the dissimilarity of the module eigengenes. The similarity cutoff was set to 0.75, and the modules with the top three high correlation coefficients were considered most relevant to metabolic features. These modules were merged for downstream analysis (Xiao et al., 2022).

For the training set, we used patient data from TCGA, whereas the validation sets consisted of data from GEO and ICGC. Univariate Cox regression using the “survival” package was employed to preliminarily filter the genes linked to the prognosis from the aforementioned metabolism-associated hub genes, with a p-value <0.001. Subsequently, LASSO regression with the “glmnet” package was used to address overfitting among the survival-related hub genes. Multivariate Cox stepwise regression was applied using the “survival” package to construct the prognostic model with the highest discriminatory ability for predicting OS. The equation for the risk score, which predicts prognosis based on metabolic hub genes, is shown in Figure 3A. The acquired formula for calculating the risk score was as follows: Risk score = (expression of DLX3 × 0.114) + (expression of HMGA2 × 0.144) + (expression of SPRR1B × 0.071) + (expression of MYEOV × 0.104) + (expression of FAM111B × 0.173).

Each patient was assigned a risk score and then categorized into high- and low-risk groups based on the median value independently within each cohort. To validate the accuracy of the model, survival curves were generated using the Kaplan–Meier (KM) method for the two groups separately in three cohorts. Additionally, the receiver operating characteristic (ROC) curve and area under the curve (AUC) values were calculated using the “survival ROC” package to evaluate the time-dependent predictive ability of the model. Furthermore, we compared our model with nine reported models from the literature in terms of their ability to predict 1-year survival, using ROC and AUC analysis (Ma et al., 2021; Ye et al., 2021; Huang et al., 2023; Zhang et al., 2023; Ren et al., 2021; Liu et al., 2021; Yan et al., 2020; Yan et al., 2022; Zhou et al., 2019).

We compared the differences between the two groups in terms of available characteristics, including age, stage, tumor site, therapies, alcohol history, and the expression of the five indicative genes in TCGA. The comparison results were presented using boxplots or heatmaps generated by the “ggpubr” and “Complex Heatmap” packages.

To explore the underlying mechanism of the risk score, we analyzed the association between risk scores and phenotypic features, including metabolism, DDR, and immune infiltration. This analysis involved comparing the differences in seven metabolic pathways, 22 immune cell infiltrations, and four types of estimates between the two different risk groups. Additionally, we analyzed the correlation between significantly different phenotypic features. The results were visualized using violin charts, beeswarm plots, or correlation heatmaps generated by the “ggplot2,” “ggbeeswarm,” and “corrplot” packages.

To further investigate the implied mechanisms of the risk score, gene set enrichment analysis (GSEA) was performed on DEGs between two risk groups. The DEGs were initially screened using the “DESeq2” package, with a |log2FC| > 1 and an adjusted FDR <0.05 as the cutoffs. The direction of regulation (up or down) for the DEGs was plotted using the “ggplot2” package. DEGs with an adjusted FDR <0.05 were then used for GSEA. Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and GSEA were performed using the “Clusterprofiler” package with reference to “org.Hs.eg.db” (Carlson M, 2019. org.Hs.eg.db: Genome wide annotation for Human. R package version 3.8.2.). An adjusted p-value <0.05 was selected as the cutoff criteria for considering significant results in GO terms, KEGG pathway, and GSEA.

The drug sensitivity of patients was predicted using the half-maximal inhibitory (IC50) calculated with the “pRRophetic” package. We selected several recommended or associated drugs for drug sensitivity analysis based on the clinical guidelines of PC. The difference between the two groups was assessed and displayed using columnar beeswarm plots created with the “ggbeeswarm” and “customLayout” packages. Additionally, we downloaded data on patient response to gemcitabine in TCGA using the “TCGAbiolinks” package, as reported by Nicolle, and validated it between the two groups (Nicolle et al., 2021).

To demonstrate the practical application, a clinical nomogram model combining the risk score and clinical parameters was constructed using the “survival” and “RMS” packages. Owing to the large number of missing values for clinical parameters in TCGA, only variables (age, stage, radiotherapy, and pharmatherapy) with a sufficient number of samples and high relevance to survival were included. The nomogram was used to evaluate the 1-, 2-, and 3-year survival rates of patients in TCGA. Calibration curves and the C-index were used to assess the consistency between the predicted and actual survival rates.

Statistical analysis was performed using GraphPad Prism (version 9.0.0 for Windows) and R software (version 4.2.1 for Windows). In GraphPad Prism, t-tests or Wilcoxon signed-rank tests were used to compare measured differences between groups, and the chi-square test was used to analyze associations between categorical variables. Other statistical analyses were primarily conducted using R software, as described earlier, including scoring each feature, identifying DEGs and hub genes, GSEA, calculating drug sensitivity, and constructing the nomogram. The correlation coefficient between phenotypic features was evaluated using Pearson’s correlation. Statistical significance was set at P < 0.05.

We obtained transcriptome and clinical data from 178 PC patients (via TCGA) in this study for constructing models and exploring mechanisms. Furthermore, we obtained an additional 159 PC cases (69 in GSE62452 and 90 in ICGC-PAAD-US) for validating the risk model. To screen DEGs in GSE62452, we collected 61 non-tumor tissues from PC patients. We identified 10,031 genes that existed in both TCGA and GSE62452, and 6,747 genes were screened as DEGs in GSE62452. For further analysis, we obtained 5,872 genes as the intersection gene dataset, and their transcriptome data were used to score metabolism, DDR, and ICI activities for each patient.

The function power analysis showed that a soft threshold of 14 was optimal for ensuring a scale-free network in the WGCNA (Figure 2A). By merging similarity modules and applying a cut-line, we identified a total of 13 modules, with each module assigned a color using the dynamic tree cut method (Figure 2C). The coexpression network was visualized in a network heatmap plot (Figure 2B). Each color module represented a gene set with closely associated expressions. The gray module represented genes that were not assigned to any other cluster. The other 12 modules were plotted in a heatmap to evaluate the association between each module and different metabolic pathways (Figure 2D). Based on the criteria (|Cor| > 0.65), yellow, black, and green-yellow modules, which collectively included 536 genes, were selected for further analysis.

To filter genes with predictive value, we performed univariate Cox regression analysis and identified 362 genes among the metabolic hub genes (p < 0.001). Based on the 1000-fold cross-validation results of LASSO regression analysis, the optimal efficacy was achieved when ten variables were selected as the target markers (Figures 2E,F). Furthermore, through multivariate Cox regression analysis, we identified five genes derived from the LASSO regression analysis that could independently construct a gene signature to predict the prognosis of PC as described in 2.4.

Patients were divided into two groups based on the median risk score for validation. The KM curve showed that the high-risk group had significantly worse OS than the low-risk group (Figure 3B), which was further verified in external validation sets from GEO and ICGC (Figures 3C,D). The ROC curves and AUC values of the three datasets demonstrated the predictive ability of the risk score from a time-dependent perspective (Figures 3E–G). In TCGA, the AUCs for 1-, 2-, and 3-year were 0.75, 0.74, and 0.80, respectively, indicating good predictive performance for PC. When compared with other reported models, our model exhibited a better ROC curve and AUC value, and the top five models are displayed in Figure 3H (Supplementary Figure S1).

The distribution of clinical and pathological information between the two groups is presented in Table 1. The analysis results showed that, except for vital status and pathological types, no significant distribution differences were observed in other characteristics. The expression levels of the five indicator genes in the high-risk group were higher than those in the low-risk group, which may explain the potential mechanisms of OS (Figure 4A).

The risk score showed strong associations with metabolism, DDR, and estimate scores but weak correlations with 22 immune cells, as revealed by the differential analysis (Figures 4B–D). The high-risk group demonstrates lower scores of stromal, immune, and estimate, while the tumor purity score is higher. In terms of metabolism, the high-risk group exhibits higher scores for TCA, DDR and metabolism of lipid, carbohydrate, amino acid, nucleotide, vitamin, but a lower score for integrated energy metabolism. To further investigate the implied mechanism of the risk score, correlation analysis revealed that the risk score had the strongest correlation with DDR rating scores (correlation coefficient: 0.41, P < 0.001), followed by nucleotide, carbohydrate, immune, TCA, and vitamin scores (Figures 4E,F). In terms of the relationship between phenotypic features, DDR showed a strong correlation with lipid, carbohydrate, amino acid, nucleotide, and TCA metabolism, particularly with the latter two features. In contrast, the immune score displayed a weak correlation with metabolism and DDR. For the internal analysis of metabolism pathways, the correlation coefficient between nucleotide and TCA was 0.72, and both were closely related to lipid, carbohydrate, and amino acid metabolism. The correlation between different factors showed a positive correlation, except for energy metabolism. This reflects that patients identified as high-risk by our model may be more closely related to changes in the DDR pathway, which requires the support of metabolic pathways such as nucleotide and TCA.

In total, 10,720 DEGs were re-identified as potential factors that could contribute to the survival differences between the two different risk groups. Among them, 705 DEGs were upregulated and 2,660 DEGs were downregulated in the high-risk group compared with those in the low-risk group (Figure 5A). GO analysis of the DEGs indicated that these genes were highly enriched in biological processes such as nucleocytoplasmic transport or replication, cellular components including mitochondria, and molecular functions related to regulating metabolism or ATP-dependent activities, as shown in Table 2 (Figure 5D). KEGG analysis and GSEA results showed that DEGs were enriched in several pathways closely related to metabolism and proliferation (Figures 5B,C,E). Nine pathways may explain the factors associated with the risk score and their effect on survival; among these pathways, five were upregulated and four were downregulated (Figures 5F,G).

IC50 values reflect the sensitivity of patients to drugs, with lower values of IC50 indicating greater sensitivity. Gemcitabine and 5-fluorouracil were the only pharmaceuticals included in the drug sensitivity analysis that were also recommended in the PC guidelines. Additionally, the IC50 values of cisplatin, olaparib, and paclitaxel were analyzed considering that oxaliplatin and albumin paclitaxel were recommended in the PC guidelines. Furthermore, the predictive IC50 values of these five pharmaceuticals to assess the predictive value of the risk score between the two groups. The results showed that, except for paclitaxel, the IC50 values of other pharmaceuticals in the high-risk group were significantly higher than those in the low-risk group (Figure 6A). In the analysis of the actual response to gemcitabine, a trend emerged suggesting that patients with lower risk scores and a higher proportion in the low-risk group exhibited better responses (Figures 6B,C).

To further validate the accuracy and usability of the risk score, a nomogram was constructed using the risk score and several clinical parameters. The total points obtained from the nomogram were used to evaluate the 1-, 2-, and 3-year OS rates in PC (Figure 6D). The nomogram displayed that younger age, administration of pharmaceutical therapy, and low-risk scores were all protective prognostic factors; these tendencies were consistent with expectations. The previous administration of radiotherapy and the later stage of the disease also showed a partial protective tendency. However, the sample size in some cohorts was insufficient to achieve statistical significance. A forest plot was used to display the predictive value of the five hub genes with hazard ratios in PC (Figure 6E). We further analyzed levels of five genes in other cancers with TCGAplot, most of them upregulated (Supplementary Figures S2A–E) (Liao and Wang, 2023). The protein expression of them on the Human Protein Atlas database may also upregulated in PC (DLX3 and SPRR1B pending analysis) (Supplementary Figures S2F,H). To verify the accuracy of the model, calibration curves were performed (Figure 6F), and the C-index was found to be 0.778 (0.729–0.828), indicating that the nomogram has good capability for predicting the OS rate of PC patients.