Abstract

Background:

There is a strong correlation between lactylation, programmed cell death, and the progression of cancer. This study aims to identify prognostic genes associated with lactylation and programmed cell death in pancreatic ductal adenocarcinoma (PDAC), providing new insights for risk stratification and therapeutic strategies.

Methods:

TCGA-PAAD, GSE62452, lactylation-related genes (LRGs), and programmed cell death-related genes (PCDRGs) were retrieved from relevant databases and references. Prognostic genes were identified through univariate Cox regression analysis, followed by random survival forest analysis for survival prediction. Subsequently, enrichment analysis, immune microenvironment analysis, drug sensitivity analysis, immunohistochemical analysis, and expression analysis of prognostic genes were conducted. Finally, the experimental verification was carried out in clinical samples.

Results:

In this investigation, two prognostic genes (HMGA1 and KIF2C) linked to lactylation and programmed cell death were identified, and a robust prognostic risk model was developed. Enrichment analysis results included Cell cycle, G2M checkpoint, Myogenesis, and Angiogenesis. Moreover, immature B cells and activated B cells demonstrated the strongest positive correlation (cor = 0.97, P < 0.001), while neutrophils and activated B cells demonstrated the strongest negative correlation (cor = −0.68, P < 0.001). Furthermore, KIF2C and HMGA1 demonstrated the strongest negative relationships with mast cells (correlation coefficients = −0.36 and −0.53, P < 0.01). Drug sensitivity analysis revealed that Sapitinib was more effective in the high-risk group (HRG), while Doramapimod was more effective in the low-risk group (LRG) (P < 0.0001). Both immunohistochemical and expression analyses of prognostic genes showed that HMGA1 and KIF2C were upregulated in PDAC patients (P < 0.05). Finally, genes in the clinical samples also showed the same expression trend.

Conclusion:

In the present investigation, two prognostic genes were identified, and subsequently, a predictive risk model was established, which may serve as a valuable reference for the clinical management of PDAC.

1 Introduction

Pancreatic ductal adenocarcinoma (PDAC) accounts for approximately 85% of all pancreatic cancer cases and is one of the lethal malignancies with an extremely poor prognosis (Jentzsch et al., 2020; Narayanan et al., 2021). Its 5-year overall survival rate remains low at about 9%–11%, making it a major cause of cancer-related death (Mizrahi et al., 2020; Hsu et al., 2022). Alarmingly, while the overall cancer mortality rate has decreased, the incidence of PDAC continues to rise (Siegel et al., 2024), with predictions showing that by 2030–2040, it will become the second leading cause of cancer death after lung cancer (Rahib et al., 2014; Rahib et al., 2021). The poor prognosis of PDAC is mainly attributed to early systemic spread, invasive local infiltration, and limited therapeutic efficacy (Connor and Gallinger, 2022; Garajová et al., 2023). The disease typically progresses through multiple stages from precancerous lesions to low-grade and high-grade dysplasia, characterized by progressive cytological atypia and genetic aberrations. Clinical manifestations often include non-specific symptoms such as loss of appetite, nausea, indigestion, back pain, and unexplained weight loss, while jaundice frequently appears as a late indicator. Although surgical resection remains the most effective treatment modality, only about 20% of patients meet the criteria for curative surgery at the time of diagnosis (Schneider et al., 2021; Garajová et al., 2023). Therefore, identifying new prognostic biomarkers is crucial for understanding the molecular mechanisms of PDAC and realizing individualized treatment strategies.

Lactate, as one of the most abundant metabolites in the circulatory system, is produced through the action of lactate dehydrogenase (LDH) on pyruvate during glycolysis (Rabinowitz and Enerbäck, 2020). Under hypoxic conditions, pyruvate is converted to lactate, allowing for the regeneration of NAD+, which is crucial for maintaining glycolytic flux (Li L. et al., 2023). Notably, even under normoxic conditions, cancer cells preferentially metabolize glucose into lactate—a phenomenon known as the Warburg effect—leading to abnormal protein lactylation modifications in the tumor microenvironment (Wang and Patti, 2023; Yu et al., 2024). Lactylation, a newly discovered post-translational modification, has emerged as a key regulator of gene expression and cell metabolism during cancer progression. In PDAC, the accumulation of lactate within the tumor microenvironment has been proven to drive histone lactylation both in vitro and in vivo, potentially affecting tumor cell behavior and immune evasion (Li F. et al., 2024).

