This article is mentioned in:

Introduction

Pancreatic adenocarcinoma (PAAD) is a highly aggressive malignancy characterized by a dismal prognosis and significant economic burden on healthcare systems worldwide (1,2). The disease is often diagnosed at advanced stages, leading to limited treatment options and a five-year survival rate of <10%. Current therapeutic strategies, including surgical resection, chemotherapy and targeted therapies, have shown limited efficacy, primarily due to the tumor's complex microenvironment and the presence of resistance mechanisms (2). Despite advancements in understanding the molecular underpinnings of PAAD, there remains a critical gap in identifying reliable prognostic biomarkers that can stratify patients based on their disease status and guide personalized treatment approaches (3). The complex tumor microenvironment and significant heterogeneity among PAAD subtypes contribute to this bleak outlook, making it imperative to identify molecular signatures that can aid in both diagnosis and treatment decision-making.

Recent advancements in machine learning have opened new frontiers in cancer biology, allowing for the identification and characterization of biomarkers that correlate with clinical outcomes (4,5). One intriguing area of research is the exploration of damage-associated molecular patterns (DAMPs), which are endogenous signals released by stressed or damaged cells, playing critical roles in inflammation and tumor progression (6,7). DAMPs, including high-mobility group (HMG) box 1, heat shock proteins and other nucleic acids, can activate innate immune responses and modulate the tumor microenvironment (8,9). Their contributions to cancer pathogenesis and prognosis have sparked interest in establishing a predictive score, termed DAMPscore, to classify cancer patients based on their distinct molecular profiles (10).

The present study focuses on the expression of DAMP-related genes in PAAD and their potential prognostic implications. Previous research has established a link between DAMPs and tumor progression (11), yet the specific role of DAMP-related gene expression in stratifying patients with PAAD remains underexplored. By employing Consensus Clustering, three distinct patient clusters were identified, and it was hypothesized that these clusters reflect varying DAMP statuses, which may correlate with disease prognosis. Notably, the present analysis revealed 141 differentially expressed genes with significant prognostic value across multiple cohorts, underscoring the relevance of DAMP-related gene expression in PAAD. The identification of six pivotal genes through Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression further highlighted the potential of the DAMPscore as a robust prognostic tool, demonstrating superior performance compared to existing models. Additionally, the functional exploration of PLEK2 suggested its involvement in key signaling pathways associated with immune response and tumor dynamics, thereby reinforcing the clinical significance of DAMP-related gene expression in PAAD.

Collectively, these findings illuminate the intricate relationship between DAMPs and PAAD, presenting a novel avenue for enhancing prognostic assessments and therapeutic strategies in this challenging malignancy.

Materials and methods

Public data acquisition and processing

The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) was used to retrieve the TCGA-PAAD cohort. The datasets GSE85916, GSE71729 and GSE62452 were sourced from the Gene Expression Omnibus (GEO) database (https://ncbi.nlm.nih.gov/gds). The ICGC_AU cohort was collected from the International Cancer Genomics Consortium (ICGC) data portal (https://dcc.icgc.org/). The E-MTAB-6134 cohort was downloaded from the ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). The mRNA data and proteomic data of PLEK2 were downloaded from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) [https://hupo.org/Clinical-Proteome-Tumor-Analysis-Consortium-(CPTAC)]. Immunohistochemical images of PLEK2 was directly downloaded from The Human Protein Atlas (HPA; http://www.proteinatlas.org/). All datasets were processed according to the methods outlined in a previous study by our group (12). In addition, the three GEO datasets (GSE85916, GSE71729 and GSE62452) were integrated into a Meta-cohort using the sva package in R (version 4.1.3) for further analyses. All datasets utilized in this study were obtained from public repositories, thereby negating the need for additional ethical approval.

Consensus clustering

Consensus clustering was performed utilizing the ‘ConsensusClusterPlus’ package in R (13). The analysis was executed with the following parameters: The maximum number of clusters (k) (MaxK)=10; the number of resampling iterations (reps)=1,000, the proportion of samples randomly selected in each resampling iteration (pItem)=0.8, the proportion of features (variables) retained in each resampling iteration (pFeature)=1, the clustering algorithm applied (clusterAlg)=‘pam’, the method for handling missing data during correlation calculations (corUse)=‘complete.obs’ and seed=123,456. The optimal number of clusters was determined based on the following criteria: i) Cumulative distribution function (CDF) Plot: The CDF plot, generated by the ConsensusClusterPlus package, was analyzed to assess the proportion of variance explained by the clustering solution for different values of k. The plot was examined for an ‘elbow point’ or stabilization in the curve, which indicates the point beyond which increasing the number of clusters does not significantly improve the clustering solution. ii) Consensus matrix: The consensus matrix, which measures the frequency with which pairs of samples are clustered together across multiple resampling iterations, was also carefully examined. A clear and stable block-like structure in the consensus matrix was sought to identify well-defined clusters.

