This study retrospectively collected transcriptomic data on pancreatic cancer from public databases after removing normal tissue specimens and specimens without clinical follow-up data, including the Cancer Genome Atlas (TCGA; https://www.cancer.gov/tcga; n = 177) as well as Pancreatic Cancer-Australia (PACA-AU; n = 91) and Pancreatic Cancer-Canada (PACA-CA; n = 234) from the International Cancer Genome Consortium (ICGC; https://www.icgc-argo.org) (Zhang et al., 2019). Due to samples from different platforms, the batch effects were removed utilizing the ComBat function of the sva package (version 3.42.0) (Leek et al., 2012). A principal component analysis (PCA) was conducted to evaluate the data before and after the removal of the batch effects. The follow-up data and clinicopathological characteristics were also collected. Moreover, single-nucleotide variant (SNV) and copy number variation (CNV) data on pancreatic cancer were retrieved from TCGA project. After reviewing the previously published literature, we collected 159 necroptosis genes, as listed in Supplementary Table S1. The GSE21501 cohort containing expression profiling and follow-up information of 101 pancreatic cancer patients was downloaded from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/), which was used as the external validation cohort (Stratford et al., 2010; Stratford et al., 2017). Supplementary Figure S1 depicted the workflow of our study.

Univariate Cox regression models were conducted between necroptosis genes and pancreatic cancer survival, and genes with p < 0.05 were determined for a consensus clustering analysis. A consensus clustering approach offers quantitative and visual stability evidence to estimate the number of unsupervised classes within a specified data set. The ConsensusClusterPlus package (version 1.58.0) adopts the consensus clustering approach, comprising consensus matrix, empirical cumulative distribution function (CDF), and delta area plots (Wilkerson and Hayes, 2010). Through implementing the consensus clustering analysis, necroptosis subtypes were clustered based on the expression values of prognostic necroptosis genes across pancreatic cancer specimens. The number of clusters k was set as 2–9, and 80% of the samples were sampled using a re-sampling method. After multiple sampling, stable and reliable unsupervised classes were found in line with the following parameters: re-samplings = 50, proportion of items to sample = 0.8, proportion of features to sample = 1, and distance = “pearson”. PCA was conducted to visualize the difference in expression levels of prognostic necroptosis genes among diverse necroptosis subtypes.

GSVA, a non-parametric and unsupervised gene set enrichment approach, can estimate the enrichment score of specific pathways or signatures in accordance with transcriptomic profiles (Hänzelmann et al., 2013). The 50 hallmarks of gene sets were retrieved from the Molecular Signatures Database (Liberzon et al., 2015). The activity of each hallmark pathway was quantified using the single-sample gene set enrichment analysis (ssGSEA) function.

The relative infiltrations of immune cell populations were estimated with the ssGSEA function derived from the GSVA package (Hänzelmann et al., 2013) on the basis of the expression values of 782 meta-genes (Charoentong et al., 2017) in pancreatic cancer specimens. The mRNA expressions of known immune checkpoints and human leukocyte antigen (HLA) molecules were quantified in each pancreatic cancer specimen.

Chen and Mellman proposed a cancer–immunity cycle to evaluate antitumor immune responses, containing seven steps: 1) release of cancer antigens, 2) cancer antigen presentation, 3) priming and activation, 4) trafficking of T cells to tumors, 5) infiltration of T cells into tumors, 6) recognition of cancer cells by T cells, and 7) killing of cancer cells (Karasaki et al., 2017). The levels of each step within the cancer–immunity cycle were quantified using the ssGSEA approach.

The gene sets of known biological processes were retrieved from Mariathasan et al. (2018), containing epithelial–mesenchymal transition (EMT1-3), immune checkpoint, antigen processing machinery, CD8 T effector, angiogenesis, pan-fibroblast TGFβ response (pan-F-TBRS), DNA damage repair, FGFR3-related genes, KEGG-discovered histones, Fanconi anemia, cell cycle, cell cycle regulators, DNA replication, nucleotide excision repair, homologous recombination, mismatch repair, and WNT target. The enrichment score of aforementioned biological processes was quantified using the ssGSEA approach.

