To obtain the data needed for this study, we used Gene Expression Omnibus (GEO), a public functional genome database. We used two words “pancreatic cancer” and “radiation” as the keywords to obtain the required items, and finally selected the series GSE179351 and GSE225767 as the data set for our main analysis.

To obtain patient gene expression and clinical data, we included cases from The Cancer Genome Atlas (TCGA) project. We collected gene expression data from PAAD cases, including 179 tumor samples and 4 normal samples. All TCGA data, mRNA expression and clinical details were manually downloaded from the website and collated by R software.

GEO2R was used to identify differentially expressed genes in the geo data sets GSE179351 and GSE225767 after RT compared to before RT. The differentially expressed genes (DEGs) in PAAD samples in TCGA database were identified by R package “edgeR”. Using |log₂FC| ≥1 and p-value <0.05 as selection criteria, the range of DEGs was determined for further analysis. Volcano mapping using R software packages “limma” and “ggplot2”. Using gene annotation and analysis resources website “Metscape” (https://www.metascape.org/) for gene function analysis of enrichment.

After obtaining 121 co-up-regulated genes and 27 co-down-regulated genes between GSE179351 and GSE225767 as candidate genes for constructing radiation-resistant gene plates, we used univariate Cox regression to select the gene that was most correlated with OS (overall survival) in PAAD patients, and the p-value was set at 0.05.

Then, using the results of univariate Cox regression, we performed proportional risk regression analysis (multivariate Cox model) by the following formula.

Where, the risk score is the product of the mRNA expression of each key predictor gene (xi) and coefficient, which is derived from multivariate Cox regression analysis. Using the survival and survminer packages in R, 178 TCGA-PAAD patients were divided into high-risk group and low-risk group. On the basis of this grouping routine, survival differences between OS and risk scores were observed between the two groups. Logistic regression model was established by using R software package “glmnet”.

Random forest model was performed by R package ggRandomForests (v2.2.1). Regarding the allocation of the training and testing sets, 148 (83%) of samples to the training set and 30 (27%) to the testing set. In the parameter, there are 183 nodes and 454 edges contained inside this representation. Number of trees was set to 1000 and minimum size of terminal node was set to 10. Other parameters were default settings.

PANC-1 cells were purchased from Wuhan Promoter Life Science & Technology Co., Ltd. (China). All cell lines were cultured at 37°C in an atmosphere of 5% CO₂ using DMEM medium containing 10% fetal bovine serum (FBS) (Gibco) and 1% penicillin/streptomycin. When cell confluence reached 70%-80%, they were subjected to ionizing radiation treatment.

Trizol was used to extract RNA from cells. After RNA extraction, the quantity and quality of the RNA were analyzed using the Qnano spectrophotometry method from Yeasen. A reverse transcription reaction was performed using 1 μg of total RNA with Takara’s reverse transcription kit (RR047A). The resulting cDNA was subjected to qPCR analysis using Hieff UNICON^® Universal Blue qPCR SYBR Green Master Mix (Yeasen, Cat #11184ES08), and quantification was performed using Yeasen 80520ES03. In all experimental replicates, all expression levels were normalized to GAPDH. The primers used were as follows:

ADAMTS12 Forward (5′-CTTTGAAGGCGGCAACAGCAGA-3′)

ADAMTS12 Reverse (5′-TCTCACAGTCTGGCAGGAAGAG-3′)

AKR1C2 Forward (5′-CCGAAGCAAGATTGCAGATGGC-3′)

AKR1C2 Reverse (5′-TTTCAGTGACCTTTCCAAGGCTG-3′)

ATP8B2 Forward (5′-CGGCTATTCCTGCAAGATGCTG- 3′)

ATP8B2 Reverse (5′-GTCCTGATAGGTGAAGCCGTTG-3′)

CCN4 Forward (5′-AAGAGAGCCGCCTCTGCAACTT-3′)

CCN4 Reverse (5′-TCATGGATGCCTCTGGCTGGTA-3′)

CTHRC1 Forward (5′-CAGGACCTCTTCCCATTGAAGC-3′)

CTHRC1 Reverse (5′-GCAACATCCACTAATCCAGCACC-3′)

GREM1 Forward (5′-TCATCAACCGCTTCTGTTACGGC- 3′)

GREM1 Reverse (5′-CAGAAGGAGCAGGACTGAAAGG-3′)

P3H3 Forward (5′-CTGAGTGTCCTGCTCTTCTACC-3′)

P3H3 Reverse (5′-ATCGGAGGATGAAGCGCTGGAT-3′)

PAPPA Forward (5′-GGAACTGAAGAGAGTGAGCCATC-3′)

PAPPA Reverse (5′-CGTCGCATTGTTCACCTTGGTC-3′)

POSTN Forward (5′-CAGCAAACCACCTTCACGGATC- 3′)

POSTN Reverse (5′-TTAAGGAGGCGCTGAACCATGC-3′)

TAFA2 Forward (5′-GATCGGAAAGGATGGAGCTGTTC-3′)

TAFA2 Reverse (5′-GCGCATGTTCAATGTCATCAGCC-3′)

GAPDH Forward (5′-GGAGCGAGATCCCTCCAAAAT- 3′)

GAPDH Reverse (5′-GGCTGTTGTCATACTTCTCATGG-3′)

Single sample gene set enrichment analysis (ssGSEA) was performed using R package “GSVA” to calculate the infiltration level of 28 immune cells. The immune gene sets were sourced from Charoentong’s study (11). Finally, we compared the expression levels of immune checkpoint molecules (CD274 and CTLA4) between the two groups of patients. TIDE was used to predict immunotherapy.

In this study, we explored the predictive value of risk score for immunotherapy and chemotherapy efficacy. Tumor immune dysfunction and rejection (TIDE) score (http://tide.dfci.harvard.edu/) is a kind of used to evaluate tumor immune escape mechanism in the immune microenvironment of tools. TIDE scores predict a patient’s response to immunotherapy, such as immune checkpoint inhibitors, by taking into account the ability of tumor cells to escape immune and the role of tumor-related immunosuppression. The drug was estimated using the “oncopredict” package in R to predict chemotherapy drug sensitivity in each patient.

All statistical analysis and graphical visualizations were performed in R (version 4.3.2). Continuous variables were compared between groups using student t test or Wilcoxon rank sum test. P<0.05 was considered statistically significant (bilateral).

To identify genes that play a key role in RT resistance, we selected two datasets, GSE179351 and GSE225767, as study objects. In these two datasets, we obtained differential genes (DEGs) in cancer tissues after RT compared with those before RT through GEO2R online analysis ( Figures 1A, B ). In the GSE179351 data set, there were 594 up-regulated genes and 480 down-regulated genes in the cancer tissues after RT compared with those before RT ( Supplementary Table S1 ). In GSE225767 data, 1033 up-regulated genes and 637 down-regulated genes were found in cancer tissues after RT compared with those before RT ( Supplementary Table S2 ). In DEGs analysis, the cutoff value was set as |log₂FC| ≥1, and p-value <0.05. We screened 148 genes with identical expression changes in both datasets as candidate genes ( Figures 1C, D , Table 1 ). To clarify the function of these genes, we used the Metascape tool website to conduct a republic KEGG enrichment analysis of this group of genes, which are mainly involved in the regulation of cell behavior, the development and regeneration of tissues and organs. These changes may affect cancer cell survival, invasion, and response to treatment, providing important clues to understanding the mechanism of action of RT in PAAD ( Figure 1E ).

In order to investigate the effect of radiation-resistant genes on the prognosis of pancreatic cancer patients, we selected PAAD patients from the TCGA database and established univariate and multivariate Cox proportional risk regression models. Among the 148 candidate genes mentioned above, a total of 37 genes were considered to be related to patients’ overall survival (OS) by single-factor Cox regression model analysis. Of these, 26 genes are considered risk factors and 11 genes are considered protective factors. Therefore, we further included these 37 genes in multivariate Cox analysis to construct the genome associated with radiation resistance. According to the analysis results, 10 genes were screened out (p<0.05). Then, multivariate Cox regression coefficient and mRNA expression levels of key genes were used to establish a risk scoring formula: risk score = (-ADAMTS12*0.05286 + AKR1C2* 0.0163542-ATP8B2 *0.04237 + CCN4*0.032409 -CTHRC1*0.00454 + GREM1* 0.008157-P3H3 *0.03781 + PAPPA* 0.129872-POSTN * 0.00178-TAFA2 *1.80003). A hazard ratio greater than 1 indicates that patients with high gene expression are more likely to develop tumor progression after RT, while a hazard ratio less than 1 indicates that the gene is a protective factor ( Figure 2A ).

According to the risk score formula, 178 patients with PAAD were divided into high-risk and low-risk groups using the optimal cutoff value determined by the “surv_cutpoint” function ( Figure 2B ). The heat map identified the expression levels of 10 genes in a single patient, 100 in the high-risk group and 78 in the low-risk group ( Figure 2C , Table 2 ). As shown in Figure 2D , patients with higher risk scores after RT are more likely to develop tumor progression and have a relatively shorter median survival. Specifically, when we constructed survival curves based on the expression of a single key gene, only ATP8B2, GREM1, and TAFA2 genes could obtain statistically significant results (p< 0.05, Supplementary Figure S1 ). However, when patients were grouped according to risk scores, the high-risk group had significantly lower survival rates at 3 and 5 years, and even beyond, compared to the low-risk group (p = 0.001) ( Figure 2D ). This suggests that the prognosis model is effective in distinguishing between high-risk and low-risk patients, and that high-risk patients have poorer survival outcomes. To further verify the reliability, we established a multivariate logistic regression analysis to evaluate its effectiveness in predicting tumor progression. The validated AUC value of 0.7664 indicates that the model performs well in differentiating between samples with different risk of tumor progression ( Figure 2E ). In addition, we applied machine learning methods for patient risk prediction, but on our dataset, the machine learning model had a low AUC value ( Supplementary Figure S2 ). Using the random forest model, we further analyzed the contribution of the above 10 genes to the risk scoring model. The results showed that among the 148 gene candidates, TAFA2 contributed the most to the prediction model, followed by ATP8B2 ( Figure 2F ).

Next, we verified the 10 genes selected above to verify their relationship with radiation resistance in pancreatic cancer. The results showed that TAFA2 and POSTN were significantly elevated in IR-resistant pancreatic cancer cells ( Supplementary Figure S3 ), suggesting that TAFA2 may play an important role in radiotherapy resistance of pancreatic cancer.

To more fully validate the predictive power of the scoring model, we downloaded and analyzed three publicly available pancreatic cancer GEO datasets (GSE28735, GSE62452, and GSE57495). In each dataset, we calculated the AUC value via the ROC curve to assess the accuracy of the model’s prediction of patient risk. The results showed that all datasets had AUC values higher than 0.77, with the GSE28735 dataset having the highest AUC value of 0.8886, indicating that the scoring model has stable and high predictive performance across multiple datasets ( Figure 3A ), supporting its potential to be widely used in diverse pancreatic cancer patient populations.

Further, we explored the association between the risk score of the scoring model and the clinical characteristics of patients, aiming to analyze the clinical significance of the score. The study found that the risk score was not significantly associated with gender or pathological stage, suggesting consistent applicability of the score to patients of different genders and stages ( Supplementary Figure S4 ). However, there was a significant correlation between risk scores and patients’ survival status and age. Specifically, patients whose survival status was death had a significantly higher risk score than those who survived, suggesting that this score may be a powerful indicator of prognosis. In addition, patients older than 60 years had significantly higher risk scores than those younger than 60 years, a finding that may reflect a more aggressive or progressive course of disease in older patients ( Figure 3B ).

We also mapped the relationship between gender, high-low risk groups, and cancer stage to visualize the interaction patterns between these variables ( Figure 3C ). In Figure 3C , the distribution of patients of different genders in high and low risk groups and cancer stages is shown in the form of Sankey charts. Although there was no significant association between gender and risk score, we could observe differences in disease stage among patients in different risk groups. Such visualization not only helps to understand the relationship between variables, but also provides a reference for the development of further individualized treatment strategies.

In summary, this study confirmed the strong predictive ability of the scoring model in pancreatic cancer patients through external validation of multiple GEO datasets, and supported the clinical application potential of the model through correlation analysis with clinical characteristics. This validation method based on multiple data sets not only enhances the robustness of the model, but also lays a foundation for its popularization in clinical practice.

Based on the matched tumor RNA-seq data from PAAD patients, we identified 933 DEGs (p.adj < 0.05 and |log₂FC| ≥ 1) between the high-risk and low-risk groups, including 348 up-regulated genes and 585 down-regulated genes ( Figure 4A ). Next, we performed GO and KEGG enrichment analysis for these differential genes.

The results of GO enrichment analysis showed that the up-regulated genes were mainly enriched in GO terms associated with immune response, cell differentiation, and digestion, processes that may be involved in tumor development and changes in the immune microenvironment of PAAD ( Figure 4B ). GO enrichment results of down-regulated genes showed that these genes were mainly related to the tissue and structural components of the extracellular matrix, basic enzyme activity, and molecular binding activity, suggesting that tumor progression may be slower in low-risk patients, and tissue remodeling and signaling activities may be less active ( Figure 4D ).

The results of KEGG enrichment analysis showed that up-regulated genes were mainly enriched in pancreatic secretion, neuroactive ligand-receptor interactions, and protein digestion and absorption pathways ( Figure 4C ). The enrichment of these pathways suggests that patients at high risk of PAAD exhibit active biological characteristics in digestion, metabolism and nerve signaling, providing support for the growth, metabolic needs and microenvironment regulation of PAAD cells, thereby promoting the invasion and metastasis of cancer cells.

In contrast, KEGG enrichment of down-regulated genes showed that these genes were mainly concentrated in pathways such as neuroactive ligand-receptor interactions, cytoskeleton of muscle cells, protein digestion and absorption, and insulin secretion ( Figure 4E ). These pathways show lower activity in low-risk PAAD patients, particularly in pathways related to nerve signaling, cytoskeleton, metabolism, and extracellular matrix. Downregulation of these pathways may limit tumor cell proliferation, migration, and nutrient acquisition, thereby slowing tumor aggressiveness and progression.

After enrichment analysis of differentially expressed genes in high and low risk groups of PAAD, we found that up-regulated genes were significantly enriched in immune response, cell differentiation and digestion. Among them, GO terms related to immune response stand out, suggesting that there may be important molecular and cellular changes in the immune microenvironment in high-risk PAAD patients. Given that the immune system plays a key role in the occurrence, development and prognosis of tumors, it is necessary to further explore the clinical significance and biological characteristics of these immune-related genes. Therefore, our next step is to focus on screening for immune-related genes in these differential genes and performing survival analyses on them to assess their impact on the prognosis of patients with PAAD.

Based on the ImmPort database, we identified 113 immune-related DEGs among the differentially expressed genes in the high-low risk group. Through univariate Cox regression analysis and Kaplan-Meier survival analysis, we further screened 11 immune-associated DEGs that were significantly associated with OS in PAAD patients. Among them, CST4, GREM1 and SLURP1 were favorable factors, while PENK, INSL5, KL, PRLR, SCG2, SLC22A17, TAFA2 and VGF were risk factors ( Figure 5A ).

In addition, to fully understand the role of immunity in PAAD progression, we also analyzed differences in immune function and immune infiltration between high and low risk groups. Immune function analysis showed that in the low-risk group, APC co-inhibition, APC co-stimulation, immune checkpoint, and T cell co-inhibition were highly active ( Figure 5B ). The high activity of these immune functions may indicate that the immune system of patients in the low-risk group achieves a balance between anti-tumor response and autoimmune protection. Enhanced APC and T cell inhibitory signaling, as well as regulation of immune checkpoints, help maintain the homeostasis of the immune microenvironment, thereby inhibiting tumor progression.

Immune infiltration analysis revealed a significant increase in CD56 dim natural killer cells and type 17 T helper cells in the high-risk group, potentially leading to a stronger pro-inflammatory response and an immune escape environment that accelerates malignant progression of tumors. In contrast, in the low-risk group, central memory CD4/CD8 T cells, effector memory CD4/CD8 T cells, eosinophils, gamma delta T cells, immature B cells, macrophages, mast cells, myeloid suppressor cells (MDSC), memory B cells, natural killer cells, natural killer T cells, plasmacytoid dendritic cells, regulatory T cells, and T follicular helpers) cell and type 1 T helper cells were more infiltrated ( Figure 5C ). These enhanced infiltrations of memory and effector immune cells, along with moderate immunomodulatory mechanisms, help suppress tumor progression and maintain anti-tumor immune surveillance.

Overall, the high and low risk groups showed significant differences in immune function and immune cell infiltration. In the high-risk group, an increase in pro-inflammatory immune cells may lead to a more aggressive tumor microenvironment; In the low-risk group, moderate immune balance and diversified immune cell infiltration may contribute to tumor suppression. These results provide important clues for understanding the immune microenvironment of PAAD and its impact on patient prognosis, and may provide a basis for personalized immunotherapy for PAAD patients.

Many studies have shown that patients with high expression levels of CD274 or CTLA4 may benefit more from immunotherapy (12, 13). Based on the above analysis of immune characteristics and tumor microenvironment in the high and low risk group of PAAD patients, we further investigated the differences of CD274 and CTLA4, two important immune checkpoint molecules, between the high and low risk groups. We observed significant differences in the expression of CD274 and CTLA4 in the high-low risk group, and the expression of CD274 and CTLA4 in the low-risk group was higher than that in the high-risk group ( Figures 6A, B ). This suggests that the low-risk group may have a relatively mild immune microenvironment compared to the high-risk group, rather than an overactivated pro-inflammatory environment. The high expression of CD274 and CTLA4 can reduce the overreaction of the immune system, thereby inhibiting the release of pro-inflammatory cytokines, and may help delay the malignant progression of tumors. At the same time, we also made TIDE predictions. The results showed that TIDE scores were higher in the low-risk group than in the high-risk group, with higher TIDE scores generally indicating a stronger immune escape capacity and a poorer response to immunotherapy ( Figure 6C ). However, the high TIDE score in patients in the low-risk group may be mainly caused by high expression of CD274 and CTLA4, and the expression of this immune checkpoint is targetable. Therefore, a high TIDE score in the low-risk group is not necessarily a marker of a malignant prognosis, but may instead mean that these patients are more sensitive to CD274 or CTLA4 inhibitors.

In addition, we assessed differences in sensitivity to multiple antineoplastic drugs in high-low risk groups ( Figure 6D ). The results showed that patients in the high-risk group had a high sensitivity to Trametinib, Dabrafenib, SCH772984, ML323, indicating that patients in the high-risk group were more sensitive to these chemotherapy agents, suggesting that PAAD patients in the high-risk group may benefit more from these drugs. In contrast, high-risk patients were insensitive to drugs such as Staurosporine, NU7441, O-3306, Rapamycin, BI-2536, GSK269962A, Fuverastine, AZ960, AZD2014, AZD1332, Rusolitinib, Uprosertib, Alpelisib, Taselisib, WA4003, I-BET-762, RVX-208, OTX015, Entospletinib, AZD5153, CDK9-5576, CDK9-5038, IGF1R-3801, JAK-8517, Carmustine, AZD5363, AZD8186, Cediranib, I-BRD9, telomerase Inhibitor IX, Uni-77, Foretinib, Pyridostatin, AMG-319, BMS-754807, and JQ1. The difference in sensitivity between different drugs further highlights the significant differences in tumor microenvironment and biology between the high and low risk groups, and also suggests potential directions in individualized treatment options.

Overall, there were significant differences in immune checkpoint gene expression, immune escape ability, drug sensitivity, and immune function and infiltrating cells in the high and low risk groups for PAAD. These differences not only deepen our understanding of the immune microenvironment of PAAD, but also provide a valuable basis for personalized immunotherapy. Future studies should further explore the practical application value of these immune features in patients with PAAD, with a view to optimizing the treatment of patients and improving the treatment effect and survival rate of PAAD.