3.1 NK cell-related immune microenvironmental crosstalk in PAAD

We identified NK cell-specific genes following the workflow outlined in Fig. 1A. Specifically, we analyzed single-cell data from the PAAD_GSE162708 dataset, performing cell clustering and differential analysis to identify genes differentially expressed across various cell populations (Fig. 1B–D). Enrichment analysis of these differential genes using Hallmark Gene-Sets revealed that NK cells are involved in a wide range of PAAD-related functions, including cellular components, development, DNA damage, immune responses, metabolic pathways, proliferation, and signaling (Fig. 1E). KEGG enrichment analysis highlighted both known functions and pathways associated with NK cells, such as Natural killer cell mediated cytotoxicity, complement and coagulation cascades, and ECM-receptor interaction. Additionally, it uncovered potentially novel roles in PAAD, such as involvement in ribosome, spliceosome, and T cell receptor signaling pathway, suggesting that NK cells may have a unique role in regulating RNA and mediating T cell signaling in PAAD (Fig. 1F). Transcription factors play a crucial role in PAAD progression [20], so we sought to identify NK cell-specific transcription factors to further elucidate their potential functions. Our results identified key transcription factors with specific expression in NK cells, including KMT2 A, MED1, BRD4, ERG, E2 F6, MAF, HEY1, H2 AFZ, STAT1, and CDK9 (Fig. 1G, H). Since intercellular communication is likely crucial in PAAD, we employed CellChat analysis to explore potential interactions between NK cells and other cell clusters. We found that Myofibroblasts, Fibroblasts, and CD8 T cells are the most closely connected clusters to NK cells (Fig. 1I). Additionally, as donors, NK cells may regulate Myofibroblasts, Fibroblasts, and CD8 T cells through the ADGRE5-CD55 interaction (Fig. 1J). Conversely, as recipients, NK cells might be strongly regulated by B cells, CD8 T cells, Endothelial cells, malignant cells, and Myofibroblasts via HLA-E-KLRC2 and HLA-E-CD94/NKG2 C interactions (Fig. 1K). These findings suggest that NK cells not only act as effectors in the traditional B cell and T cell immune system but may also regulate the functions of these cells.

Fig. 1

NK cell-related immune microenvironmental crosstalk in PAAD. A A detailed flowchart outlines the process for identifying and validating the prognostic model associated with NKDEGs. B The cellular map depicts the distribution and abundance of different cell subgroups (B, CD8 T, Endothelial, Fibroblasts, Malignat, Mast, Mono/Macro, Myofibroblasts, NK cells) within PAAD. C The pie chart shows the precise counts of each cell subgroup (B, CD8 T, Endothelial, Fibroblasts, Malignat, Mast, Mono/Macro, Myofibroblasts, NK cells) in PAAD. D The bar graph illustrates the proportional distribution of each cell subgroup (B, CD8 T, Endothelial, Fibroblasts, Malignat, Mast, Mono/Macro, Myofibroblasts, NK cells) among individual PAAD patients. E The heatmap presents the HALLMARK gene sets regulated by each cell subgroup (B, CD8 T, Endothelial, Fibroblasts, Malignat, Mast, Mono/Macro, Myofibroblasts, NK cells). F The heatmap displays the KEGG gene sets regulated by each cell subgroup (B, CD8 T, Endothelial, Fibroblasts, Malignat, Mast, Mono/Macro, Myofibroblasts, NK cells). G The heatmap highlights key transcription factors that are differentially expressed across various cells within PAAD. H The dot plot reveals significantly expressed transcription factors in NK cells of PAAD, including KMT2 A, MED1, BRD4, ERGE2 F6, MAF, HEY1H2 AFZ, STAT1, and CDK9. I Using CellChat, the interaction probability between specific cell groups and other cell groups is visualized. J Using CellChat, the interaction probability between NK cells as donors and specific gene pairs in other cells is illustrated. K Using CellChat, the interaction probability between NK cells as recipients and specific gene pairs in other cells is depicted

3.2 Construction of NKDEGs prognostic model of PAAD based on multiple machine learning models

