Transcriptome files and clinicopathological data of patients with PC were obtained from the Cancer Genome Atlas Program (TCGA) (https://tcga-data.nci.nih.gov/tcga/), Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), and the International Cancer Genome Consortium (ICGC) (https://www.icgc.org) databases. After removal of patients without complete clinical information and outcome status, as well as follow-up of fewer than 30 days, 176 PC patients from the TCGA pancreatic adenocarcinoma (TCGA-PAAD) cohort were retained as a training set. Ensembl was converted into gene symbol, and median value was kept when a genes had multiple gene symbols. The validation set contains 63 samples from the GSE57495 cohort and 215 patients of the PACA-CA cohort from the ICGC database. When multiple gene symbols appear or multiple probes appear for a gene, the median is taken as the gene expression value. A total of 134 human NKGs were downloaded from the ImmPort (https://www.immport.org/resource) database.

The prognostic NKGs were identified via univariate Cox regression analysis and were used to perform consensus clustering of PC patients. Consensus clustering analysis was conducted using the “ConsensusClusterPlus” R package to determine subgroups of PC patients based on the prognostic NKGs (Wilkerson and Hayes, 2010). The best classification was determined using the partition around medoids (PAM) algorithm and 1-Pearson correlation distance, with 500 bootstraps.

The DEGs among NKGs-derived clusters were screened out using “limma” package according to the false discovery rate (FDR) < 0.05 and |log2 [fold change (FC)]| > log2 (2) (Ritchie et al., 2015). The univariate and the least absolute shrinkage and selection operator (LASSO) Cox regression analysis were adopted to identify and filter prognosis-related NKGs, respectively. Finally, by choosing the optimal penalty parameter lambda correlated with the minimum 10-fold cross-validation, multivariate Cox regression analysis was then implemented to establish the prognostic signature. The formula for the risk signature was as follows: risk score = ∑ β i × E x p i . Where the βi represents the coefficient and Expi represents the normalized expression level of a gene. Two risk groups (high and low) were generated by a threshold of zero, and K–M analysis was conducted to compare overall survival (OS) differences between the high- and low-risk groups. The receiver operating characteristic (ROC) analysis was performed to estimate the predictive accuracy of the risk score.

GSEA was performed to analyze the differences in specific gene sets using the “GSVA” R package (Hänzelmann et al., 2013). The hallmark gene sets from the Molecular Signatures Database (MSigDB), the inflammation-related gene sets (Liu et al., 2020), and the angiogenesis-related gene set (Masiero et al., 2013) were used to be analyzed. These pathways with the FDR <0.05 was considered to be significant. Functional enrichment analysis included Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) (biological process (BP), cellular component (CC), and molecular function (MF)) analysis was performed on DEGs in clusters using WebGestaltR package (Yu et al., 2012).