Programmed cell death (PCD) encompasses multiple genetically regulated pathways, such as apoptosis, autophagy, ferroptosis, pyroptosis, and necroptosis. These five represent some of the 18 characterized forms. In PDAC, tumor cells exhibit significant differences in ferroptosis sensitivity, with some cells escaping ferroptosis by upregulating antioxidant systems. This suggests that inducing ferroptosis could be a novel therapeutic strategy (Li G. et al., 2024). Interestingly, lactate plays a dual role in cell death: it can inhibit apoptosis by creating an acidic microenvironment, thereby enhancing tumor cell survival, while, under specific conditions, lactate accumulation may paradoxically induce ferroptosis (Zhao et al., 2020; Yang Z. et al., 2023). Despite these emerging findings, the mechanistic interaction between lactylation and programmed cell death in PDAC remains incompletely elucidated, and no research has yet established a prognostic model for PDAC that integrates these two biological processes.

This study used public PDAC datasets to identify differentially expressed genes associated with lactylation and programmed cell death via differential expression analysis. Prognostic genes were identified using univariate and random forest algorithms, followed by the construction and validation of a predictive risk model. Clinical feature analysis, functional enrichment analysis, immune correlation analysis, tumor mutation burden, and drug sensitivity analysis were also performed. Additionally, the prognostic genes were subjected to chromosomal localization and immunohistochemistry analysis. Finally, experimental verification confirmed that the expression of prognostic genes was consistent with the bioinformatics results. This study aims to provide a new theoretical basis for the clinical prognostic evaluation and the formulation of treatment strategies for PDAC.

2 Materials and methods2.1 Data collection

The training set of PDAC patients was derived from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), including 143 tumor (PDAC) and four para-cancerous control tissue samples. Of these, 142 patients had survival information. Meanwhile, the clinical data and somatic mutation data were also collected. The validation set GSE62452 (GPL6244), which encompassed tumor tissue specimens from 65 PDAC patients with survival information, was derived from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The visit took place on 26 February 2025. Furthermore, a total of 332 lactylation-related genes (LRGs) were obtained from reference (Huang et al., 2023) (Supplementary Table S1), and 1,548 programmed cell death-related genes (PCDRGs) were obtained from reference (Qin et al., 2023) (Supplementary Table S2).

2.2 Differential expression analysis

To identify differentially expressed genes (DEGs) between PDAC and control samples within the training set, the “DESeq2” package (v 1.38.0) (Love et al., 2014) was employed for differential expression analysis, and DEGs were screened based on P < 0.05 and |log2 fold change (FC)| >0.5. Subsequently, the “ggplot2” package (v 3.4.1) (Gustavsson et al., 2022) was utilized to construct a volcano plot for visualizing the DEGs.

2.3 Identification and functional characterization of candidate genes

To identify the DEGs linked to lactylation and programmed cell death in PDAC, the “ggvenn” package (v 0.1.9) (Chen and Boutros, 2011) was used to intersect DEGs, LRGs, and PCDRGs; the intersection genes were used as candidate genes. Next, to explore the functional pathways related to candidate genes in PDAC, gene set enrichment analysis (GSEA) was performed. The reference gene set “c2.cp.kegg.v7.5.1.symbols” was carefully chosen from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb). First, based on all samples in the training set, Spearman correlation analyses were performed between each candidate gene and all the remaining genes in the training set, separately, using the “psych” package (v 2.2.9) (Orifjon et al., 2023) to obtain the correlation coefficients for all samples in the training set. Subsequently, genes were arranged in descending order according to these coefficients. The sorted data were then utilized to conduct GSEA (|normalized enrichment score (NES)| > 1 and P < 0.05) using the implementation of the “clusterProfiler” package (v 4.2.2) (Yu et al., 2012).

2.4 Determination of prognostic genes, development, and validation of a prognostic risk model

Within the samples of PDAC with survival information in the training set, univariate Cox regression analysis was utilized to analyze the candidate genes by the “survival” package (v 3.5.3) (Kumar et al., 2020) (hazard ratio (HR) ≠ 1, P < 0.2). Then, Random Survival Forest (RSF) analysis was conducted for the outcomes of the proportional hazards (PH) assumption test (P > 0.05) using the “randomForestSRc” package (v 3.2.3) (Ishwaran et al., 2008) to obtain risk scores for each patient.