Construction of DAMPscore in PAAD

Differentially expressed genes (DEGs) among the subgroups identified through consensus clustering were extracted using the ‘limma’ package in R. Following this, these DEGs were subjected to univariate Cox regression analysis. Genes that demonstrated a prognostic P<0.05 in both the TCGA-PAAD and Meta-cohort were incorporated into a LASSO regression model applied to the Meta-cohort. The genes selected were then analyzed using a multivariate Cox regression model. The key genes and their respective coefficients were utilized to calculate the risk score for each patient using the following formula: Score=0.12012 × ADM + 0.05077 × HMGA2 + 0.10670 × EREG + 0.16538 × SLC16A3 + 0.06207 × LAMC2 + 0.27480 × IGF2BP2. The DAMPscore for patients in each cohort was computed using the following formula: DAMPscore=(score-Min)/abs(Max), in accordance with previous research findings (12,14).

Weighted gene co-expression network analysis (WGCNA)

WGCNA was performed using R software, as detailed in previous studies by our group (15,16). In summary, transcriptional expression data from each dataset were utilized to construct a gene co-expression network, following the exclusion of genes and samples with excessive missing values. The correlation between nodes was assessed by creating an adjacency matrix using the following formula:

where Sij represents the co-expression similarity, indicating the Pearson correlation coefficient between two distinct genes, i and j. xi and xj denote the transcriptional expression values for genes i and j, while aij reflects the correlation strength between these two genes. A scale-free R2 threshold of 0.9 was established for selecting the appropriate soft-threshold β. Subsequently, one-step network construction and module detection methods were applied, with a minimum module size of 200 and a mergeCutHeight set at 0.25 for module merging. Finally, module-trait associations were evaluated to pinpoint modules that were significantly related to the DAMPscore. In addition, the definitions and calculations for module eigengenes, gene significance and module significance followed the methodology described in a prior study (17).

Enrichment analysis

