Exploring prognostic genes related to lactylation and programmed cell death in pancreatic ductal adenocarcinoma: a comprehensive study combining bulk transcriptomics and experimental verification

TL;DR

This study identifies two genes linked to lactylation and programmed cell death in pancreatic cancer, offering a new risk model for better treatment strategies.

Contribution

The study introduces a novel risk model based on lactylation and programmed cell death-related genes in pancreatic cancer.

Findings

HMGA1 and KIF2C are identified as prognostic genes associated with lactylation and programmed cell death in PDAC.

A risk model was developed showing strong correlations with immune cell infiltration and drug sensitivity differences between risk groups.

Both genes were upregulated in PDAC patients and confirmed in clinical samples.

Abstract

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.

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.

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.

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.

Infographic summarizing a multi-step analysis pipeline for pancreatic cancer that includes candidate gene identification using PDAC versus control comparison and Venn diagram, prognosis risk modeling with gene selection and survival curves, and mechanism analysis involving gene set enrichment analysis, immune infiltration heatmap, and drug sensitivity box plots.

Figures

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).

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.

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.

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.

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).

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.

Tables

TABLE 1. Primer Sequence.

Primer	Sequence
HMGA1 F	GCA​TCC​GCA​TTT​GCT​ACC​AG
HMGA1 R	TCT​CAG​TGC​CGT​CCT​TTT​CC
KIF2C F	TCC​GTG​TCA​GCC​ATC​AAG​AG
KIF2C R	CAG​GCA​AAC​AGT​CGG​GTA​CT
Internal reference-GAPDH F	ATG​GGC​AGC​CGT​TAG​GAA​AG
Internal reference-GAPDH R	AGG​AAA​AGC​ATC​ACC​CGG​AG

Full text

Introduction

1

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.

Materials and methods

2

Data collection

2.1

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).

Differential expression analysis

2.2

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 |log_2_ 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.

Identification and functional characterization of candidate genes

2.3

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).

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

2.4

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.

Clinicopathological characteristics analysis

2.5

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.

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

2.6

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 log_2_FC values. Subsequently, the genes were arranged in descending order by log_2_ fold change log_2_FC. 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.

Analysis of immune microenvironment

2.7

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).

Drug susceptibility prediction and gene mutation landscape analysis

2.8

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 (IC_50_) 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).

Prognostic gene analysis

2.9

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).

Reverse transcription quantitative PCR (RT-qPCR)

2.10

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).

Statistical analysis

2.11

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.

Results

3

There were five candidate genes affecting the development of PDAC in different ways

3.1

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).

A total of two genes that impact PDAC prognosis were discovered

3.2

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.

A reliable prognostic risk model was established

3.3

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.

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

3.4

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).

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

3.5

The lower IC_50_ 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.

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

3.6

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.

Discussion

4

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 histone-activating modifiers, thereby promoting FGF19 expression, which, in turn, drives autocrine and paracrine tumor growth (Chia et al., 2023). The latest evidence further indicates that elevated HMGA1 is associated with poor prognosis and plays an immunosuppressive role in the tumor microenvironment (Song et al., 2025), and single-cell analysis identifies HMGA1 as a marker of a subpopulation of metastatic cancer cells (Yang D. et al., 2023). Our study expands this understanding by demonstrating that HMGA1 is co-expressed with lactylation-related genes and negatively correlated with mast cell infiltration, suggesting its potential role in lactylation-mediated immune evasion. The consistency between our bioinformatics predictions and experimental verification results highlights the robustness of the findings and positions HMGA1 as a promising therapeutic target.

Similarly, KIF2C, a member of the kinesin superfamily, encodes a key regulator of spindle dynamics and chromosome segregation during mitosis, playing an essential role in maintaining genomic stability (Kreis et al., 2024). The dysregulation of KIF2C is closely associated with tumor progression; its elevated expression can promote tumor cell migration, invasion, and chemoresistance, and impair DNA damage repair (Li R.-Q. et al., 2024). In clear cell renal cell carcinoma, KIF2C exacerbates tumor grade, stage, metastasis status, and a poor prognosis by activating the JAK2/STAT3 signaling axis (Deng et al., 2023). Importantly, KIF2C is significantly overexpressed in PDAC tissues and cell lines (ASPC-1, MIA-PACA-2), driving tumor growth and metastasis (Huang et al., 2023), and its expression level is negatively correlated with patient survival rate (Xiong et al., 2022). This study validated these findings, confirming that KIF2C expression is low in normal pancreatic exocrine cells and undetectable in endocrine cells, but significantly upregulated in PDAC. Notably, KIF2C expression is strongly negatively correlated with mast cell infiltration and associated with the enrichment of cell cycle and G2/M checkpoint pathways in high-risk patients, consistent with its established roles in mitotic regulation and genomic instability. These findings indicate that KIF2C is a key driver of PDAC invasiveness and a potential biomarker for risk stratification.