To elucidate the role of NKDEGs in PAAD, we first analyzed all differentially expressed genes in PAAD tissues using the TCGA-PAAD dataset, and presented the expression levels of the top 50 genes across all samples (Fig. 2A). Additionally, we performed univariate Cox regression analysis and identified 29 NKDEGs strongly associated with PAAD prognosis (Fig. 2B). Leveraging these NKDEGs, we developed accurate prognostic models for PAAD using 101 different machine learning algorithms, and validated these models with data from the PACA-AU dataset in the ICGC database. The results indicated that NKDEGs demonstrated robust prognostic prediction capability for PAAD patients in both the TCGA-PAAD and ICGC validation datasets, with the CoxBoost + Ridge model achieving the highest predictive accuracy (mean AUC of 0.69) (Fig. 2C, Table S1). Consequently, we selected the CoxBoost + Ridge model to construct our prognostic model, extracting weights for 11 specific NKDEGs to calculate the risk scores for each PAAD sample (Table S2). Survival curves revealed that the risk scores effectively predicted PAAD patient outcomes in both training and validation sets (Fig. 2D, E). Moreover, ROC curves confirmed that the NKDEGs prognostic model has excellent predictive accuracy (Fig. 2F, G).

Fig. 2

Establishment of prognostic model of NKDEGs in PAAD. A The heatmap illustrates the expression levels of NKDEGs associated with PAAD prognosis in both PAAD samples and paired adjacent non-tumor samples. B The univariate Cox regression analysis highlights the NKDEGs related to PAAD prognosis. C The NKDEGs prognostic model was developed using 101 machine learning models (training set: TCGA, validation set: ICGC). D The survival curves for the TCGA cohort are shown based on risk scores derived from the NKDEGs prognostic model. Based on patient risk scores assessed by the model, the graphs indicate that the survival prognosis for the high-risk group (red line) is significantly worse than that for the low-risk (blue line). E The survival curves for the ICGC cohort are depicted based on risk scores derived from the NKDEGs prognostic model. Based on patient risk scores assessed by the model, the graphs indicate that the survival prognosis for the high-risk group (red line) is significantly worse than that for the low-risk (blue line). F The ROC curve for the NKDEGs prognostic model risk scores in the TCGA cohort is presented. G The ROC curve for the NKDEGs prognostic model risk scores in the ICGC cohort is shown

3.3 Clinical correlation analysis of NKDEGs prognostic model for PAAD

Next, we analyzed the clinical relevance and potential clinical utility value of the NKDEGs prognostic model. Univariate Cox regression analysis revealed that age, grade, and risk score are all prognostic factors for PAAD. Notably, the hazard ratio (HR) for the risk score (162.964; 95% CI: 41.361–642.075) was significantly higher than that of other clinical pathological factors (Fig. 3A). Multivariate Cox regression analysis indicated that only age and risk score served as independent prognostic factors for PAAD (Fig. 3B). Concordance analysis demonstrated that the risk score accurately predicted the survival of PAAD patients across different time frames, outperforming other clinicopathological factors (Fig. 3C). Subsequently, we incorporated the significant variables, age and risk score, from the multivariate Cox regression analysis into a nomogram. This nomogram provides a straightforward and precise method for predicting the survival of PAAD patients, thereby offering valuable guidance for clinical decision-making (Fig. 3D, E). Clinical relevance analysis showed that PAAD patients with high-risk scores were more likely to have advanced tumor stage and higher T classification (Fig. 3F, G). There was also an apparent association between high-risk scores and high grade, although it did not reach statistical significance due to limited sample size (p = 0.054) (Fig. 3F, G). Furthermore, we assessed the applicability of the NKDEGs prognostic model across various clinical subgroups. The results indicated that the NKDEGs prognostic model exhibited strong predictive capability for PAAD patient survival across different age, gender, grade, and stage groups (Fig. 3H–K), highlighting the model’s consistency within diverse clinical subpopulations.

Fig. 3