Gene Set Enrichment Analysis (GSEA) was performed for patients with pancreatic cancer using the ‘clusterProfiler’ package in R. The results of the GSEA were visualized with the ‘GseaVis’ package in R (18). The analysis utilized the hallmark gene sets database (https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp#H).

Cell culture

AsPC-1 cells were sourced from the Cell Bank of the Chinese Academy of Sciences, whereas PANC-1 cells were obtained from the American Type Culture Collection. Prior to experimentation, both cell lines were verified to be free of mycoplasma contamination. The cells were cultured in RPMI-1640 medium (Invitrogen; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; ExCell), and were incubated in a humidified atmosphere with 5% CO2 at 37°C.

Design, synthesis and transfection of small interfering RNA (siRNA)

The coding DNA sequences of the target genes were obtained from the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/) and siRNAs were designed for each gene using the DSIR design website (http://biodev.extra.cea.fr/DSIR/DSIR.php). The synthesis was commissioned from Shanghai Jima Pharmaceutical Technology Co., Ltd. (detailed information can be found in Table SI). Cells were seeded in 6-well plates and transfected when they reached ~80% confluency. SiRNAs and plasmids were transfected into the cells using the Lipo8000™ transfection reagent (cat. no. C0533; Biyuntian) according to the manufacturer's instructions.

Cell Counting Kit-8 (CCK-8) assay

Cells were seeded in a 96-well plate, transfected with lentivirus for 72 h and then subjected to continuous selection with puromycin (5 µg/ml) for one week. After that, the cells were expanded and subjected to the CCK-8 assay. A working solution was prepared by mixing cell culture medium with CCK-8 reagent (cat. no. C0038; Beyotime Institute of Biotechnology) at a 10:1 ratio at 37°C in the dark for 1 h. After discarding the residual medium, 100 µl of the working solution was added to each well and the cells were incubated at 37°C for 1 h prior to measurement with a multifunctional microplate reader (SynergyHTX; Gene Company Ltd.).

Transwell and invasion assay

For the migration (Transwell) assay, 1×104 cells were seeded in 100 µl of RPMI-1640 medium in the upper chamber of Transwell® Permeable Supports with a polycarbonate membrane (8 mm pore size; Corning, Inc.). In the lower chamber, 800 µl of RPMI-1640 medium supplemented with 30% FBS was added. For the invasion assay, 1×105 cells were seeded in 200 µl RPMI-1640 medium on a Matrigel® layer in the upper chamber of Corning® BioCoat™ Matrigel® Invasion Chambers (8 mm pore size; Corning, Inc.), with 800 µl RPMI-1640 medium containing 20% FBS added to the lower chamber. In both assays, cells were incubated at 37°C with 5% CO2 for 24 h, followed by fixation in 4% paraformaldehyde at 25°C for 20 min and staining with Giemsa Stain solution (product no. 32884; Sigma-Aldrich; Merck KGaA) at 25°C for 20 min. The cells that had transgressed to the bottom of the membrane were visualized under a light microscope and quantified by counting the number of cells in seven randomly selected fields at 200-fold magnification.

RNA extraction and reverse transcription-quantitative (RT-q)PCR

Total RNA was extracted using TRIzol and subsequently quantified with a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific, Inc.). cDNA synthesis was performed according to the manufacturer's instructions (cat. no. A5003; GoScript™ Reverse Transcriptase from Promega Corp.). RT-qPCR was conducted as previously described in earlier studies by our group (19,20). The primer sequences used in the experiments are listed in Table SI. Results were normalized to the housekeeping gene beta-actin, which remains stable across all samples.

Statistical analysis

Data analysis and visualization were performed using various R packages (version 4.1.3), including ‘tidyverse’, ‘ggplot2’, ‘ggrepel’, ‘survival’, ‘survminer’, ‘survivalROC’, ‘timeROC’, ‘rms’, ‘msigdbr’, ‘scales’, ‘dplyr’, ‘ggalluvial’, ‘VennDiagram’, ‘patchwork’, ‘org.Hs.eg.db’, ‘CompareC’, ‘enrichplot’, ‘glmnet’, ‘aplot’, ‘corrplot’ and ‘pec’. The Pearson correlation method was utilized for correlation analyses, while prognosis assessment was conducted using the Kaplan-Meier method with the log-rank test. ROC curve analysis was performed by using the timeROC package in R (version 4.1.3). Comparison between two groups or among three groups was conducted by R software using the Wilcox rank-sum test or ANOVA. P<0.05 was considered to indicate statistical significance.

Results

Unsupervised clustering of DAMP-related genes in PAAD

Drawing upon the expression of 40 DAMP-related genes (Table SII) (21,22), Consensus Clustering was employed to stratify patients with PAAD within the TCGA-PAAD and Meta-cohort datasets. The optimal parameter for further analysis was determined to be k=3, based on the consensus matrix evaluated from k=2 to 10 (Fig. 1A, B, D and E; Fig. S1A and B). In the TCGA-PAAD cohort, 176 patients with PAAD were categorized into three clusters: C1, C2 and C3, with the C3 cluster exhibiting the shortest median overall survival (OS) (Fig. 1C). A similar trend was observed in the Meta-cohort, where significant OS differences were noted among the three clusters, with C3 displaying the worst prognosis (Fig. 1F).

Figure 1. Consensus clustering of patients with PAAD based on DAMP-related genes. (A) The CDF curve of consensus clustering for k=2 to 10 in the TCGA-PAAD cohort. (B) Heatmap depicting consensus clustering solution (k=3) for DAMP-related genes in patients with PAAD in the TCGA-PAAD cohort. (C) Kaplan-Meier curves of OS in the C1, C2 and C3 clusters in the TCGA-PAAD cohort. (D) CDF curve of consensus clustering for k=2 to 10 in the Meta-cohort. (E) Heatmap depicting consensus clustering solution (k=3) for DAMP-related genes in patients with PAAD in the Meta-cohort. (F) Kaplan-Meier curves of OS in the C1, C2 and C3 clusters in the Meta-cohort. (G) Transcriptional expression of DAMP-related genes in the C1, C2 and C3 clusters in the TCGA-PAAD cohort. (H) Transcriptional expression of DAMP-related genes in the C1, C2 and C3 clusters in the Meta-cohort. *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001. ns, no significance; TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma; CDF, cumulative distribution function; DAMP, Damage-Associated Molecular Pattern; OS, overall survival.