Functional enrichment analysis revealed that high-risk patients were characterized by abnormal activation of cell cycle and G2/M checkpoint pathways, which was not only associated with lower recurrence-free survival and overall survival, but also predicted enhanced sensitivity to adjuvant chemotherapy and potential resistance to immune checkpoint inhibitors (Mao et al., 2022). This finding is consistent with the established roles of HMGA1 and KIF2C in promoting mitotic progression and genomic instability. In addition to cell cycle dysregulation, our GSEA results also highlighted the enrichment of neuroactive ligand-receptor interaction pathways, reflecting the critical role of neural infiltration in PDAC metastasis. Tumor cells promote distant organ dissemination by utilizing neural pathways (Guo et al., 2023; Bhattacharjee and Ghosh, 2025). Furthermore, upregulation of G2/M checkpoint components (including PLK1) is associated with MEK inhibitor resistance (Zhao et al., 2025), suggesting a potential therapeutic vulnerability in high-risk patients. Notably, PDAC progression is also closely related to angiogenesis; tumor-associated macrophages (TAMs) can promote increased vascular density (Yang et al., 2021). These pathway alterations collectively reveal a synergistic oncogenic program driven by HMGA1 and KIF2C, integrating metabolic reprogramming, cell-cycle dysregulation, and immune-escape mechanisms.

The PDAC tumor microenvironment plays a critical role in disease progression, and the pattern of immune cell infiltration significantly impacts clinical outcomes (Werba et al., 2023). Our analysis revealed 22 types of differentially infiltrated immune cells, among which the immune cell infiltration score in the high-risk group was significantly lower. High B cell infiltration is associated with improved PDAC survival (Delvecchio et al., 2022); we observed a strong positive correlation between immature and activated B cells, suggesting a coordinated adaptive immune response. Conversely, N2 polarized neutrophils, as the dominant population in the PDAC microenvironment, drive tumor progression by promoting proliferation, angiogenesis, tissue remodeling, and immune suppression (Mantovani et al., 2021). Neutrophils undergo a transition from TAN-3 to TAN-0 phenotype during transendothelial migration into the tumor (Wang L. et al., 2023) and enhance metastatic ability by secreting TNF-α and TGF-β1 (Lianyuan et al., 2020). We observed a strong negative correlation between neutrophils and activated B cells, reflecting the antagonistic effect between innate and adaptive immunity in PDAC. Notably, both HMGA1 and KIF2C showed the strongest negative correlation with mast cell infiltration. Mast cells in the PDAC tumor microenvironment have been confirmed to promote immune suppression and fibrous stromal hyperplasia (Ma et al., 2024). In melanoma, mast cells hinder tumor cell homing and metastasis by inhibiting HMGA1 secretion from tumor cells (Benito-Martin et al., 2023), suggesting a potential regulatory axis. Our study proposes that tumors with high HMGA1/KIF2C expression may be associated with reduced mast cell infiltration or altered mast cell functional states, potentially influencing immune surveillance. This requires further verification through co-culture systems, flow cytometry, or spatial transcriptomics to elucidate the causal relationship and therapeutic implications.

In summary, this study successfully identified HMGA1 and KIF2C as prognostic genes and established a robust risk-stratification model integrating lactylation and programmed cell death. Through multi-dimensional analysis, the roles of these two genes in cell cycle regulation and immune microenvironment remodeling were elucidated, providing new insights into the mechanism of PDAC progression. However, several limitations still need attention. First, the conclusions of this study are based on retrospective analysis of public databases, lacking validation from prospective cohorts or functional experiments. Subsequent studies should combine in vitro and in vivo models to dissect the causal relationship of HMGA1 and KIF2C in lactylation-mediated tumorigenesis and programmed cell death evasion. Second, although the bioinformatics predictions are consistent with the published literature, the causal relationship between HMGA1/KIF2C and specific immune cell subsets (such as mast cells and neutrophils) still needs verification using flow cytometry, co-culture systems, or spatial transcriptomics. Furthermore, the molecular bridge connecting lactylation modification and programmed cell death pathways remains unclear. Integrating proteomics with CRISPR-based functional screening will help clarify the lactylation substrates and programmed cell death regulatory factors regulated by HMGA1 and KIF2C. Finally, the clinical application value of this model needs validation in independent multicenter cohorts, and its predictive value for immune therapy and targeted therapy responses awaits prospective evaluation. At the same time, the number of clinical samples used for qRT-PCR validation in this study was relatively limited. Future studies should expand the sample amount to further improve statistical power and verify the robustness of the results. It should be noted that this study primarily relied on transcriptomic data analysis and has not yet obtained proteomic evidence, such as identification of specific lactation sites or quantitative detection of lactation levels, to directly verify the lactation modification of HMGA1 and KIF2C. Therefore, future studies should combine proteomic analysis or related experiments to further confirm the direct lactation mechanism.

Conclusion

5

In summary, this study not only systematically integrated the lactylation and programmed cell death-associated genes HMGA1 and KIF2C into a PDAC prognostic evaluation system for the first time but also revealed their potential immune regulatory and pathway interaction mechanisms. These findings provide new insights into the precise typing, prognostic assessment, and development of targeted immune combination therapy strategies for PDAC.