Clinical correlation analysis of NKDEGs prognostic model in PAAD. A Single-factor Cox regression analysis examines the relationship between individual variables(age, gender, tumor grade, cancer stage, and risk scores) and overall survival (OS) in PAAD patients. B Multi-factor Cox regression analysis assesses whether risk scores and other clinicopathological factors can serve as independent prognostic factors for PAAD. C The C-index curve compares the concordance index (C-index) of risk scores with other clinical variables (age, gender, tumor grade, and stage). Risk scores demonstrate the highest predictive power for PAAD prognosis over time. D The nomogram integrates age and risk scores to predict individual probabilities of 1-, 3-, and 5-year survival. E Displays the agreement between predicted and observed survival rates for 1, 3, and 5 years, validating the model’s predictive accuracy with a high C-index (0.708). F The heatmap visualizes the distribution of clinical (age, gender) and pathological (tumor grade, stage, T, M, N categories) features across low- and high-risk score groups. Notable clustering patterns emerge, showing differences in clinical characteristics based on risk levels. G The bubble plot highlights the proportions and correlations between risk score groups (high vs. low) and clinical variables such as age, gender, grade, and tumor-node-metastasis (TNM) classification. H Kaplan–Meier survival curve stratifies PAAD patients younger than 65 and older than 65, showing significantly worse OS in the high-risk group for both age categories (p < 0.001 and p = 0.006, respectively). I Kaplan–Meier survival curve separates patients by gender (male and female), demonstrating that high-risk patients consistently have poorer OS across both groups (p < 0.001 and p = 0.002). J Kaplan–Meier survival curve analyzes tumor grade (G1 - 2 vs. G3 - 4), with high-risk scores predicting worse OS in both categories (p < 0.001 and p = 0.036). K Kaplan–Meier survival curve examines tumor stage (Stage I–II vs. Stage III–IV), indicating that high-risk patients have worse OS, with significant differences in both early- and late-stage disease (p < 0.001 and p = 0.023). *p < 0.05, **p < 0.01, ***p < 0.001

3.4 Characterization and analysis of NKDEGs prognostic model genes

We further analyzed the genes included in the NKDEGs prognostic model to gain a better understanding of the role of NK cells in PAAD and why NKDEGs are effective in predicting PAAD patient outcomes. Our in-depth analysis of the GSE162708 single-cell dataset revealed that IL32, CLEC2B, UCP2, EVL, and SATB1 are significantly upregulated in NK cells, whereas PSAP, SEZ6L2, TSPAN8, EFNB2, ITGA6, and ARHGAP29 are significantly downregulated in NK cells (Fig. 4A). To validate our findings, we selected IL32, the most representative NK-related differentially expressed gene, for experimental verification due to its highest fold change in NK cells (log2 FC = 1.56) and predominant expression in this cell type. qRT-PCR confirmed that IL32 expression was significantly elevated in the NK cell line NK- 92 compared to PAAD cell lines, including BXPC- 3, PANC- 1, SW1990, ASPC1, and COLO357 (Figure S1 A). Functional assays using IL32 knockdown and overexpression in NK- 92 cells co-cultured with PANC- 1 demonstrated that IL32 overexpression promoted PAAD cell proliferation, whereas knockdown suppressed it (Figure S1B). These findings provide compelling evidence that IL32, highly expressed in NK cells, plays a crucial role in PAAD progression. Survival analysis indicated that high expression of IL32, CLEC2B, TSPAN8, EFNB2, ITGA6, and ARHGAP29 correlates with poorer survival in PAAD patients, while the remaining genes are associated with better survival outcomes (Fig. 4B). We also investigated the correlations among the NKDEGs prognostic model genes to assess their potential synergistic effects. Our findings showed that SATB1 and ARHGAP29 have the highest correlation (0.53), while TSPAN8 and EFNB2 exhibit correlations with almost all other NKDEGs prognostic model genes, suggesting they may function as hub genes (Fig. 4C). Additionally, using GeneMANIA, we characterized the specific interactions and functional relationships among the NKDEGs prognostic model genes. The analysis revealed that these genes primarily regulate pathways related to the regulation of cell activation and positive regulation of leukocyte activation, which further underscores the potential role of NK cells in PAAD.

Fig. 4

Characteristics of NKDEGs prognostic model genes. A A comprehensive reanalysis of the single-cell dataset GSE162708 showcases the expression patterns of NKDEGs prognostic model genes in various cell clusters. B Survival analysis curves visually represent the significant influence of NKDEGs prognostic model genes on the outcomes for PAAD patients. C The gene correlation analysis reveals the interrelationships among the NKDEGs prognostic model genes. D Insights derived from the Genemania database analysis illustrate the functional regulation of NKDEGs prognostic model genes within the context of PAAD. *p < 0.05, **p < 0.01, ***p < 0.001

3.5 Functional enrichment in different risk score groups for PAAD

