The study flowchart is presented in Fig. 2. The present study utilized the RNA-sequencing (RNA-seq) dataset of PAAD from the TCGA database (https://portal.gdc.cancer.gov/), along with its corresponding mutation dataset (available on the TCGA official website as of June 30, 2025). The RNA-seq dataset included a total of 181 samples, comprising 4 non-tumor (11A) samples and 177 tumor (01A) samples, which were used for prognostic gene exploration and model construction. Additionally, the mutation dataset contained 180 samples for mutation landscape analysis. From Gene Expression Omnibus, the GSE28735 dataset (11A=45, 01A=45; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28735) was employed for differential expression analysis and validation of the predictive models constructed using the TCGA dataset.

Differential expression analysis of mRNAs was performed using the limma package in R software (https://bioconductor.org/packages/limma). Genes with an adjusted P<0.05 and a log₂(fold change) >1 or <-1 were considered significantly differentially expressed. The identified differentially expressed genes (DEGs) were subsequently visualized using volcano and heat maps. Functional enrichment analyses were systematically performed on both upregulated and downregulated gene sets using Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) pathway analyses.

A univariate Cox analysis was first performed to identify genes associated with prognosis (P<0.05). LASSO regression was then applied to refine the selection of key genes influencing patient outcomes. A prognostic model was developed using the training cohort (TCGA dataset), incorporating these genes and their corresponding coefficients. Patients were stratified into high- and low-risk groups based on the median risk score. Survival analyses were performed for both groups, and the predictive accuracy of the model was evaluated.

Nomograms were constructed to visualize the results of Cox regression. A scoring system was developed based on the magnitude of the regression coefficients for all independent variables, allowing for the calculation of a total score for each patient to estimate the probability of clinical outcomes. The nomograms were used to predict the 1-, 3- and 5-year mortality rates in patients (TCGA dataset). Receiver operating characteristic (ROC) curves were then generated to assess the predictive performance of the model.

The predictive model was further validated using the GSE28735 cohort. Patients were categorized into high- and low-risk groups based on the median risk score. Survival analyses were performed, and the accuracy of the model and its value as an independent prognostic indicator were assessed.

Immune infiltration and immune function analyses were performed using the ssGSEA algorithm (https://bioconductor.org/packages/GSVA/). The ‘ESTIMATE’ software package (https://bioconductor.org/packages/estimate/) was used to compute immune and stromal scores, with each sample assigned a StromalScore, ImmuneScore and ESTIMATEScore. Box plots were generated to compare immune scores between high- and low-risk groups. To evaluate the potential efficacy of immunotherapy, the Tumor Immune Dysfunction and Exclusion (TIDE) score for each sample was calculated using the TIDE database (http://tide.dfci.harvard.edu/). A scatter plot illustrating the correlation between TIDE scores and risk groups was generated to assess their relationship.

Mutation data for pancreatic cancer were obtained from the TCGA database using the ‘MAfTools’ package (https://bioconductor.org/packages/maftools), with ‘Mutect2’ selected as the mutation data type. Tumor mutation burden (TMB) for each sample was separately calculated for high- and low-risk groups using the TMB function within the ‘MAfTools’ package.

The present study used paired pancreatic ductal models for primary screening: PANC-1 cells (representing classical PDAC) and HPDE-6 cells, an immortalized normal pancreatic ductal epithelial line. PANC-1 was selected due to its extensive use in PDAC research (>4,400 PubMed-indexed studies as of June 4, 2025) and its genomic profile, which harbors key driver mutations (KRAS and TP53) frequently identified in patient-derived PAAD cohorts. This mutational status was specifically documented for PANC-1 by Deer et al (43). The PANC-1 cell line (HLCL-023) and the HPDE-6 cell line (H1-3201) were purchased from iCell Bioscience, Inc., and OriCell Therapeutics. PANC-1 was cultured in DMEM (Provider: Hyclone,) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. HPDE-6 was maintained in the iCell Primary Epithelial Cell Culture System, supplemented with 5 ml Primary Epithelial Cell Culture Additive (iCell Bioscience, Inc.), 2% FBS and 5 ml penicillin/streptomycin. HPDE-6 cells and PANC-1 cells were removed from the liquid nitrogen tank and placed in a 37°C water bath with rapid shaking to thaw. After thawing, the cells were transferred to centrifuge tubes, and 2 ml of corresponding complete medium was added to each of them. The supernatant was removed by centrifugation, 2 ml of complete medium was re-added, gently pipetted and mixed, and the content was transferred to a six-well plate for incubation in a 37°C incubator. Cell passaging was performed when the cell density grew to 80–90%. The passaging ratio was 1:3. When the cell density of the passage reaches ~95%, reverse transcription-quantitative PCR (RT-qPCR) experiments were performed on each cell.

Total RNA was extracted from HPDE-6 cells and PANC-1 cells using Trizol Plus (catalog number 15596026; Takara Bio, Inc.; iCell Bioscience, Inc.), and cDNA (reverse transcribed at 42°C for 15 min, then at 85°C for 5 min, and then placed on ice to cool) was synthesized using a first-strand cDNA synthesis kit (catalog number #K1622; Thero Fisher Scientific, Inc.). The RT-qPCR reaction mixture was prepared according to the standard protocol, consisting of 8.2 µl ddH₂O, 10.0 µl SYBR Green Master Mix (catalog number 04913914001; Roche Diagnostics GmbH), 0.4 µl forward primer, 0.4 µl reverse primer and 1.0 µl cDNA template, making a final volume of 20.0 µl. The components were thoroughly mixed in a PCR reaction tube. Primers for RT-qPCR were designed and synthesized using Primer Premier 5.0 software (PREMIER Biosoft) and are listed in Table SI. β-actin was used as the internal reference. Amplification was performed under the following conditions: Initial denaturation at 95°C for 5 min, followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec. The instrument automatically analyzed the results. After PCR amplification, the qPCR system was calibrated using the negative control to set the threshold and baseline, determining the Ct value for each sample. The validity of the Ct value was confirmed through the melting curve analysis. The results were then exported, and differences in target gene expression between the control group and each experimental group were analyzed using the 2^−ΔΔCq method (44).

Data were analyzed using R software (version 4.3.1; The R Foundation). The following packages were used for analysis: The limma package was used for variance analysis (Linear Models and Empirical Bayes); the survival package (https://cran.r-project.org/package=survival) was used for Kaplan-Meier curve plotting (survfit function) and cox regression analysis (coxph function); the GLMNet package (https://cran.r-project.org/package=glmnet) was used for regularization analysis (GLMNet function) in LASSO regression; pROC (https://cran.r-project.org/package=pROC) and timeROC (https://cran.r-project.org/package=timeROC) packages were used for ROC analysis and timeROC analysis (ROC function and timeROC function); the RMS package (https://cran.r-project.org/package=rms) was used for nomogram analysis and calibration curve analysis (CPH function and calibrate function); the GSVA package (https://bioconductor.org/packages/GSVA) for single sample Gene Set Enrichment Analysis immune infiltration and immune function analysis (gsva function); the estimate package was used to calculate tumor purity, the immune score and the stromal score (estimateScore function); and the MafTools package was used to construct a mutation landscape analysis (TMB function). Statistical analyses for validation experiments were performed using GraphPad Prism 9 software (Dotmatics). All quantitative data are expressed as mean ± standard error of the mean, based on results from ≥3 independent experiments. P<0.05 was considered to indicate a statistically significant difference.

Based on the criteria of log₂(foldchange) >1 or <-1, DEGs were identified by comparing PAAD tissues with normal tissues. A total of 244 upregulated and 170 downregulated DEGs were detected. Volcano and heat maps were then generated to visualize these findings (Fig. 3A and B).

Subsequent GO and KEGG enrichment analyses were performed on both upregulated and downregulated DEGs to assess differences in functional pathways. Upregulated DEGs demonstrated significant enrichment in pathways characteristic of pancreatic cancer, especially those involving ECM remodeling. Regarding biological processes, these genes were primarily associated with ECM organization and extracellular structure organization. At the cellular component level, enrichment was observed in the collagen-containing ECM and basement membrane. Molecular function analysis indicated involvement in ECM structural constituent and integrin binding (Fig. 3C). Conversely, downregulated DEGs were associated with normal pancreatic physiological activities. Biological processes included digestion and the disruption of cellular anatomical structure in another organism. Cellular components were associated with platelet α granule and brush border membrane. Molecular functions involved serine-type peptidase activity and serine hydrolase activity (Fig. 3E). KEGG pathway analysis further revealed several key signaling pathways: Upregulated DEGs were enriched in ECM-receptor interaction and cytoskeletal regulation in muscle cells, whilst downregulated DEGs were enriched in pancreatic secretion and protein digestion and absorption pathways (Fig. 3D and F).

To further evaluate genes associated with prognosis, a univariate Cox analysis was performed on 414 DEGs in the TCGA dataset, identifying 18 genes with significant prognostic relevance. The results of the univariate Cox analysis for these 18 genes are presented in a forest plot (Fig. 4A). LASSO regression analysis was then performed (Fig. 4B and C), and multifactorial Cox regression analysis was used to calculate risk scores, resulting in a final prognostic model comprising 12 genes. Patients were categorized into high- and low-risk groups based on the median risk score. The results revealed that patients in the high-risk group had significantly worse prognoses than those in the low-risk group (P<0.05; Fig. 4D). Moreover, the ROC curves demonstrated that the area under the curve (AUC) values for predicting patient prognosis at 1, 3, and 5 years were 0.812, 0.832 and 0.858, respectively. These values exceeded those of clinical characteristics such as sex, age and disease stage (Fig. 4E and F).

Given the association between risk score and patient prognosis, a Cox regression analysis was performed to assess the relationship between the risk score and other common clinical parameters (Fig. 5A). A nomogram was then constructed to estimate the 1-, 3- and 5-year overall survival (OS) of patients with PAAD (Fig. 5B). Additionally, the AUC values for several clinical factors used to predict 1-, 3- and 5-year OS were assessed. The results revealed that the nomogram exhibited good prognostic predictive ability, though its performance was not as strong as the riskScore clinical feature (Fig. 5C). Furthermore, calibration curve analysis demonstrated that the predicted values of the nomogram closely aligned with actual observed survival outcomes, indicating high accuracy (Fig. 5D).

Similar results were observed in the external validation cohort GSE28735. In this cohort, patients in the high-risk group had significantly worse prognoses than those in the low-risk group (Fig. 6A). Due to the absence of survival data beyond 3 years, ROC curves were analyzed for risk scores at 1, 2 and 3 years. The AUC for these time points were 0.845, 0.921 and 0.960, respectively (Fig. 6B). Univariate and multivariate Cox analyses confirmed that the risk score functioned as an independent prognostic factor (Fig. 6C).

The immune microenvironment serves a crucial role in determining patient prognosis (45). In the present study, the low-risk group exhibited significantly elevated immune infiltration compared with the high-risk group in myeloid-derived suppressor cells (MDSCs), macrophages and regulatory T cells (Fig. 7A). Additionally, the T-cell co-stimulation immune function pathway demonstrated significantly increased activity in the low-risk group compared with in the high-risk group (Fig. 7B). Furthermore, immune-related scores, namely the ESTIMATE, immune and stromal scores, were significantly higher in the high-risk group than in the low-risk group (Fig. 7C).

The present study also evaluated differences in immunotherapy responses between the two groups (Fig. 7D). The results indicated that TIDE scores were significantly lower in the low-risk group compared with in the high-risk group, suggesting a lower likelihood of tumor immune escape and a greater potential benefit from immunotherapy. Moreover, TIDE analysis revealed a significant positive correlation between the TIDE score and the risk score (Fig. 7E). As the risk score increased, the TIDE score increased, indicating a higher likelihood of tumor immune escape and a reduced response to immunotherapy. These findings suggest that the risk score may serve as a predictive biomarker for ICI efficacy, aiding in the clinical selection of patients who are most likely to benefit from immunotherapy whilst minimizing overtreatment in non-responsive individuals.

Gene mutations serve a notable role in tumor progression and patient outcomes (46). In the present study, the five most frequently mutated genes in both the high- and low-risk groups were KRAS, TP53, ‘mothers against decapentaplegic homolog 4’, cyclin dependent kinase inhibitor 2A and titin (TTN) (Fig. 8). These findings highlight the need for further investigation into targeted therapies against high-frequency mutations (such as KRAS and TP53) in conjunction with risk-stratification strategies to optimize precision treatment approaches.

To validate the aforementioned findings, RT-qPCR analysis was performed on the pancreatic cancer cell line, PANC-1, and the normal pancreatic cell line, HPDE-6. The results demonstrated significant differences in the expression levels of prognostic genes between the two cell lines (Fig. 9). Compared with the HPDE-6 group, the PANC-1 group demonstrated significantly upregulated expression of chromogranin B, dehydrogenase/reductase 9, endothelial cell specific molecule 1, fibronectin type III domain containing 1, glutaminase 2, 15-hydroxyprostaglandin dehydrogenase, matrix metallopeptidase 9 (MMP9), metallothionein 1M and TTN, whilst insulin growth factor-like family member 2, kynureninase and solute carrier family 16 member 10 were significantly downregulated.