Analysis of the transcriptional expression of the majority of immunogenic cell death (ICD)-related genes in the TCGA-PAAD cohort revealed that the C3 cluster was characterized by notably high expression of NT5E, which is associated with impaired ICD (21), alongside low expression of several genes that promote ICD (21,22), such as TLR4, P2RX7 and FPR1 (Fig. 1G). By contrast, the C1 cluster exhibited the highest expression levels of most ICD-facilitating genes, including IL33, FPR1 and TLR7. A comparable expression pattern for these genes was also observed among the three clusters in the Meta-cohort (Fig. 1H). Given the observed differences in prognosis and the expression profiles of ICD-related genes across these three clusters, it may be hypothesized that the C1 cluster reflects a pro-DAMP status and the C3 cluster an anti-DAMP status, while the C2 cluster occupies an intermediate position between the two.

Construction of DAMPscore in PAAD

To further elucidate the key factors influencing the release of DAMP-related molecules in PAAD, differential gene analysis was we conducted on the aforementioned C1 and C3 clusters. As illustrated in Fig. 2A, a total of 1,444 upregulated and 1,030 downregulated genes from the TCGA-PAAD cohort fulfilled the criteria of an absolute fold change of ≥1.5 and an adjusted P<0.05. Similarly, in the Meta-cohort, 763 upregulated genes and 663 downregulated genes met the same threshold (Fig. 2B). Among the differential genes identified from both cohorts, 962 common DEGs were found (Table SIII; Fig. 2C). These common DEGs were subjected to univariate Cox analyses within both the TCGA-PAAD dataset and the Meta-cohort, with 141 genes demonstrating significant prognostic significance (P<0.05) in both datasets (Fig. 2D). To identify the pivotal genes that differentiate between the C1 and C3 clusters, LASSO Cox regression analyses were performed (Fig. 2E and F). Following the determination of the optimal value of λ, six genes were selected for further multivariate Cox regression analysis, which allowed for the calculation of their respective coefficients (Fig. 2G). The six identified genes included (SLC16A3, LAMC2, IGF2BP2, HMGA2, EREG and ADM). Utilizing the transcriptional expressions and corresponding coefficients of these nine genes, the DAMPscore was derived for each sample.

Figure 2. Construction of Damage-Associated Molecular Pattern score in PAAD. (A) Volcano plot of DEGs between the C1 and C3 cluster from the TCGA-PAAD cohort. (B) Volcano plot of DEGs between the C1 and C3 cluster from the Meta-cohort. (C) Venn diagram of common DEGs in the TCGA-PAAD and Meta-cohort. (D) Venn diagram of DEGs with significant prognostic P<0.05 in the TCGA-PAAD and Meta-cohort. (E and F) The Least Absolute Shrinkage and Selection Operator Cox regression model was used for hub gene selection. (E) The tuning parameter (λ) was calculated based on the partial likelihood deviance with ten-fold cross validation. (F) An optimal log λ value shown by the vertical black dot-lines in the plots. The six signature genes were identified according to the best fit profile. (G) Coefficients of the six hub genes. TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma; DEG, differentially expressed gene; FC, fold change; adj.P.Val, adjusted P-value.

Prognostic value of DAMPscore in PAAD

As depicted in Fig. 3A, within the training cohort (Meta-cohort), patients with PAAD presenting a high DAMPscore, based on the median DAMPscore threshold established in the cohort, exhibited markedly shorter median OS compared to those with a low DAMPscore (P<0.0001). Similarly, in the three independent GEO datasets (GSE85916, GSE71729 and GSE62452) that were not integrated, using the median DAMPscore as a cutoff, patients with a higher DAMPscore demonstrated progressively worse prognoses (Fig. 3B-D). Furthermore, the prognostic value of the DAMPscore was validated across three external datasets (E-MTAB-6134, ICGC-AU and TCGA-PAAD). As illustrated in Fig. 3E-G, applying the median DAMPscore as a cutoff across all cohorts revealed that patients with a lower DAMPscore in these three cohorts experienced significantly prolonged median OS.