The therapeutic responses to known chemotherapy agents (cisplatin, gemcitabine, docetaxel, and paclitaxel) were estimated using the pRRophetic package (Geeleher et al., 2014). Through construction of the ridge regression model on the basis of the Genomics of Drug Sensitivity in Cancer (GDSC) pharmacogenomics database (www.cancerRxgene.org) (Yang et al., 2013) and transcriptomic data, the half-maximal inhibitory concentration (IC50) of each chemotherapeutic agent was calculated across pancreatic cancer specimens.

Utilizing the maftools package (version 2.10.0) (Mayakonda et al., 2018), SNV data were analyzed and visualized on the basis of the mutation annotation format (MAF) of pancreatic cancer. GISTIC2.0 (Mermel et al., 2011) was implemented to analyze copy number amplification and deletion.

A differential expression analysis was implemented between any two necroptosis subtypes utilizing linear models for the microarray data (limma; version 3.50.0) package (Ritchie et al., 2015). Genes with adjusted p-value<0.05 were screened, and necroptosis-relevant genes were determined following the intersection of differentially expressed genes.

Utilizing the clusterProfiler package (version 4.2.0) (Yu et al., 2012), a functional enrichment analysis of necroptosis-relevant genes was carried out, comprising Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses.

Prognostic necroptosis–relevant genes with p < 0.05 were determined through univariate Cox regression models, which were ranked by using the randomForestSRC package (version 2.14.0), following number of replication = 100, number of step = 5, Monte Carlo iteration number = 100, and genes with relative importance ˃ 0.4. Thereafter, a necroptosis-relevant risk model was generated on the basis of the expression of the most important genes and regression coefficients from a multivariate Cox regression model. Following calculation of the necroptosis-relevant risk score of each pancreatic cancer patient, high- and low-risk groups were separated in accordance with the median value of risk score. Receiver operator characteristic (ROC) curves at 1-, 3-, and 5-year survival were conducted to evaluate the predictive reliability of the necroptosis-relevant risk model in pancreatic cancer survival. Using the GEPIA web tool (http://gepia.cancer-pku.cn/), the expression of genes in the necroptosis-relevant risk model was examined in pancreatic cancer (n = 179) and normal tissues (n = 171).

In the nomogram, the line length indicates the degree of influence of a specific variable and diverse values of this variable on outcomes. After univariate and multivariate Cox regression models, a nomogram was generated on the basis of independent prognostic factors through the rms package (version 6.2-0), showing the intuitive and effective results of the risk model. Calibration curves were utilized to validate the predictive accuracy of the nomogram-predicted survival probabilities for 1-, 3-, and 5-year survival.

R software (version 3.6.1) was implemented for data processing. Univariate and multivariate Cox regression analyses were conducted, and the hazards ratio (HR) and p-value were calculated to evaluate the correlations of variables with pancreatic cancer survival. Kaplan–Meier curves and log-rank test were depicted for the survival difference between groups. The difference between two groups was compared with student’s t-test or Wilcoxon test, while comparison between three groups was presented through the Kruskal–Wallis test. A correlation analysis was carried out via Pearson’s or Spearman’s test. The C-index was calculated for estimating the prediction performance through the survival package. p < 0.05 indicated statistical significance.

We retrospectively collected transcriptomic data on pancreatic cancer from TCGA, PACA-AU, and PACA-CA cohorts. The batch effects of integrated data were eliminated for subsequent analyses, which were visualized through PCA (Figures 1A,B). Among 159 necroptosis genes, 18 genes (SPATA2, AIFM1, SLC25A4, BCL2, SPATA2L, TYK2, SMPD1, STAT5B, SLC25A6, USP21, STAT4, VPS4A, RIPK1, PLA2G4C, IL33, CAMK2B, MAPK10, and BAX) were protective factors of pancreatic cancer prognosis, while 14 genes (TNFRSF10B, HSP90AA1, BIRC3, TNFRSF10A, CHMP4C, CASP8, FADD, CAPN2, GLUD1, PYGL, BIRC2, CAPN1, CHMP2B, and IFNA13) were risk factors of prognosis, as depicted in Figure 1C. These prognostic necroptosis genes were utilized for the consensus clustering analysis. When k = 5, pancreatic cancer samples were clearly separated into five clusters (Figure 1D). Figure 1E depicted the CDF when k takes different values, and we found that when k = 5, CDF reached the approximate maximum, indicative of cluster stability. Figure 1F showed the relative change in CDF of k compared to k-1. When k = 6, CDF only slightly decreased, so 5 was the appropriate value of k. Ultimately, five necroptosis subtypes were identified for pancreatic cancer, namely, C1 (n = 187), C2 (n = 135), C3 (n = 121), C4 (n = 21), and C5 (n = 38). PCA also confirmed the reliability of necroptosis subtypes (Figure 1G). The survival analysis demonstrated the remarkable survival difference among necroptosis subtypes (Figure 1H). The C1 subtype had the worst survival outcomes, followed by C3, C2, C4, and C5. Figure 1I depicted the prominent expression difference of prognostic necroptosis genes among diverse subtypes. The accuracy and reliability of necroptosis subtypes were confirmed in the TCGA cohort (Supplementary Figures S2A-F).

Further analysis was conducted to uncover the mechanisms underlying five necroptosis subtypes. In Figure 2A, tumorigenic pathways (hedgehog signaling, KRAS, angiogenesis, glycolysis, etc.) were remarkably activated in C1 and C2 subtypes, contributing to an undesirable prognosis. C4 and C5 subtypes presented the relatively high infiltrations of immune cells, while C2 was characterized by low infiltrations of immune cells (Figure 2B). Most immune checkpoints were markedly downregulated in C4 and C5 subtypes, while their upregulations were found in C2 (Figure 2C). Tumors can evade T-cell responses through losing the major histocompatibility complex (MHC)/HLA class I and II molecules (Godfrey et al., 2018). In Figure 2D, we observed the loss of HLA class I and II molecules in C4 and C5. Differently, C2 displayed the prominent activation of HLA molecules, followed by relatively modest expression in C1 and C3. C4 and C5 subtypes presented the relatively low levels of almost all steps within the cancer–immunity cycle in comparison to other subtypes; meanwhile, the C2 subtype had the highest activation of each step (Figure 2E). Similarly, CD8 T effector and antigen processing machinery, immune checkpoint, and stromal activation (EMT1-3) were relatively downregulated in C4 and C5 (Figure 2F); C2 had relatively high levels of immune and stromal activation pathways; and C1–3 presented the enhanced cell cycle progression (cell cycle, cell cycle regulators, DNA replication, etc.). Overall, five necroptosis subtypes had diverse immunogenic features.

We compared the sensitivity to known chemotherapy agents (cisplatin, gemcitabine, docetaxel, and paclitaxel) in five necroptosis subtypes. As depicted in Figure 3A, C4 had the highest IC50 values of cisplatin, gemcitabine, docetaxel, and paclitaxel, while C5 presented the lowest IC50 values of aforementioned chemotherapy agents, indicating that C4 presented the highest probability of chemotherapy resistance while C5 was the most sensitive to these chemotherapy agents. We also investigated that C1 and C3 had relatively higher tumor mutation burden (TMB) than other subtypes (Figure 3B). KRAS (53%) and TP53 (53%) were the most frequent mutant genes. The widespread copy number amplification (Figure 3C) and deletion (Figure 3D) occurred in pancreatic cancer specimens. Among five necroptosis subtypes, C2 and C4 had the relatively decreased fractions of genome altered (FGAs), as depicted in Figure 3E. Moreover, we noted the relatively lowered copy number amplification and deletion in C2 and C4 in comparison to other subtypes (Figure 3F). Aforementioned data uncovered the difference in tumor mutations among necroptosis subtypes.

Through the intersection of DEGs (adjusted p < 0.05) between any two subtypes, we determined 591 necroptosis-relevant genes (Supplementary Table S2). Their biological significance was further analyzed through GO and KEGG enrichment analyses. In Figure 4A, necroptosis-relevant genes were remarkably linked with regulation of protein localization to the membrane. Moreover, they had prominent associations with tumorigenic pathways (p53 signaling pathway and cell senescence) and immune pathways (Th17 cell differentiation, PD-L1 expression, PD-1 checkpoint pathway in cancer, Th1 and Th2 cell differentiation, etc.), indicative of the critical roles of necroptosis-relevant genes in pancreatic cancer progression (Figure 4B). Among all necroptosis-relevant genes, 207 displayed significant correlations to pancreatic cancer prognosis (Supplementary Table S3). Using the random forest approach, we determined the most important genes with the relative importance > 0.4 (Figures 4C,D). A multivariate Cox regression model was constructed in line with the following formula: risk score = 0.119399555 * MYEOV expression + (−0.258345687) * HDAC4 expression + 0.26238863 * TLDC1 expression + (−0.395042137) * PITPNA +0.175544976 * FNDC3B expression + 0.338675676 * HMGXB4 expression + (−0.150557275) * BAX expression. Following the calculation of the risk score, all patients were separated into high- and low-risk groups (Figure 4E). The high-risk group had more dead patients relative to the low-risk group (Figure 4F). The survival analysis demonstrated the survival advantage of low-risk patients (Figure 4G). The difference in expression of aforementioned genes between groups is visualized in Figure 4H.

Uni- and multivariate Cox regression analyses uncovered that age- and necroptosis-relevant risk models were both independently associated with pancreatic cancer survival (Figures 5A,B). AUCs at 1-, 3-, and 5-year survival were separately 0.714, 0.724, and 0.757, indicative of the reliability of the necroptosis-relevant risk model in predicting survival outcomes (Figure 5C). To facilitate the clinical application of the necroptosis-relevant risk model, we generated a nomogram following integration of age (Figure 5D). Calibration curves demonstrated the predictive accuracy of this nomogram in pancreatic cancer survival (Figures 5E–G). In addition, we also performed a stratified analysis and demonstrated that the risk model can serve as an independent prognostic factor without consideration of the impact of age (Supplementary Figures S3A,B).

The robustness of the necroptosis-relevant risk model was verified in the GSE21501 cohort. In accordance with the same formula, we computed the necroptosis-relevant risk score of each pancreatic cancer patient in the external cohort (Figure 6A). As expected, high-risk patients had poorer survival outcomes than low-risk patients (Figure 6B). AUCs at 1- and 3-year survival were separately 0.69 and 0.71 (Figure 6C), demonstrating the excellent performance in predicting prognosis. Compared with the existing prognostic models constructed by Chen et al. (2021), Xiao et al. (2022), and Zhang et al. (2022), the necroptosis-relevant risk model had a higher C-index (Figure 6D), indicating the advantage of this model in predicting prognosis. We also examined the expression of genes in the necroptosis-relevant risk model using the GEPIA web tool. BAX, FNDC3B, HDAC4, HMGXB4, MYEOV, and TLDC1 displayed upregulated expressions in pancreatic cancer than normal tissues (Figure 6E).

Compared with other necroptosis subtypes, C1 presented a relatively higher necroptosis-relevant risk score, followed by C3 (Figure 7A), indicating the heterogeneity in the risk score among diverse necroptosis subtypes. Further analysis was conducted to evaluate the correlations of necroptosis-relevant risk score with tumor immunogenicity. In Figure 7B, as the necroptosis-relevant risk score increased, the infiltrations of immune cells gradually decreased, indicative of the negative correlations of the necroptosis-relevant risk score with immune cell infiltrations. Moreover, we noted that the necroptosis-relevant risk score was negatively linked with the expression of immune checkpoints and HLA molecules (Figures 7C,D). The aforementioned data indicated the role of the necroptosis-relevant risk model in tumor immunogenicity of pancreatic cancer.

In Figure 8, we noted that the necroptosis-relevant risk score presented a significantly positive correlation to the release of cancer cell antigens but displayed significantly negative correlations to recruiting of B cells, CD4 T cells, dendritic cell, macrophages, T cells, Th17 cells, Treg cells, and killing of cancer cells, indicative of the remarkable interactions of the necroptosis-relevant risk score with the cancer–immunity cycle. Moreover, necroptosis-relevant risk score was negatively linked with CD8 T effector, and angiogenesis but was positively associated with pan-F-TBRS, FGFR3-related genes, EMT2, KEGG discovered histones, Fanconi anemia, cell cycle, cell cycle regulators, DNA replication, DNA damage repair, nucleotide excision repair, homologous recombination, and mismatch repair, indicative of the mechanisms underlying the necroptosis-relevant risk score.

Further analysis was conducted to evaluate the correlations between the necroptosis-relevant risk score and chemosensitivity of pancreatic cancer. No remarkable differences of the IC50 values of cisplatin and gemcitabine were noticed between high- and low-risk groups (Figures 9A,B). But the high-risk group presented the prominently reduced IC50 values of docetaxel and paclitaxel relative to the low-risk group (Figures 9C,D). This indicated that high-risk patients were more likely to be sensitive to docetaxel and paclitaxel.