Additionally, patients were categorized into a high-risk group (HRG) and a low-risk group (LRG) based on the optimal threshold derived from PDAC samples with survival information. The “survival” package (v 3.5.3) (Kumar et al., 2020) was used to plot the Kaplan-Meier (K-M) survival curve (P < 0.05, log-rank test). Furthermore, the risk score distribution maps and survival state distribution maps were plotted using the “ggplot2” package (v 3.4.1) (Gustavsson et al., 2022). Next, the receiver operating characteristic (ROC) curve of the prognostic risk model was plotted using the “survivalROC” package (v 1.0.3.1) (Zheng et al., 2021). The area under the curve (AUC) was calculated to evaluate the model’s effectiveness (AUC >0.6). Finally, the model’s reliability was validated in the validation set using the above method.

2.5 Clinicopathological characteristics analysis

Among the samples with survival information in the training set, the “pheatmap” package (v 1.0.12) (Wang et al., 2023c) was used to draw a heat map to show the distribution of risk scores in clinicopathological characteristics (age, gender, stage T, N, and tumor stage) and the expression of prognostic genes in HRG and LRG.

2.6 Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

To examine the biological functions and pathways among patients in distinct risk groups, the following analysis was conducted. The “c2.cp.kegg.v7.5.1.symbols.gmt” gene set was acquired from the MSigDB as the reference gene set. Within the samples of PDAC with survival information in the training set, the “DESeq2” package (v 1.38.0) (Love et al., 2014) was utilized to scrutinize the disparities between HRG and LRG to obtain the corresponding genes and log2FC values. Subsequently, the genes were arranged in descending order by log2 fold change log2FC. The “clusterProfiler” package (v 4.2.2) (Yu et al., 2012) was utilized to carry out GSEA (|NES| >1 and P < 0.05).

GSVA was conducted in HRG and LRG using the “GSVA” package (v 1.46.0) (Hänzelmann et al., 2013). The “h.all.v7.5.1.symbols.gmt” gene set was obtained from the MSigDB as the reference gene set. The ssGSEA scores of the signaling pathways were evaluated. The “limma” package (v 3.58.1) (Ritchie et al., 2015) was utilized to compare the differences in scores between groups. Generally, |t| >2 and P < 0.05 were considered to indicate a significant difference.

To quantify lactylation activity, a score was calculated via ssGSEA using a curated gene set (including LDHA/B, SIRT1/2, EP300, and GTPSCS) and compared between risk groups. This activity score was then correlated with HMGA1 and KIF2C expression to evaluate their relationship with the tumor’s lactylation-associated metabolic landscape.

2.7 Analysis of immune microenvironment

To understand the immune infiltration of different samples in the training set, the ssGSEA algorithm was employed to analyze the distribution of 28 distinct immune cells (Jiang et al., 2022) in both HRG and LRG samples. Subsequently, the Wilcoxon test was harnessed to investigate the disparities in the infiltration levels of these 28 immune cells between HRG and LRG samples (P < 0.05), and the results were visualized using the “ggplot2” package (v 3.4.1) (Gustavsson et al., 2022). Furthermore, Spearman correlation was performed to investigate the connections between the differential immune cells, as well as differential immune cells and prognostic genes (|correlation coefficient (cor)| >0.3 and P < 0.05).

Additionally, to fully understand the immune infiltration landscape of PDAC, the stromal score, immune score, and ESTIMATE score were computed using the “estimate” package (v 1.0.13) (Wang Y. et al., 2022) in the training set, which included tumor samples with survival information. And differences in scores between HRG and LRG were assessed by the Wilcoxon test, with P < 0.05. Moreover, Spearman correlations were carried out to investigate the relationships between the risk scores and ESTIMATE scores, immune scores, and stromal scores (|cor| > 0.3 and P < 0.05) using the “psych” package (v 2.2.9) (Orifjon et al., 2023).

2.8 Drug susceptibility prediction and gene mutation landscape analysis

In the tumor samples with survival information in the training set, PDAC chemotherapeutic agents were retrieved from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org). Moreover, the half-maximal inhibitory concentration (IC50) values for the HRG and LRG were estimated using the “pRRophetic” package (v 0.5) (Wang et al., 2023b) to evaluate the drug susceptibility. The sensitivities of the HRG and LRG were contrasted by means of the Wilcoxon test, with P < 0.05.