Figure 3. Prognostic evaluation of the DAMPscore in PAAD. (A-G) Kaplan-Meier curves of overall survival in patients with PAAD with high and low DAMPscore in datasets (A) from the Meta-cohort, (B) GSE85916, (C) GSE71729, (D) GSE62452, (E) E-MTAB-6134, (F) ICGC-AU and (G) TCGA-PAAD. (H) The 1-, 3- and 5-year AUCs of DAMPscore in the TCGA-PAAD, E-MTAB-6134, ICGC-AU and Meta-cohort. DAMP, Damage-Associated Molecular Pattern; TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma; AUC, area under the receiver operating characteristic curve.

To further evaluate the prognostic predictive capability of the DAMPscore in patients with PAAD at the baseline, a ROC curve analysis was performed and the corresponding area under the curve (AUC) values for the 1-, 3- and 5-year survival rates were calculated. In the Meta-cohort, the AUC values for the 1-, 3- and 5-year DAMPscore were 0.65, 0.70 and 0.77, respectively. In the TCGA-PAAD cohort, the AUC values for the same time intervals were 0.71, 0.73 and 0.72. For the E-MTAB-6134 cohort, the AUC values were 0.76, 0.66 and 0.62, while in the ICGC-AU cohort, they were 0.71, 0.70 and 0.89 (Fig. 3H). These findings collectively suggest that DAMPscore demonstrates a substantial prognostic performance in PAAD.

Comparison of DAMPscore and other prognostic models

PAAD exhibits subtype-specific prognostic heterogeneity. Bailey et al (23) established four molecular subtypes: i) Squamous (TP53/KDM6A mutations, TP63ΔN activation, endodermal gene hypermethylation); ii) progenitor (FOXA2/3-PDX1-MNX1 signatures); iii) immunogenic (immunosuppressive pathway enrichment); and iv) ADEX (KRAS-driven exocrine/endocrine differentiation via NR5A2/RBPJL and NEUROD1/NKX2-2). Crucially, the squamous subtype independently predicts poor survival (23). While the DAMPscore distribution showed non-significant inter-subtype variation (P=0.274; Fig. S1C), its gradient pattern suggests a complementary value to molecular subtyping for precision oncology applications. To further assess the prognostic utility of the DAMPscore in patients with PAAD, the concordance index (C-index) for the DAMPscore in comparison to other clinical characteristics, including T stage, N stage and age, was calculated. As depicted in Fig. 4A, the DAMPscore achieved the highest C-index among the evaluated metrics. In addition, 25 previously published mRNA-related prognostic signatures were included for comparative analysis (Table SIV). These signatures are linked to various factors such as gemcitabine resistance, cancer-associated fibroblasts, DNA and RNA methylation, pyroptosis, immune cell infiltration, mitophagy, invasion and other biological processes (24–29). Beyond its moderate C-index in the TCGA-PAAD cohort (Fig. 4B), the DAMPscore exhibited superior performance in the E-MTAB-6134 (Fig. 4C), GSE62452 (Fig. 4D), GSE85916 (Fig. 4E) and Meta-cohort (Fig. 4F). Furthermore, within the ICGC-AU and GSE71729 cohorts, the DAMPscore's C-index ranked third among all prognostic models assessed (Fig. 4G and H).

Figure 4. Comparison of DAMPscore and other factors in prognostic prediction. (A) The C-index (95% confidence interval) of DAMPscore and clinical features in the TCGA-PAAD, E-MTAB-6134, ICGC-AU and GSE62452 cohorts. (B-H) The C-index (95% confidence interval) of DAMPscore and 25 previously published mRNA signatures in the (B) TCGA-PAAD, (C) E-MTAB-6134, (D) GSE62452, (E) GSE85916, (F) Meta-cohort, (G) ICGC-AU and (H) GSE71729. *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001. ns, no significance; DAMP, Damage-Associated Molecular Pattern; TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma.

Identification of DAMPscore-related hub genes in PAAD