To investigate the unique characteristics of various risk score groups and clarify the potential role of NK cells in PAAD regulation, we conducted a differential gene expression analysis using the TCGA-PAAD dataset. This analysis identified 58 upregulated genes and 210 downregulated genes (Fig. 5A), with the top fifty differentially expressed genes highlighted (Fig. 5B). Subsequent Gene Ontology (GO) analysis suggested that NK cells may play a role in key cellular molecular functions and biological processes in PAAD, including membrane potential regulation, signal release, and hormone secretion (Fig. 5C, D). We then performed KEGG enrichment analysis to identify pathways differentially represented across risk score groups. This analysis revealed significant pathways such as Neuroactive Ligand-Receptor Interaction and the MAPK signaling pathway (Fig. 5E). Finally, Gene Set Enrichment Analysis (GSEA) was used to assess the importance of these enriched pathways. The results showed that, in the high-risk score group, pathways like aminoacyl-tRNA biosynthesis and the cell cycle were prominently enriched (Fig. 5F). Conversely, in the low-risk score group, the enriched pathways included primary immunodeficiency and steroid hormone biosynthesis (Fig. 5F). This comprehensive analysis provides valuable insights into the molecular mechanisms underlying different risk profiles in PAAD and the involvement of NK cells.

Fig. 5

Functional enrichment analysis of different risk score groups. A The volcano plot displays the genes with differential expression across high and low risk patients. Genes highlighted in red have a fold change greater than 2 and an FDR less than 0.05, while those in green have a fold change less than − 2 and an FDR below 0.05. B The heatmap illustrates the expression levels of the top fifty differentially expressed genes among different risk score groups in PAAD samples. C The bubble plot highlights the pathways from the GO enrichment analysis, showing their proportions across different risk score groups. Key pathways include those involved in immune response regulation, extracellular matrix organization, and cellular metabolic processes. D The circular diagram provides a comprehensive view of GO enrichment analysis, categorizing pathways into three domains: Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF). Each section indicates the number of associated genes, with statistical significance (adjusted p values) represented by a color gradient. This visualization emphasizes the interconnected roles of pathways, including those involved in immune system processes, signal transduction, and tumor microenvironment remodeling. E The bubble plot summarizes the results of KEGG pathway enrichment analysis for different risk score groups. Notable pathways include the MAPK signaling pathway, neuroactive ligand-receptor interaction, and insulin secretion, which are closely linked to tumor growth, metastasis, and metabolic dysregulation in PAAD. F GSEA analysis identifies the key pathways enriched, including epithelial-mesenchymal transition (EMT) and immune-related processes in PAAD patients. G The box plot compares the tumor mutation burden (TMB) between high- and low-risk groups. High-risk patients exhibit a significantly elevated TMB (p = 0.033), suggesting that genetic instability correlates with higher risk scores. H Pearson correlation analysis reveals the relationship between risk scores and tumor mutation burden in PAAD patients. The positive correlation (R = 0.18, p = 0.007) indicates that higher risk scores are associated with increased TMB, suggesting a potential link between genetic mutations and risk stratification. I Analysis of somatic mutation data highlights the mutation rates of key genes in PAAD patients. The color-coded mutations (e.g., missense mutations, frameshift deletions) highlight the mutation frequencies, with high-risk patients exhibiting higher overall mutation rates. J The Kaplan–Meier survival curve evaluates the prognostic impact of TMB alone. Patients are divided into high- and low-TMB groups, with high TMB associated with worse overall survival (p = 0.009). K The Kaplan–Meier curve combines TMB and risk score groups to assess their joint impact on survival. Four subgroups are defined: TMB-high/risk-high, TMB-low/risk-high, TMB-high/risk-low, and TMB-low/risk-low. Patients in the TMB-high/risk-high group exhibit the poorest survival, while the TMB-low/risk-low group has the best prognosis (p < 0.001). *p < 0.05, **p < 0.01, ***p < 0.001