Furthermore, the mutation frequencies of genes in HRG and LRG were analyzed using the “maftools” package (v 2.14.0) (Mayakonda et al., 2018), and the top 20 mutation frequencies were visualized using a waterfall plot. Based on the somatic mutation numbers of PDAC patients, the “maftools” package (v 2.14.0) (Mayakonda et al., 2018) was utilized to calculate the tumor mutational burden (TMB) scores of the patients. The Wilcoxon test was employed to contrast the disparities in TMB scores between the HRG and LRG (P < 0.05). Meanwhile, Spearman correlation was performed to investigate the relationship between the risk and TMB scores (|cor| >0.3 and P < 0.05) using the “psych” package (v 2.2.9) (Orifjon et al., 2023).

2.9 Prognostic gene analysis

To map the positions of prognostic genes on chromosomes, the “RCircos” package (v 1.2.2) (Hu et al., 2014) was implemented.

To study the expression distribution of proteins encoded by prognostic genes in tumor and control tissues, immunohistochemical analysis of prognostic genes in PDAC tissues and control tissues was performed using the Human Protein Atlas (HPA) database (https://www.proteinatlas.org). Additionally, the expression variations of prognostic genes were explored using the Wilcoxon test in PDAC and control samples from the training set (P < 0.05).

2.10 Reverse transcription quantitative PCR (RT-qPCR)

A total of six PDAC tissues and six control tissue specimens were collected. The six pairs of control and PDAC samples were obtained from the clinic at the First Affiliated Hospital of Nanchang University. The study was approved by the Ethics Committee with Ethical Number 2022 (4–027). RNA concentration was detected using the NanoPhotometer N50, and reverse transcription was performed using a cDNA synthesis kit (HP All-in-one qRT Master Mix II RT203-Ver.1), and primers were synthesized by a biological company (Table 1). The RT-qPCR experiments were performed with GAPDH as the internal reference gene, and the expression levels of the prognostic genes were calculated using the 2−ΔΔCT method. “Graphpad Prism” (v 10.1.2) (Al-Rawi et al., 2023) was employed to plot and calculate the P value. Differences between PCR experimental categories were obtained through a t-test (P < 0.05).

PrimerSequenceHMGA1 FGCA​TCC​GCA​TTT​GCT​ACC​AGHMGA1 RTCT​CAG​TGC​CGT​CCT​TTT​CCKIF2C FTCC​GTG​TCA​GCC​ATC​AAG​AGKIF2C RCAG​GCA​AAC​AGT​CGG​GTA​CTInternal reference-GAPDH FATG​GGC​AGC​CGT​TAG​GAA​AGInternal reference-GAPDH RAGG​AAA​AGC​ATC​ACC​CGG​AG2.11 Statistical analysis

Bioinformatics analyses were conducted using the R programming language (v 4.2.2). The disparities between the two groups were evaluated via the Wilcoxon test, with a significance level of P < 0.05 considered statistically significant.

3 Results3.1 There were five candidate genes affecting the development of PDAC in different ways

RGs were intersected, and five intersection genes were obtained as candidate gene. There were 4,431 DEGs between the tumor and control samples in the training set, including 2,363 upregulated genes and 2,068 downregulated genes (Figure 1A). Furthermore, to identify DEGs associated with lactylation and programmed cell death, 4,431 DEGs, 332 LRGs, and 1,548 PCDenes (HMGA1, KIF2C, GAPDH, HDAC1, and SOD1) (Figure 1B). Subsequently, GSEA was carried out for candidate genes. Functional enrichment analysis revealed that HMGA1 (87 pathways), KIF2C (102), GAPDH (101), HDAC1 (93), and SOD1 (114) are involved in diverse biological processes (Figures 1C–G; Supplementary Tables S3–S7). Notably, HMGA1 and KIF2C shared enrichment in 69 pathways, and several lactate metabolism–related pathways, such as KEGG_GLYCOLYSIS_GLUCONEOGENESIS and KEGG_GLYCEROPHOSPHOLIPID_METABOLISM, were significantly upregulated. To clarify the biological relevance of GAPDH and SOD1 in the lactylation-associated context, we found that GAPDH showed a strong positive correlation with LDHA (r = 0.73), the key enzyme for lactate production (Supplementary Figure S1A) Furthermore, single-gene GSEA confirmed that both GAPDH and SOD1 are significantly enriched in lactylation-related metabolic pathways (Supplementary Figure S1B). Overall, the comprehensive analysis has offered a wealth of information that could guide further research efforts aimed at deciphering the underlying mechanisms of PDAC and devising innovative therapeutic approaches.

Identification of candidate genes and functional enrichment analysis. (A) Volcano plot showing differentially expressed genes (DEGs) between PDAC tumor samples and control samples in the training set. Red dots represent upregulated genes, blue dots represent downregulated genes, and gray dots represent genes with no significant difference. (B) Venn diagram showing the intersection of DEGs, lactylation-related genes (LRGs), and programmed cell death-related genes (PCDRGs). (C–G) Gene Set Enrichment Analysis (GSEA) results for the five candidate genes. The top five enriched KEGG pathways are shown for HMGA1 (C), KIF2C (D), GAPDH (E), HDAC1 (F), and SOD1 (G).

3.2 A total of two genes that impact PDAC prognosis were discovered

Univariate Cox regression analysis was executed based on five candidate genes. Altogether, two genes were obtained (HR ≠ 1, P < 0.2) (Figure 2A) (Emura et al., 2019). Next, these two genes (HMGA1 and KIF2C) successfully satisfied the criteria of the PH hypothesis test (P > 0.05) (Figure 2B). Subsequently, the RSF model was constructed based on the above two prognostic genes to obtain the risk scores for each patient. It could be seen from the residual plot that when ntree was 29, the error rate of the model was the lowest; therefore, the parameters of the model were set as ntree = 29 and maximum split node mtry = 5 (Figure 2C).

Identification of prognostic genes and construction of a Random Survival Forest model. (A) Forest plot showing univariate Cox regression analysis results for the five candidate genes. (B) Proportional hazards (PH) assumption test results for HMGA1 and KIF2C. (C) Random Survival Forest (RSF) model optimization plot showing the relationship between the number of trees (ntree) and prediction error rate.

3.3 A reliable prognostic risk model was established

Within the tumor samples with survival information in the training set, patients were segmented into HRG (n = 60) and LRG (n = 82) based on the optimal threshold value (40.52855). Then it was found that there were substantial disparities in survival probability between HRG and LRG (P < 0.0001) (Figure 3A); as risk scores increased, survival time decreased, and more deaths occurred (Figure 3B). The AUC value (1-(0.741), 2-(0.726), and 3-(0.758) years) indicated that this prognostic risk model functioned effectively in predicting the survival status (Figure 3C). Then, in the validation set, patients were segmented into HRG (n = 21) and LRG (n = 44) in accordance with the optimal threshold value (43.94293). There were also substantial disparities in survival probability between HRG and LRG (P = 0.00013) (Figures 3D,E). The AUC value (1-(0.723), 2-(0.746), and 3-(0.638) years) also indicated it performed well in predicting the survival status (Figure 3F). The results as mentioned above demonstrated that the prognostic model exhibited remarkable robustness.

Construction and validation of the prognostic risk model. (A) Kaplan-Meier survival curves comparing the high-risk group (HRG) and low-risk group (LRG) in the training set. (B) Risk score distribution (upper panel), survival time distribution (middle panel), and survival status (lower panel) of patients in the training set. Patients are ordered by increasing risk score. As risk scores increased, survival time decreased, and mortality rate increased. (C) Time-dependent receiver operating characteristic (ROC) curves for the prognostic model in the training set. The area under the curve (AUC) values for 1-, 2-, and 3-year overall survival were 0.741, 0.726, and 0.758, respectively, indicating good predictive performance. (D–F) Validation of the prognostic model in the GSE62452 validation set. Kaplan-Meier curves (D), risk score and survival status distribution (E), and time-dependent ROC curves (F) confirmed the robustness of the model (P = 0.00013; AUC: 0.723, 0.746, and 0.638 for 1-, 2-, and 3-year survival). (G) Heatmap showing the distribution of clinicopathological features and prognostic gene expression across risk groups. HMGA1 and KIF2C showed higher expression in the HRG compared to the LRG.

The distribution of clinicopathological features in the two risk groups was demonstrated (Figure 3G). The results of the heatmap showed that the prognostic genes were lowly expressed in the LRG and highly expressed in the HRG; most patients were in Stage II, Stage T3/T4, and N1; the majority of patients were over 40 years old, while there was no obvious distribution pattern in terms of gender.

3.4 Differences in enrichment pathways and immune microenvironment were strongly associated with different prognostic populations

The GSEA results based on the KEGG gene set, included the following pathways: Cell cycle, Neuroactive ligand receptor interaction, Chemokine signaling pathway, Cytokine-cytokine receptor interaction, and DNA replication (Figure 4A; Supplementary Table S8). The upregulated pathways of GSVA included DNA repair, MYC targets V2, and G2M checkpoint, while the downregulated pathways included Myogenesis, KRAS signaling DN, and Angiogenesis (Figure 4B; Supplementary Table S9). These findings furnished a robust basis for the ongoing investigation of the molecular mechanisms underlying PDAC.

Functional enrichment analysis and immune microenvironment characteristics. (A) Gene Set Enrichment Analysis (GSEA) based on KEGG gene sets comparing HRG and LRG. (B) Gene Set Variation Analysis (GSVA) showing differentially enriched Hallmark pathways between HRG and LRG. (C) Heatmap displaying the infiltration levels of 28 immune cell types in HRG and LRG samples. Color intensity represents the ssGSEA enrichment score for each immune cell type. (D) Box plots comparing the infiltration levels of 22 significantly different immune cell types between HRG and LRG. Most immune cells showed lower infiltration scores in HRG (**P < 0.01, Wilcoxon test). (E) Correlation heatmap showing the relationships among differentially infiltrated immune cells. Immature B cells and activated B cells showed the strongest positive correlation (cor = 0.97, P < 0.001), while neutrophils and activated B cells showed the strongest negative correlation (cor = −0.68, P < 0.001). (F) Correlation analysis between prognostic genes (HMGA1 and KIF2C) and immune cells. Both genes showed the strongest negative correlation with mast cells (cor = −0.36 and −0.53 for KIF2C and HMGA1, respectively, P < 0.01). (G) Violin plot comparing ESTIMATE scores (stromal score, immune score, and ESTIMATE score) between HRG and LRG. All three scores were significantly lower in HRG (***P < 0.001, Wilcoxon test). (H) Correlation analysis between the risk score and tumor microenvironment scores (StromalScore, ImmuneScore, and ESTIMATEScore). Red indicates positive correlation, and blue indicates negative correlation. ***P < 0.001.

Figure 4C illustrates the infiltration levels of immune cells within the HRG and the LRG. The outcomes showed that 22 types of immune cells had substantial disparities between the HRG and the LRG. The outcomes revealed that most cells had lower scores in the HRG (P < 0.01) (Figure 4D). Among them, immature B cells and activated B cells demonstrated the strongest positive correlation (cor = 0.97, P < 0.001), while neutrophils and activated B cells demonstrated the strongest negative correlation (cor = −0.68, P < 0.001) (Figure 4E). Furthermore, KIF2C and HMGA1 demonstrated the strongest negative relationship with mast cells (cor = −0.36, −0.53, P < 0.01) (Figure 4F; Supplementary Table S10). Furthermore, there was a notable disparity in the scores of the immune microenvironment in two risk groups (P < 0.001) (Figure 4G). Moreover, there was a significant negative correlation between the risk score and the ESTIMATE score, immune score, and stromal score (cor = −0.35, −0.32, −0.37, P < 0.001) (Figure 4H). The above outcomes suggested that prognostic genes might influence the immune cell infiltration of PDAC, and this might provide a reference for the clinical management of PDAC.

The lactylation-associated score was significantly higher in the low-risk group than in the high-risk group, indicating a distinct lactylation-related metabolic context between risk subgroups (Supplementary Figure S2A). Moreover, correlation analysis showed that HMGA1 expression was negatively correlated with the lactylation-associated score (r = −0.38, P < 0.001) (Supplementary Figure S2B).

3.5 The patients' sensitivity to drugs and gene mutation landscape affected the treatment of PDAC patients

The lower IC50 values indicated that the drugs could achieve a significant therapeutic effect at a lower dose and reduce toxicity and side effects. A total of 127 different drugs were identified (P < 0.05) (Supplementary Table S11). Among the top 20 drugs with significant differences, Sapitinib and Lapatinib were more effective in the HRG, while AZ6102 and Doramapimod were more effective in the LRG (P < 0.0001) (Figure 5A).

Drug sensitivity analysis and mutation landscape. (A) Box plots showing the top 20 drugs with significant differences in predicted IC50 values between HRG and LRG. Lower IC50 values indicate higher drug sensitivity. Sapitinib and Lapatinib were more effective in HRG, while AZ6102 and Doramapimod were more effective in LRG (****P < 0.0001, Wilcoxon test). (B,C) Waterfall plots displaying the top 20 genes with the highest mutation frequencies in HRG (B) and LRG (C). (D) Violin plot comparing tumor mutational burden (TMB) scores between HRG and LRG. TMB was significantly higher in HRG (***P < 0.001, Wilcoxon test). (E) Scatter plot showing the positive correlation between risk scores and TMB scores (cor = 0.35, P = 4.5e-05, Spearman correlation).

The top 20 genes exhibiting the highest mutation frequencies in the HRG and LRG were presented (Figures 5B,C). The obtained results clearly demonstrated that the missense mutations emerged as the most frequently occurring mutation type. Notably, the mutation rates of KRAS and TP53 not only surpassed 50% in both the HRG and LRG, but also stood out as the most common mutations. This indicated their significant role and potential influence in the genetic alterations associated with these groups, potentially affecting various biological processes and disease progression. Furthermore, there was a marked discrepancy in TMB was observed between the HRG and LRG (P < 0.001), and TMB exhibited a notable positive relationship with risk scores (cor = 0.35, P = 4.5e-05) (Figures 5D,E). This provided a strong basis for clinicians to formulate personalized treatment plans for patients with different risk stratifications.

3.6 The prognostic genes located on autosomes showed significant differences in expression between tumor samples and control samples

The chromosomal localization map revealed that two prognostic genes were located on autosomes: KIF2C on chromosome one and HMGA1 on chromosome 6 (Figure 6A). Both the immunohistochemical analysis of prognostic genes and the expression analysis of prognostic genes in the training set showed that HMGA1 and KIF2C were upregulated in PDAC patients (P < 0.05) (Figures 6B,C). The RT-qPCR results showed that HMGA1 and KIF2C were still prominently overexpressed in PDAC samples (P < 0.01) (Figure 6D). This was in accordance with the outcomes of the bioinformatics analysis, indicating that the results of the bioinformatics analysis were reliable.

Chromosomal localization and expression validation of prognostic genes in PDAC. (A) Chromosomal localization map showing the distribution of prognostic genes on autosomes. KIF2C is located on chromosome 1, and HMGA1 is located on chromosome 6. (B) Immunohistochemical analysis of prognostic genes. (C) Expression analysis in the training set demonstrating upregulated expression of HMGA1 and KIF2C in PDAC patients compared to control samples. (D) RT-qPCR validation showed significant overexpression of HMGA1 and KIF2C in PDAC samples compared with normal controls (P < 0.01), confirming the reliability of the bioinformatics analysis results.

4 Discussion

PDAC remains one of the most lethal malignancies, typically characterized by an immunosuppressive tumor microenvironment, in which lactate-driven lactylation modification can promote immune evasion by regulating myeloid-derived suppressor cells (Saloman et al., 2016; Peng et al., 2024). Recent studies indicate that programmed cell death pathways are closely linked to metabolic reprogramming in PDAC (Chen et al., 2021; Zhang et al., 2025). Still, the mechanistic interactions between lactylation and programmed cell death remain urgently in need of exploration. This study integrates multi-dimensional bioinformatics analyses of TCGA and GEO databases to systematically identify two prognostic genes, HMGA1 and KIF2C, closely associated with lactylation and programmed cell death. It constructs a robust risk-stratification model and elucidates their roles in immune microenvironment remodeling, functional pathway regulation, and drug sensitivity, thereby providing new insights for precision medicine in PDAC.

This study identified HMGA1 and KIF2C as genes significantly upregulated in PDAC with crucial prognostic significance. HMGA1 encodes a non-histone chromosomal structural protein belonging to the high-mobility group protein A family, which is overexpressed in most malignancies and associated with tumor invasiveness, chemoresistance, and therapeutic response (Wang L. et al., 2022). Functionally, HMGA1 coordinates oncogenic programs through multiple mechanisms: it promotes lipid synthesis in colorectal cancer (Zhao et al., 2024), enhances the efficacy of palbociclib via PI3K/mTOR signaling in intrahepatic cholangiocarcinoma (Li Z. et al., 2023), and upregulates the pentose phosphate pathway in esophageal squamous cell carcinoma (Liu et al., 2024). Crucially, in PDAC, HMGA1 directly binds the FGF19 promoter and recruits his