To identify hub genes potentially involved in DAMPs release and immune response induction, WGCNA was employed to select gene modules that exhibit a strong correlation with both the DAMPscore and CD8 T-cell infiltration across the TCGA-PAAD (Fig. 5A), GSE71729 (Fig. 5B) and E-MTAB-6134 (Fig. 5C) cohorts. The analyses revealed that the blue module in TCGA-PAAD (r=0.72, P=7×10−29), the black module in GSE71729 (r=0.74, P=3×10−22) and the green-yellow module in E-MTAB-6134 (r=0.85, P=1×10−80) were significantly positively correlated with DAMPscore. In addition, these modules displayed a notable negative correlation with CD8 T-cell infiltration (Fig. 5A-C).

Figure 5. Identification of DAMPscore-related hub genes in PAAD. (A-C) Correlation heatmap of module genes and DAMPscore or CD8 T-cell infiltration in the (A) TCGA-PAAD, (B) GSE71729 and (C) E-MTAB-6134 cohorts. The correlation coefficient varies from −1 to 1 as the color turns from blue to red gradually. (D-F) Scatterplot of the correlation coefficient between the selected module and DAMPscore in the (D) TCGA-PAAD, (E) GSE71729 and (F) E-MTAB-6134 cohorts. (G) Venn plot of hub genes of the selected module in the TCGA-PAAD, GSE71729 and E-MTAB-6134 cohorts. (H) Correlation analyses between the expression of selected hub genes and DAMPscore in the TCGA-PAAD, GSE71729 and E-MTAB-6134 cohorts. ME, module eigengene; DAMP, Damage-Associated Molecular Pattern; TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma; cor, correlation coefficient.

Furthermore, the gene module membership correlation (cor.geneModuleMembership) and gene trait significance correlation (cor.geneTraitSignificance) were computed for the blue module in TCGA-PAAD (Fig. 5D), the black module in GSE71729 (Fig. 5E) and the green-yellow module in E-MTAB-6134 (Fig. 5F) in order to screen for hub genes. The thresholds were set to |cor.geneTraitSignificance| >0.5 and |cor.geneModuleMembership| >0.5. This approach identified a total of 394 genes in the blue module of TCGA-PAAD, 114 genes in the black module of GSE71729 and 38 genes in the green-yellow module of E-MTAB-6134 that met the established criteria (Fig. 5G).

A total of eight hub genes were subsequently identified, including (ADAM6, HK2, KRT19, LAMB3, LAMC2, PLEK2, S100A11 and S100A16), all of which exhibited a strong positive correlation with the DAMPscore across the three datasets (Fig. 5G and H).

Exploration of the function of PLEK2 in PAAD

Among the eight core genes, LAMB3, LAMC2 and PLEK2 showed the highest correlations with the DAMP score (Fig. 5H). While the roles of LAMB3 and LAMC2 in pancreatic cancer have been extensively studied (30–32), research regarding the functional involvement of PLEK2 in this context remains limited.

Initially, the expression levels of PLEK2 in PAAD were assessed. Transcriptional analyses across three independent datasets indicated that PLEK2 expression was significantly higher in tumor tissues compared to adjacent normal tissues (Fig. 6A). Furthermore, by integrating the mRNA data and proteomic data from the CPTAC, a strong positive correlation between PLEK2 mRNA levels and protein expression in pancreatic cancer was found (R=0.71, P<2.2×10−16; Fig. 6B).

Figure 6. Oncogenic function of PLEK2 in PAAD. (A) Transcriptional analyses of PLEK2 in tumor and normal tissues in the GSE62165, GSE71729 and Meta-cohort datasets. (B) Correlation analysis between the mRNA and protein level of PLEK2 in PAAD. (C) Immunohistochemical staining result of PLEK2 in pancreas and PAAD tissues deposited in the Human Protein Atlas. (D-F) Kaplan-Meier curves of overall survival in the PLEK2-stratified subgroups in the (D) TCGA-PAAD, (E) GSE85916 and (F) E-MTAB-6134 cohorts. (G-I) Enrichment analyses of PLEK2 stratified subgroups. (J) Quantitative PCR analyses of PLEK2 in AsPC-1 and PANC-1 cells after knockdown of PLEK2. (K and L) Cell Counting Kit-8 analyses of (K) AsPC-1 and (L) PANC-1 cells after knockdown of PLEK2. Transwell (M) invasion and (N) migration assay of AsPC-1 and PANC-1 cells after knockdown of PLEK2 (magnification, ×200). **P<0.01; ***P<0.001; ****P<0.0001. Adj., adjusted; TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma; Ctrl, control; si-PLEK2, small interfering RNA targeting PLEK2; OD, optical density.