Mutations in critical genes within tumor cells are strongly linked to the progression and prognosis of PAAD [21]. We sought to investigate whether the functionality of NK cells in PAAD is associated with these mutations. Our findings revealed that patients in the high-risk score group exhibited a higher tumor mutation burden (TMB), with a positive correlation between risk scores and TMB (Fig. 5G, H). Further analysis showed that individuals with high-risk scores had notably higher somatic mutation rates in essential genes, such as KRAS (72% vs. 47%) and TP53 (65% vs. 49%), compared to those with low-risk scores (Fig. 5I). Additionally, PAAD patients with elevated TMB experienced significantly poorer prognoses than those with lower TMB (Fig. 5J). When evaluating prognosis using both TMB and risk scores, patients with both high TMB and high-risk scores had the poorest outcomes, while those with low TMB and low-risk scores had the best outcomes. Patients with either high TMB or high-risk scores alone had intermediate prognoses (Fig. 5K). These results suggest a potential link between NK cell functionality and somatic mutations in PAAD, and indicate that combining the NKDEGs prognosis model with TMB could improve the prediction of patient outcomes.

3.6 Assessment of immune microenvironment in different risk score groups of PAAD

NK cells, as essential immune components, are likely crucial in regulating the immune microenvironment of pancreatic ductal adenocarcinoma (PAAD). An analysis of the immune microenvironment reveals that patients with high-risk scores in PAAD tend to have a lower ESTIMATE score (p = 0.051), though there are no significant differences in StromalScore and ImmuneScore (Fig. 6A). Regarding immune cell infiltration, a significant negative correlation exists between risk scores and NK cell infiltration, with various B cells and T cells displaying similar trends to NK cells (Fig. 6B, C). Using the CIBERSORT algorithm, we identified a significant increase in activated NK cells accompanied by a concurrent depletion of resting NK cell subsets in high-risk PAAD patients. Despite the elevated NK cell levels observed in this group, their prolonged activation—likely induced by sustained exposure to immunosuppressive cytokines (e.g., TGF-β, IL- 10) and metabolic stress within the tumor microenvironment (TME)—may drive aberrant inflammatory responses and contribute to an unfavorable prognosis. Additionally, the responses of T cells and B cells exhibited an inverse trend relative to the overactivation of NK cells (Fig. 6D). Further immune correlation analysis revealed a marked negative correlation between activated NK cells and the responses of B cells and T cells, indicating potential intercellular communication, consistent with findings in Fig. 1 (Fig. 6E). Additionally, we explored the potential regulatory role of NKDEGs prognostic model genes in intercellular communication and identified SEZ6L2, CLEC2B, and EFNB2 as key genes potentially regulating NK cell interactions with other cells (Fig. 6F).

Fig. 6

Assessment of the immune microenvironment in different risk score groups. A The violin plot illustrates the StromalScore, ImmuneScore, and ESTIMATE scores across different risk score groups. B The heatmap shows the distribution and proportions of immune cells across various risk score groups, as determined by the CIBERSORT algorithm. C The bubble plot reveals the correlation between risk scores and immune cell infiltration in PAAD samples. D The box plot displays the abundance of immune cells across different risk score groups, as calculated using the CIBERSORT algorithm. E The heatmap highlights the correlations between various immune cells and functions in PAAD. F The heatmap illustrates the regulatory interactions of NKDEGs prognostic model genes with different immune cells and functions in PAAD. G The box plot presents the expression levels of immune checkpoint genes across different risk score groups. H Survival analysis evaluates the response and efficacy of immune therapy across different risk score groups within the IMvigor210 immune therapy cohort. I The box plot shows the variation in risk scores among different immune therapy response cohorts. J The association between TCGA immune subtypes and NKDEGs prognostic subtypes is examined. *p < 0.05, **p < 0.01, ***p < 0.001

We next evaluated the potential of the NKDEGs prognostic model in guiding clinical immunotherapy. Our analysis revealed that several immune checkpoint genes are upregulated in high-risk patients, suggesting that these individuals with PAAD might be more responsive to immune checkpoint inhibitors (Fig. 6G). Additionally, drug sensitivity analysis indicated that high-risk patients might respond better to targeted therapy with BI- 2536, while those in the low-risk group are more likely to benefit from AZD5363, Afuresertib, and BMS- 754807 (Figure S2 A). In the IMvigor210 immunotherapy cohort, we confirmed that high-risk patients experience notably poorer outcomes from immunotherapy (Fig. 6H). Furthermore, we observed a trend indicating that patients with immune therapy responses (CR/PR) tend to have lower risk scores compared to non-responders (SD/PD), although this trend did not reach statistical significance due to sample size limitations (Fig. 6I). Moreover, the TCGA official team has classified PAAD into four subtypes: Wound Healing (Immune C1), IFN-gamma Dominant (Immune C2), Inflammatory (Immune C3), and TGF-beta Dominant (Immune C6), which has become a significant reference [22]. Our findings demonstrate that the NKDEGs prognostic model can effectively differentiate among these immune subtypes, particularly showing a significant increase in Immune C1 patients within the high-risk group, and a notable decrease in Immune C3 and C6 patients (Fig. 6J).