The relative proportion of immune cells was calculated using the CIBERSORT algorithm (https://cibersort.stanford.edu/), which performs cell type enrichment analysis from gene expression data for 22 immune cells. The “ESTIMATE” R package was applied to estimate and extrapolate the fraction of stromal and immune cells in tumor samples (Yoshihara et al., 2013). The expression levels of the immune checkpoints were compared in different groups. To predict the chemosensitivity of osteosarcoma patients to several common anti-cancer drugs (methotrexate, paclitaxel, cisplatin, and doxorubicin), we adopted the “pRRophetic” R package to infer the half-maximal inhibitory concentration (IC₅₀) values by constructing the ridge regression model based on Genomics of Drug Sensitivity in Cancer (GDSC) (www.cancerrxgene.org/) cell line expression spectrum and gene expression profiles (Geeleher et al., 2014).

The decision tree model was applied to classify subgroups based on clinicopathogicial features and risk scores by using the “rpart” R package (https://cran.r-project.org/web/packages/rpart/index.html). The independent prognostic factors of OS for PC were identified by univariate and multivariate Cox regression analysis. A nomogram integrating the risk signature and independent prognostic clinicopathological factors was constructed in the TCGA cohort by the “rms” R package (https://cran.r-project.org/web/packages/rms/index.html). The calibration curves were utilized to evaluate the prediction accuracy between the predicted 1-, 2- and 3-year OS probabilities and the actual observations. The discriminate ability of the nomogram was assessed by time-dependent ROC curves. The decision curve analysis (DCA) was conducted to test the clinical utility of the nomogram using the “rmda” R package (https://cran.r-project.org/web/packages/rmda/index.html).

Tumor mutation burden (TMB) is was determined as the number of somatic indels and base substitutions per million bases in the coding region of the genome detected. Gene mutation data of PC patients were downloaded from the TCGA database and TMB was calculated using the “maftools” package (Mayakonda et al., 2018) as previously described (Chalmers et al., 2017).

12 PCD patterns were acquired from previous (Zou et al., 2022). Altogether, 580 apoptosis genes, 52 pyroptosis genes, 87 ferroptosis genes, 367 autophagy genes, 14 cuproptosis genes, 9 parthanatos genes, 15 entotic cell death genes, 101 necroptosis genes, 8 netotic cell death genes, 7 alkaliptosis genes, 220 lysosome-dependent cell death genes, and 5 oxeiptosis genes were collected. Based on the expression data of above gene sets, ssGSEA analysis was conducted on tumor samples using the R package GSVA.

The R software (v3.6.3) was used for statistical analyses. Wilcoxon test compared differences between two groups. Survival differences were compared using K–M curves with a Log-rank test. The Cox proportional hazard model was performed to estimate the β regression coefficient, hazard ratios, p-value, and their corresponding 95% confidence interval for each of the selected risk predictors. a nomogram was constructed with the “rms” package in R. The C-index and calibration curve with the bootstrap method were used to evaluate the prediction performance of the nomogram. A p-value <0.05 was deemed to be a statistical significance.

The flowchart is shown in Supplementary Figure S1. To obtain molecular subtypes of PC based on NKG, we first identified 32 NKGs that were significantly associated with the prognosis of PC (p < 0.05, Figure 1A). Notably, positive correlations among the expression of the 32 NKGs were observed in Figure 1B. Subsequently, consensus clustering of the 32 NKGs generated three stable clusters (C1, C2, and C3) in the TCGA-PAAD cohort (Figures 1C–E). Survival analysis demonstrated that the C3 cluster had a favorable prognosis whereas the C1 cluster had a poorer prognosis (Figure 1F). The individuals in the PACA-CA cohort were also divided into three clusters, which exerted similar prognosis characteristics as the clusters in the TCGA-PAAD cohort (Figure 1G). Among the 32 NKGs, the risk genes were generally overexpressed in the C1 cluster, and the protective genes were mainly elevated in the C3 clusters (Figure 1H).

We compared defined three clusters with the molecular subtypes derived from a pan-cancer study and immune signatures (Thorsson et al., 2018). As shown in Figure 2A, the C1 cluster presented with a higher TMB, aneuploidy, homologous recombination defects, and loss of heterozygosity. Meanwhile, a significantly higher proportion of immune signature-derived C3 subtype in our defined C3 subtype was observed (Figure 2B). The immune signature-derived C3 subtype was characterized by the overexpression of TH17 and Th1 genes, a low to moderate proliferation rate of tumor cells, and lower levels of aneuploidy and overall somatic copy number alterations. Meanwhile, the immune signature-derived C3 subtype showed a better prognosis than other subtypes, which is consistent with our defined C3 cluster showing the best prognosis, as shown in Figure 1F. The gene mutations in each cluster were compared and the top 20 genes with a lower p-value were illustrated in Figure 2C. Most mutations were present in KRAS, TP53, and SMAD4, accounting for 75.3%, 28.2%, and 19.7%, respectively. It was noticed that the C1 cluster with a poor prognosis had more gene mutations.

GSEA was performed to elucidate the pathway features in each cluster by using the Hallmark candidate gene sets. As shown in Figure 3A, the C1 cluster was significantly enriched in 38 pathways in the TCGA cohort. Generally, the activated pathways mainly included cell cycle-related pathways, such as E2F_TARGETS, G2M_CHECKPOINT, MYC_TARGETS_V1, whereas the inhibited pathways primarily contained INFLAMMATORY_RESPONSE, COMPLEMENT, and INTERFERON_GAMMA_RESPONSE. Similar results were also observed in the PACA-CA cohort. In addition, we compared the pathway characteristics between clusters (Figures 3B–D). It revealed that PC patients with the 3 subtype had activated immune pathways, such as cell cycle-related pathways, indicating that the 32 NKGs might play vital roles in the regulation of cell cycle and TME.

Furthermore, we assessed the relative abundance of 22 immune cells in the TCGA-PAAD and PACA-CA cohorts using the CIBERSORT algorithm. As shown in Figures 4A, C, significant differences among three clusters were observed for several immune cell types, such as CD8+T cells and activated CD4⁺ memory T cells. Meanwhile, we observed a significantly higher immune score in the C3 cluster than in other clusters (Figures 4B, D), indicating that the C3 cluster had a higher immune infiltration. In addition, we investigated the 7 inflammation-related metagenes clusters in the three molecular subtypes. As a result, 6 of the 7 metagenes clusters were significantly differently expressed among subtypes, except for interferon (Figure 4E). Overall, the C1 cluster presented with a higher inflammation activity than other clusters. Meanwhile, we also observed a higher enrichment score of LCK and MHC-II, and STAT1 in the C1 cluster than the other two clusters in the PACA-CA cohort (Figure 4F). The ssGSEA analysis of 12 PCD patterns indicated that 10 PCD patterns had obviously differences among 3 subtypes, and in general, C1 or C3 subtype had higher ssGSEA scores (Figure 4G).

Immunotherapy achieved favorable therapeutic effects in various cancers and immune checkpoint genes (ICG) play vital roles in these processes. Therefore, we evaluated the expression of ICGs among clusters and found an elevated expression of PD-1, PD-L1, and CTLA4 in the C3 cluster, as shown in Figure 5A. Meanwhile, we assessed the capability of clusters in predicting immunotherapy response using the T cell inflamed GEP score and observed a higher score in the C3 cluster than in other clusters (Figure 5B). INF-γ is a cytokine that plays a key role in immune regulation and anticancer immunity (Zhang et al., 2017), therefore, we calculated the ssGSEA score of the GOBP_RESPONSE_TO_INTERFERON_GAMMA gene set and found a significantly higher score of INF-γ response in the C3 cluster (Figure 5C). In addition, we also observed a higher CYT score in the C3 cluster than in other clusters (Figure 5D), which was used to reflect cytotoxic effects. Moreover, our data showed that the C1 cluster was more sensitive to cisplatin, gemcitabine, and erlotinib.

A total of 294 DEGs among clusters were identified, as shown in Figures 6A–C. Enrichment analysis on the DEGs was performed and the results showed that the C3 cluster contained DEGs that were significantly associated with immune-related pathways (Figure 6D). Univariate COX analysis showed that 122 of the 293 DEGs were significantly associated with the prognosis of PC (p < 0.01), including 84 risk genes and 38 protective genes (Figure 7A). Subsequently, lasso COX regression was adopted to compress the gene number and found 9 candidate genes when lambda = 0.0666 (Figures 7B, C). Finally, four genes were identified after stepwise multivariate regression analysis on the 9 candidate genes and were used to construct a prognosis model (Figure 7D), RiskScore = +0.306*MET+0.299*EMP1-0.225*NGFR+0.182*MYEOV. The four gene expressions were negatively correlated with methylation level (Supplementary Figure S2). The risk score was calculated for each patient in the TCGA cohort and was used to divided the patient into the high and low group (Figure 8A). ROC analysis demonstrated a favorable predictive capability in forecasting the 1-, 3-, and 5-year survival rates (Figure 8B). Survival analysis showed a significantly difference in prognosis between the high and low groups (Figure 8C). In addition, we evaluated the robustness of the prognosis model in the PACA-CA and GSE57495 cohorts, which had similar results as the TCGA cohort (Figures 8D–G).

The correlations between risk score and clinicopathological characteristics were analyzed in the TCGA and PACA-CA cohorts, and the results found significant associations between risk score and grade, but not stage, age, and gender (Figures 9A, C). Meanwhile, the risk score was significantly different among the three clusters, which manifested by a higher risk score in the C1 cluster and a lower risk score in the C3 cluster (Figures 9B, D). In addition, K-M curves showed that the risk score exhibited a favorable capability in the prognostic prediction of PC in sub-populations with specific clinicopathological features (Figures 9E, F).

As shown in Figure 10A, we observed a significantly difference in the relative abundance of four immune cells, including naive B cells, CD8 T cells, monocytes, and M0 macrophages, between the high and low groups in the TCGA cohort. The correlations between risk score and immune cells were illustrated in Figure 10B. In addition, a higher immune score was observed in the low group than the high group, indicating a higher immune infiltration in the low group (Figure 10C). The ssGSEA scores on each pathway were calculated for individuals and were compared between two risk groups. The results demonstrated that the High group was significantly associated with cell cycle-related pathways (Figures 10D, E).

As shown in Figure 11A, we observed a significantly higher T cell inflamed GEP score in the Low group as compared with those in the High group. Our data also revealed a significantly higher response to IFN-γ and cytolytic activity in the Low group, when compared with the high group (Figures 11B, C). In addition, we found elevated expression of PD-1 and CTLA4, but not PD-L1, in the low group (Figure 11D), suggesting potential differences in immunotherapy response between the two risk groups. Chemotherapy sensitivity in different risk groups was analyzed and found that the patients in the high group were more likely to be sensitive to gemcitabine, cisplatin, and erlotinib, as shown in Figure 11E. In addition, four of 12 PCD patterns had increased ssGSEA score in low group, while 3 PCD had higher ssGSEA score in high group (Figure 11F, left). Furthermore, we analyzed the correlation between RiskScore, four model genes and 12 PCD patterns, and there were different degrees of correlation with each other (Figure 11F, right).

As shown in Figure 12A, a decision tree was constructed based on the risk score and clinicopathological features and generated four groups (Lowest, Low, Mediate, High) using three parameters (risk score, N stage, age). Survival analysis demonstrated significant differences in prognosis among the four groups (Figure 12B, p < 0.001). The correlations between the decision tree-derived groups and risk groups were illustrated in Figures 12C, D. Univariate regression analysis showed that T stage, N stage, age, and risk score was associated with the prognosis of PC, and three of them (N stage, age, and risk score) were identified as independent risk factors via multivariate regression analysis (Figures 12E, F). Therefore, a nomogram was built using the three factors (Figure 12G). It was observed that the predicted values were close to the observed values in terms of the 1-, 2, and 3-year OS (Figure 12H), indicating that the nomogram had good prediction performance. In addition, a decision curve was used to evaluate the reliability of the model, and it was observed that the risk signature and nomogram model had a higher standardized net benefit as compared with other clinicopathological features (Figure 12I).