Immunohistochemical analysis from the HPA revealed that PLEK2 is predominantly expressed in islet tissues of normal pancreatic tissue, while it is nearly undetectable in exocrine gland tissues (Fig. 6C). By contrast, PLEK2 expression was significantly increased in pancreatic cancer tissues (Fig. 6C). Using the median expression level of PLEK2 as a cutoff, it was observed that higher levels of PLEK2 expression were associated with poorer prognosis across three independent PAAD datasets (Fig. 6D-F).

Enrichment analysis indicated that tissues with elevated PLEK2 expression were primarily enriched in hallmark signaling pathways, including interferon alpha response (Fig. 6G), G2M checkpoint (Fig. 6H) and glycolysis (Fig. 6I). The interferon alpha response is closely related to immune cell infiltration (33), suggesting that PLEK2 may play a role in DAMP release. The G2M checkpoint is associated with cell proliferation, while glycolysis significantly impacts tumor growth, invasion and treatment resistance (34).

Subsequently, PLEK2 expression was knocked down in AsPC-1 and PANC-1 cells, and RT-qPCR confirmed the success of the interference of PLEK2 expression (Fig. 6J). CCK-8 assays indicated that silencing PLEK2 in PAAD cells promoted cell proliferation (Fig. 6L). In addition, Transwell assays demonstrated that knockdown of PLEK2 significantly impaired the migration and invasion capabilities of PAAD cells (Fig. 6M and N).

Discussion

The role of DAMPs in the context of PAAD has garnered increasing attention in recent years due to their potential to influence tumor behavior and patient outcomes (7,9). The present study aimed to evaluate the prognostic value of DAMP-related scores derived from a cohort of 40 DAMP-associated genes (21,22). Through meticulous bioinformatics analyses, patients with PAAD were stratified into three distinct clusters using Consensus Clustering: C1, C2 and C3. This classification underscores the dichotomy within the tumor microenvironment, wherein C1 is hypothesized to represent a pro-DAMP status, C3 an anti-DAMP status and C2 embodying an intermediate state. Notably, the C3 cluster was defined by markedly elevated expression of NT5E - a gene linked to impaired ICD - coupled with suppressed expression of key ICD-promoting genes, including TLR4, P2RX7, and FPR1 (Fig. 1G). Conversely, the C1 cluster demonstrated the highest transcriptional activity of genes associated with ICD facilitation, such as IL-33, FPR1 and TLR7. Mechanistically, NT5E (in concert with CD39) catalyzes the conversion of extracellular ATP into adenosine, a metabolite known to suppress antitumor immunity through immunosuppressive pathways (35,36). Furthermore, the release of HMG box 1 (HMGB1) - a hallmark of ICD - requires sequential permeabilization of nuclear and plasma membranes, representing a terminal post-mortem event. Extracellular HMGB1 acts as a potent pro-inflammatory mediator by engaging immune cell surface receptors, including TLR4 (37). Similarly, P2RX7 contributes to ICD induction via the purinergic receptor-inflammasome-IL-1 signaling axis (21,38). These findings collectively underscore the functional divergence of TME clusters in modulating immunogenic pathways and antitumor responses. This stratification not only reflects the complexity of DAMP biology in PAAD but also raises important implications for personalized therapeutic strategies directed towards modulating the immune landscape.

To further dissect the biological underpinnings of these clusters, differential gene expression analyses between the C1 and C3 groups were conducted. The results revealed that 141 DEGs had significant prognostic implications across both TCGA-PAAD and Meta-cohort datasets. This extensive overlap in prognostic genes indicates the potential of using DAMP-related scores to predict clinical outcomes in patients with PAAD. Employing LASSO Cox regression helped hone in on 6 pivotal genes, SLC16A3, LAMC2, IGF2BP2, HMGA2, EREG and ADM, which were able to differentiate between clusters with contrasting DAMP statuses. The subsequent computation of DAMPscores, informed by these genes, paves the way for a refined prognostic tool with clinical applicability.