3.7 Novel molecular classification based on 11-NKDEGs

The therapeutic plateau of current PAAD immunotherapies, evidenced by a universal 100% non-durable response rate, necessitates novel stratification paradigms accounting for immune heterogeneity [23]. Leveraging our 11-NKSG prognostic signature, we implemented a consensus clustering framework (ConsensusClusterPlus) to delineate molecular subtypes. Cluster stability optimization via Monte Carlo cross-validation identified k = 4 as the inflection point, demonstrating minimal cluster disassociation index and maximum proportional area change in cumulative distribution function analysis (Fig. 7A, B). The sample distribution across different k values is presented in Fig. 7C, with the consistency matrix heatmap showing consistent blue shading when k = 4 (Fig. 7D). Kaplan–Meier survival analysis delineated a striking dichotomy: the C4 subtype exhibited superior median overall survival, while C1-C3 shared comparable mortality trajectories (Fig. 7E). A Sankey diagram illustrated the connection between the novel molecular subtypes and the prognosis model categories (Fig. 7F). Both PCA and tSNE analyses clearly differentiated PAAD patients according to the newly identified molecular subtypes (Fig. 7G, H). Further, we examined whether the molecular subtyping approach could distinguish PAAD patients based on their immune characteristics. Immune microenvironment analysis revealed that the C1 group exhibited elevated ESTIMATE, Stromal, and ImmuneScores, while patients in the C4 group had lower scores, and those in the C2 and C3 groups showed intermediate values (Fig. 7I). Tumor Purity exhibited an inverse trend (Fig. 7I). Regarding immune cell infiltration, the C4 group exhibited the most pronounced NK cell infiltration and T cell activation, suggesting that the potent anti-tumor effects of NK cells might be a key factor in their favorable prognosis (Fig. 7J). Conversely, C1 group patients showed extensive activation of various immune cells, which might be beneficial for tumor killing but could also contribute to the poor prognosis due to excessive inflammation (Fig. 7J). In contrast, C3 group patients exhibited suppressed immune cells, representing an “immunologically desert” state with a lack of anti-tumor immune response (Fig. 7J). The C2 group showed no distinct immune cell characteristics, which might indicate other molecular mechanisms contributing to their poor prognosis (Fig. 7J). Furthermore, we analyzed the expression of immune checkpoint genes across different molecular subtypes. Patients in each PAAD molecular subtype could potentially benefit from targeted immune checkpoint inhibitors based on their characteristic gene expression profiles (Fig. 7K). We also identified potential targeted drugs suitable for each molecular subtype. C1 group patients were most likely to benefit from GSK343, C2 group patients from Sapitinib, C3 group patients from AZD6738, and C4 group patients from Linsitinib (Figure S3).

Fig. 7

Novel molecular subtyping for identification of PAAD based on NKDEGs. A Delta area curves for varying numbers of classifications. B Cumulative Distribution Function (CDF) curves for different classification numbers, where each curve represents the model stability at different K values. C The heatmap displays the distribution of PAAD patients across different K values. D The consistency clustering plot shows the clustering results when K equals 4. E Survival curves illustrate the prognosis of PAAD patients with different molecular subtypes. F The Sankey diagram illustrates the correspondence between different molecular subtypes of PAAD patients and various risk score groups. G PCA analysis reveals the distribution of samples across different PAAD molecular subtypes. H t-SNE analysis shows the distribution of samples across different PAAD molecular subtypes. I The box plot demonstrates the ESTIMATE score, StromalScore, ImmuneScore, and tumor purity across different PAAD molecular subtypes. J The heatmap displays immune cell infiltration status across different PAAD molecular subtypes, based on various algorithms. K The box plot shows the expression levels of immune checkpoint genes across different PAAD molecular subtypes. *p < 0.05, **p < 0.01, ***p < 0.001