We downloaded RNA-seq data as well as clinicopathological and mutation data of GC, from the TCGA database (https://portal.gdc.cancer.gov/), which included 375 GC tissues and 32 adjacent normal tissues. The TCGA cohort was used for model development, and the GSE84437 dataset that contained 433 GC samples was obtained as an external validation cohort. Patients with complete follow-up and clinical information were included in this study. The basic characteristics of the patients in the TCGA and GEO datasets are shown in Table 1. Additionally, we obtained 159 necroptosis-associated genes from previous literature (Supplementary Table S1) [12].

To identify prognosis-related NAGs in GC, we performed univariate Cox regression analysis to screen for prognostically relevant NAGs. The Least Absolute Shrinkage and Selection Operator (LASSO) was employed to avoid overfitting and get rid of those tightly correlated genes. Subsequently, candidate NAGs were employed to construct a prognostic signature. The risk score for each patient is calculated as follows: R i s k   s c o r e = r e g r e s s i o n   c o e f f i c i e n t ( g e n e i ) × e x p r e s s i o n   v a l u e ( g e n e i )

We stratified patients into high-risk groups and low-risk groups according to the median risk score, and the difference in OS, cancer-specific survival (CSS), and progression-free survival (PFS) between the high-risk and low-risk groups was assessed using the Kaplan-Meier (K-M) curves. Additionally, the risk score of each sample was reordered, and a risk curve and a survival status-related scatterplot were shown as a result. The area under the receiver operating characteristic (ROC) curve (AUC) was calculated to assess the accuracy of the model. Patterns of gene expression in the defined patient groups were examined via principal component analysis (PCA) with the “scatterplot3d” package in R.

The external cohort GSE84437 was employed to verify the necroptosis-related prognostic model. Stratifying patients into low- and high-risk groups were done using the same cut-off value of the training set. Then, Kaplan Meier, ROC, and PCA were performed.

We analyze the association of risk scores with clinicopathological traits, including age, gender, grade, and pathologic stage. Then, univariate and multivariate Cox regression analysis was performed to identify important predictive clinical variables for the development of genome-clinicopathologic nomogram to predict individual survival probability for GC patients. The nomogram was constructed by using R package “rms” and the accuracy of the nomogram was assessed via ROC and calibration curves. Moreover, we also investigated whether the risk score could affect the OS of patients in distinct clinical subgroups by log-rank test.

The infiltrating score of immune cells in the TIME was calculated by Single-sample Gene Set Enrichment Analysis (ssGSEA), including 16 types of infiltrating immune cells. Moreover, the “estimate” package was used to evaluate the infiltration of immune cells and stromal cells in tumor tissues and to infer the tumor purity.

The tumor immune dysfunction and exclusion (TIDE) score was calculated to predict inhibitory responses to PD-1 and CTLA4 immune checkpoint in GC patients [13]. The immunophenoscore (IPS) is calculated based on the expression of various important immune molecules, including immune regulators, MHC molecules, suppressor cells, and effector cells, and can well reflect the response rate to immune checkpoint inhibitors (ICI) [17]. Moreover, the “pRRophetic” R package was employed to calculate the predicted half-maximal inhibitory concentration (IC50) of commonly applied chemotherapy drugs for obtaining the drug sensitivities of GC patients.

There is growing evidence that patients with high TMB have an acceptable response to immunotherapy and are associated with a good prognosis. We evaluated the tumor somatic mutations presented in high- and low-risk patients separately using the “maftools” R package. We further performed TMB variance analysis and correlation analysis for different risk groups. Kaplan-Meier analysis was performed to compare survival differences between low and high TMB. Subsequently, risk score and TMB were combined to perform a survival analysis to determine if there are differences in patients between different groups. We performed a combined survival analysis of risk score and TMB to explore whether the prognostic value of risk score was influenced by TMB status. Moreover, the correlation between the risk score and MSI was explored.

To explore the signaling pathways that necroptosis-related signatures may be involved in regulation, differentially expressed genes (DEGs) between two risk subgroups were retrieved (adjusted p < 0.001 and |log₂FC| ≥ 2). Subsequently, Gene Ontology (GO) and functional annotation of Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed using “clusterProfiler” R package.

All statistical analyses were completed by R software 4.1.0. The Spearman correlation analyses were applied to determine the correlation. Kaplan-Meier method was used for survival analysis of different risk subgroups using a two-sided log-rank test. The Wilcoxon signed-rank test was employed to compare difference of the two groups. p < 0.05 was the threshold of significance.

In univariate Cox regression analysis, a total of 10 NAGs were evidently associated with OS (Figure 1A). To minimize overfitting, cross-validation was performed in the LASSO regression (Figures 1B,C). Finally, six NAGs were employed to construct the prognostic signature. Risk score = GLUD2* 0.0562 + MAPK10*0.1662 + CHMP4C*(-0.0033) + IFNA4*0.7025 + IFNA14* 0.2007 + IFNB1*0.0054. Patients were categorized into low- and high-risk subgroups and the Kaplan-Meier plot revealed that the OS, CSS, and PFS of low-risk patients was higher than those of high-risk patients (Figures 1D–F). ROC curves presented the AUCs were 0.704 (3 years) and 0.724 (5 years), which manifested that the prognostic signature had a satisfactory predictive efficiency (Figure 1G). Besides, the PCA demonstrated that the risk score could distinguish the two risk subgroups (Figure 1H). The scatter plot risk and score distribution plot suggested that the high-risk subgroup had higher risk scores and worse survival times than that of the low-risk subgroup (Figure 1I).

The predictive ability of the signature was well reproduced in an external validation set (GSE84437). The K-M curve revealed that patients in low-risk group had significantly favorable prognoses (p < 0.05, Figure 2A). The AUC value of ROC curve revealed that the risk score could predict the OS rate of GC patients to some extent (Figure 2B). PCA indicated that the patient distribution of the two risk groups was in two directions (Figure 2C). The risk curve and survival status-related scatterplot also demonstrated the same tendency of prognosis in two risk groups as the K-M curve (Figure 2D).

We investigated the association between risk score and clinicopathological parameters in the TCGA training cohort. The prognostic signature observed statistically obvious higher values in G3 than in G1-2 (p = 0.001; Figure 3A). Similarly, a higher risk score was observed in stage III-IV than in stage I-II (p = 0.02; Figure 3B). To further explore the prognostic role of this signature, we performed Cox regression analyses in the TCGA cohort to select independent prognostic indicators. The results suggested that the risk score (HR: 8.137, 95% CI: 3.577–18.511, p < 0.001), age (HR: 1.035, 95% CI: 1.017–1.053, p < 0.001) and stage (HR: 1.661, 95% CI: 1.328–2.076, p < 0.001) were independent prognostic factors for OS (Table 2). Then, a nomogram based on three identified prognostic variables was built in the TCGA cohort to predict 3- and 5-year OS rates, respectively (Figure 3C). The ROC curves were subsequently depicted in Figure 3D, with 3-, and 5-year AUCs of 0.72 and 0.71, respectively. Furthermore, favorable consistency between prediction based on nomogram and actual observed outcomes of 3- and 5-year OS rates were illustrated in the calibration plots, respectively (Figure 3E). Additionally, we evaluated the OS of patients in two risk groups among distinct clinical subgroups, including age, gender, tumor grade, and TNM stage. As is shown in Supplementary Figure S1, patients with high risk had worse survival probabilities than those of patients with low risk in all distinct clinical subgroups.

Based on the ssGSEA algorithm, most immune infiltration cells had higher scores in the high-risk group, including CD8⁺ T cells, B cells, DCs, immature DCs (iDCs), neutrophils, plasmacytoid dendritic cells (pDCs), T helper cells, mast cells, tumor-infiltrating lymphocytes (TIL), T follicular helper cells (Tfh), and regulatory T (Treg) cells (Figure 4A). Moreover, the score of type-Ⅱ interferon (IFN) response was higher in the high-risk group, but scores of MHC-class Ⅰ were higher in the low-risk group (Figure 4B). The ESTIMATE algorithm was employed to score the TIME. The violin plot shows the ratio of TIME scores (stromal, immune, and ESTIMATE scores) between two risk groups. The TIME score of the high-risk subgroup was higher than that of the low-risk subgroup (Figure 4C).

The TIDE algorithm was employed to identify GC patients who can benefit from immunotherapy. The results demonstrated that TIDE scores were significantly lower in the low-risk subgroup than in the high-risk subgroup, suggesting a better response to ICI immunotherapy (Figure 4D). Furthermore, the immunogenicity of two risk groups was analyzed by IPS analysis. The ips_ctla4_pos_pd1_neg, ips_ctla4_neg_pd1_pos, and ips_ctla4_pos_pd1_pos scores were higher in the low-risk group (Figures 4E–G), suggesting that low-risk patients have a better response for immunotherapy. These results suggest that the signature can predict the response of GC patients to immunotherapy.

To explore whether risk scores predict chemotherapy response in GC patients, we compared IC50 levels in two groups of chemotherapy drugs or inhibitors. The results are shown in Figure 5, that Dasatinib, Pazopanib, Axitinib, and Rapamycin may be potential to treat high-risk patients, while Gefitinib and Metformin may suitable for patients in the low-risk group (Figure 5).

There is growing evidence that patients with high TMB have an acceptable response to immunotherapy and are associated with a good prognosis. A lower somatic frequency (85.29%) was frequently observed in the high-risk subgroup compared to the low-risk subgroup (95.93%) (Figures 6A,B). TMB score was obviously lower in high-risk group (p < 0.001; Figure 6C) and risk score was negatively correlated with the TMB score (R = −0.43, p < 0.001; Figure 6D). As a TMB survival curve, we found that high TMB patients have a better prognosis (p = 0.03; Figure 6E). A combined analysis of risk score and TMB found that the GC patients with high TMB in high-risk group had significantly shorter OS relative to low-risk patients (p = 0.001; Figure 6F). And the results in patients with low TMB were consistent, but not statistically different (p = 0.196; Figure 6G). In addition, multivariate Cox regression analysis indicated that the risk score was independent of TMB in predicting OS (Table 2). Additionally, the MSI status in the two groups was also analyzed. As shown in Figure 6H, the proportion of MSI-H and MSI-L in the low-risk subgroup was higher than that in the high-risk subgroup. Risk scores were significantly different among the different microsatellite status groups (Figure 6I).

KEGG pathway analyses and GO annotation were performed based on 1671 differentially expressed genes (FDR <0.001 and |logFC| > 2) between two risk subgroups. The most enriched terms in GO are extracellular matrix organization, receptor ligand activity, and collagen-containing extracellular matrix (Figure 7A); while those in KEGG with higher enrichment include PI3K-Akt signaling pathway, neuroactive ligand-receptor interaction, calcium signaling pathway, and focal adhesion (Figure 7B).