The ability of the DAMPscore to predict OS was corroborated through robust ROC curve analyses. Furthermore, its superior performance in multiple datasets reinforces the necessity for larger, multi-center studies to validate its clinical utility. In a comparative analysis of the DAMPscore with 25 existing prognostic models, it was observed that the DAMPscore exhibited improved performance across most cohorts. This positions it as a potentially valuable addition to the existing arsenal of prognostic assessments in PAAD. Other studies have also found that prognosis models constructed based on DAMP or ICD-related genes can effectively guide the prognosis and personalized treatment of malignant tumors such as melanoma (29,39–41). The present study further expands the application scope of prognosis models developed based on DAMP-related genes.

Delving deeper into the biological relevance of the identified hub genes through WGCNA, it was noted that several of these genes, such as ADAM6, HK2 and S100, showed a positive correlation with the DAMPscore across datasets, suggesting their integral role in tumor biology and DAMP modulation. In particular, the role of PLEK2 warrants further attention. PLEK2 is a member of the pleckstrin family of proteins that primarily plays roles in cell signaling and cytoskeletal dynamics. Several studies have found that PLEK2 is emerging as an important player in tumorigenesis, with roles in promoting cell proliferation (42), enhancing metastatic potential (43) and facilitating immune escape (44) within cancers. Enrichment analyses revealed that PAAD tissues overexpressing PLEK2 were closely linked to hallmark signaling pathways, including the interferon alpha response, G2M checkpoint and glycolysis. The Interferon alpha response is particularly intriguing, as it serves as an indicator of immune cell infiltration and participates in the occurrence of ICD (45,46), hinting that PLEK2 may be a significant player in DAMP release and the subsequent inflammatory response within the tumor microenvironment. In advanced investigations, PLEK2′s functional relevance was confirmed through knockdown experiments in AsPC-1 and PANC-1 PAAD cell lines. The findings revealed that silencing PLEK2 expression promoted cell proliferation, impairing migration and invasion capabilities. This dual role in cell proliferation and invasive behavior reinforces its potential as a therapeutic target within the DAMP framework.

This study presents several noteworthy findings; however, it is essential to also acknowledge its limitations. First, the reliance on retrospective data from the TCGA-PAAD and Meta-cohort datasets may introduce biases inherent to such datasets, including selection bias and confounding factors that could affect the generalizability of the present results. In addition, while the DAMPscore demonstrated promising prognostic capabilities, its performance varied across different cohorts, indicating that external validation in larger, independent populations is necessary to confirm its robustness. Methodologically, this can be achieved by quantifying expression levels of DAMPscore-associated genes in patient tumor specimens via RNA sequencing or RT-qPCR, followed by correlation analyses with longitudinal prognostic data. Of note, the clinical implementation of multigene assays such as the 21-gene recurrence score (Oncotype DX®) in breast cancer (47) - which informs prognosis and guides therapeutic decision-making-provides a validated paradigm. This precedent offers critical methodological and translational insights for further investigating the clinical utility of the DAMPscore in PAAD. Furthermore, the functional exploration of PLEK2, although revealing significant insights into its role in PAAD, was limited to in vitro experiments, which may not fully resemble the complex tumor microenvironment in vivo. Future studies should aim to address these limitations by incorporating prospective data collection, validating the DAMPscore across diverse populations and exploring the functional implications of PLEK2 in more comprehensive in vivo models.

In conclusion, the present study validated the prognostic value of the DAMPscore in PAAD and elucidated its relationship with underlying molecular mechanisms driving tumor behavior. The framework established not only advances the understanding of DAMP biology in cancer but also unlocks possibilities for novel therapeutic strategies aimed at modulating the tumor immune response. The interplay between DAMPs, immune status and tumor metabolism represents a fertile ground for future investigations that may ultimately lead to enhanced therapeutic approaches in the challenging landscape of pancreatic cancer.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

Funding: No funding was received.

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

The study was designed by WD and JZ. In vitro experiments were conducted by HW, GL and JZ. Data analysis was carried out by WD, JL and ZS. Bioinformatics analysis was conducted by JL and ZS. WD and JZ confirm the authenticity of all raw data. The manuscript was drafted by WD, JL and HW, and was revised by all authors All authors have read and approved the final version of the study.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References