Identification of immunogenic cell death-related prognostic genes in gastric cancer

Introduction

Gastric cancer (GC), ranking as the fifth most prevalent malignancy and the third leading cause of cancer-related mortality worldwide (1), arises from synergistic interactions between environmental carcinogens [e.g., Helicobacter pylori infection, Epstein-Barr virus (EBV), smoking, dietary nitrosamines] and genetic predisposition (2). Current therapeutic paradigms are stage-dependent: localized tumors undergo radical resection with perioperative chemotherapy, while advanced or metastatic disease is managed with platinum-fluoropyrimidine combinations, taxanes (docetaxel/paclitaxel), or irinotecan-based regimens (3,4). Nevertheless, >60% patients present with advanced-stage disease due to nonspecific symptomatology, culminating in a median overall survival (OS) <12 months (5). These unmet clinical needs underscore the imperative need to identify robust prognostic genes, refine early diagnostic modalities, and develop mechanism-driven therapies.

The GC tumor microenvironment (TME) forms a spatially organized, immunosuppressive niche composed of diverse cellular components [e.g., cytotoxic T cells, tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs)], stromal elements [e.g., cancer-associated fibroblasts (CAFs), extracellular matrix (ECM)], and soluble mediators [e.g., transforming growth factor-beta (TGF-β), interleukin-10 (IL-10)] (6,7). Key metabolic perturbations, including hypoxia and oncometabolite accumulation (e.g., lactate), drive immune cell plasticity towards tolerogenic states, thereby underpinning resistance to chemo- and immunotherapy (8-10). Immunogenic cell death (ICD), a regulated cell death modality triggered by chemotherapeutics (e.g., oxaliplatin) or physical stressors (e.g., photodynamic therapy), is characterized by spatiotemporal release of damage-associated molecular patterns (DAMPs), such as surface-exposed calreticulin (CRT), secreted high mobility group box 1 (HMGB1), and extracellular ATP (11). These DAMPs orchestrate dendritic cell (DC) maturation, antigen cross-presentation, and clonal expansion of tumor-specific CD8⁺ T cells, effectively transforming dying tumor cells into in situ vaccines that prime durable immunological memory (12-14). Emerging evidence highlights the dualistic roles of ICD in GC pathogenesis. While ICD induction via oxaliplatin or oncolytic viruses enhances tumor immunogenicity and synergizes with immune checkpoint inhibitors (ICIs) by activating antitumor immunity (15,16), paradoxically, key ICD-related mediators such as CRT can exert context-dependent pro-tumorigenic effects. Specifically, CRT promotes metastatic progression through G9a-mediated H3K9 dimethylation at the E-cadherin promoter, driving epithelial-mesenchymal transition (EMT) and facilitating GC dissemination (17). This functional dichotomy extends to microRNA (miRNA) regulation: miR-637 upregulation exacerbates endoplasmic reticulum (ER) stress-induced apoptosis by suppressing CRT, yet, simultaneously disrupts CRT-mediated prosurvival signaling, unveiling a delicate balance between ICD-driven immunogenicity and CRT-facilitated metastatic adaptation (18). Despite these advances, the prognostic utility of ICD-related biomarkers remains controversial, necessitating multi-omics integration to resolve heterogeneity in immune contexture and therapeutic responses (15).

Single-cell RNA sequencing (scRNA-seq) has gradually emerged as a non-biased cutting-edge technology. The core idea of the reported single-cell integration studies for prognosis or classification is to analyze the heterogeneous subpopulation characteristics of cancer cells, immune cells, and stromal cells in the tumor and its microenvironment, and to combine bulk transcriptome data to construct or optimize tumor classification and prognosis models, thereby revealing potential mechanisms related to tumor progression and treatment response (19-23). In contrast, this study, based on the screening of GC-specifically expressed cells, further combined weighted gene co-expression network analysis (WGCNA) to screen out module genes highly related to ICD, thereby obtaining prognostic genes. This design is distinct from previous single-cell integrated prognosis or stratification studies, and clearly highlights the unique contribution of this study-focusing on the specific biological process of ICD rather than generalized immune characteristics. Furthermore, the prognostic model was then constructed using prognostic genes to evaluate the survival probability of GC patients. The identification of these novel prognostic genes not only enhances the diagnostic accuracy but also provides potential therapeutic opportunities for GC patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-aw-831/rc).

Methods

Data source

The training set for this study was derived from The Cancer Genome Atlas (TCGA)-GC dataset, accessed through the UCSC Xena database (https://xenabrowser.net/datapages/). This dataset consisted of 410 GC and 36 control samples, with OS information available for 383 GC samples. Additional validation dataset (GSE62254) and scRNA-seq dataset (GSE183904) were both acquired from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE62254 dataset was sequenced by GPL570 platform, containing 300 GC tissue samples with survival information, while GSE183904 dataset, sequencing by GPL24676 platform, included 26 GC and 10 control tissue samples. A total of 32 ICD-related genes (ICDRGs) were collected from the published literature (24). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

ScRNA-seq analysis

The data in GSE183904 dataset were filtered using “Seurat” package (version 5.1.0) (25) with following standards: removing genes detected in less than 200 cells, nFeature_RNA ≥6,000, nCount_RNA ≥80,000, and cells with ≥20% mitochondria. After data normalization, vst method was utilized to extract top 2,000 high variable genes for following analysis. Principle component analysis (PCA) was performed to reduce dimension and select principle components (PCs). FindNeighbors and FindClusters functions in “Seurat” package were conducted to process unsupervised cluster analysis with a resolution of 0.1. Uniform manifold approximation and projection (UMAP) analysis was used to cluster cells and visualize the results. Among cell clusters, cell subpopulations were annotated relied on marker genes in the literature (26,27). Potential doublets in scRNA-seq data were identified and removed using “DoubletFinder” package (version 2.0.4) (27).

Cell communication and pseudo-time trajectories analyses of key cells

To investigate the correlation of cell subpopulations, cell communication analysis was processed by “CellChatDB” package (version 1.6.1) (28). Furthermore, the proportion of cell subpopulations was evaluated in GC and control samples. AUCell algorithm was used to calculate ICDRGs-related score of samples in GSE183904 dataset. The cells with the significant differential score between GC and control samples, and exhibiting the highest expression in GC samples, were identified as key cells. Based on obtained key cells, pseudo-time trajectories analysis was processed by “Monocle” package (version 2.26.0) (29), and the expression of ICDRGs in the key cells differentiation process was evaluated. Differential expression analysis was then performed to obtain differentially expressed genes (DEGs) in key cells with adj.P<0.05.

WGCNA

Based on ICDRGs, single-sample gene set enrichment analysis (ssGSEA) scores of each GC samples were computed. GC samples were split into high and low scoring cohorts using surv_cutpoint function. WGCNA was conducted to determine modules that showed a significant correlation with ssGSEA scores, followed by the identification of genes associated with modules. The outlier samples in the TCGA-GC dataset were eliminated to ensure analysis accuracy. A soft threshold value was selected with R² close 0.9 and mean connectivity close to 0 to assure that gene interactions were maximally consistent with scale-free distributions. After that, the least gene number was set as 250 in each module. Spearman analysis was conducted to assess correlation between ssGSEA scoring groups and obtained modules. The module showing the highest correlation was termed as the key module, and genes within this module were considered as key module genes.

Identification and analysis of the differentially expressed-ICDRGs (DE-ICDRGs)

Intersection of DEGs and key module genes yielded DE-ICDRGs. To investigate molecular functions and effect mechanisms, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using “clusterProfiler” package (version 4.6.2) (30) with adj.P<0.05. The default statistical test method used for the enrichment analysis was the Hypergeometric test, and the BH method was employed to correct for multiple testing of P values.

Selection of the prognostic genes for GC

Univariate Cox regression analysis was conducted on the DE-ICDRGs to identify GC survival-related genes with hazard ratio (HR) ≠1 and P<0.05. Additionally, least absolute shrinkage and selection operator (LASSO) regression analysis was carried out using the “glmnet” package (version 4.1-8) (31) to pinpoint prognostic genes. Expression of prognostic genes was then compared with the Wilcoxon test (P<0.05).

Construction and validation of the prognostic model

Relied above obtained prognostic genes in TCGA-GC dataset, a prognostic model was constructed. Risk score of GC samples was calculated by following formula:

Riskscore=∑i=0nβi×xi[1]

where βi represented LASSO coefficient of each prognostic gene, and xi denoted expression level of each gene. GC samples were categorized into high and low-risk groups based on the median value of risk score. Survival curves were then drafted by ranking risk score. Efficacy of prognostic model was evaluated by a receiver operator characteristic (ROC) curve, drafted by “survivalROC” package (version 0.4) (32). The survival difference of GC samples in high- and low-risk cohorts was displayed by Kaplan-Meier (KM) curve, compared by log-rank test (P<0.05). Besides, performance of prognostic model was evaluated in the validation set.

Establishment of the nomogram

To enhance prediction accuracy, risk score combined with clinical features were included in univariate (HR ≠1, P<0.05) and multivariate Cox analyses (HR ≠1, P<0.20) to screen independent prognostic factors. Subsequently, nomogram was created to predict GC patients’ survival probability. Each factor had a score, and sum of score was total points. The higher total point represented to the lower survival probability. Furthermore, accuracy of nomogram was assessed by calibration and ROC curves.

Immune microenvironment analysis

The “CIBERSORT” package (version 0.1.0) (33) was utilized to research proportions of 22 immune cells in different groups, visualized by “ggplot2” package (version 0.1.0) (34). Then, the difference of immune cells was compared by Wilcoxon test in distinct cohorts (P<0.05). In addition, correlation of immune cells and prognostic genes was also explored. Tumour immune dysfunction and exclusion (TIDE) scores of GC samples were calculated in TIDE website, compared between two cohorts via Wilcoxon test (P<0.05). Finally, the ImmuneScore, StromalScore, and EstimateScore were calculated by “estimate” package (version 1.0.13) (35), and compared by Wilcoxon test (P<0.05).

Gene set enrichment analysis (GSEA) and drugs sensitivity analysis

Between two risk cohorts, differential expression analysis was conducted, and log₂fold change (FC) value was the ranking thresholds. GSEA was then performed through “clusterProfiler” package with background gene set “c2.cp.kegg.symbols.gmt” and “h.all.v2024.1.Hs.symbols.gmt”, downloading in GSEA website. The screening criteria of GSEA were |NES| >1 and adj.P<0.05. To evaluate sensitivity for chemotherapeutics, half maximal inhibitory concentration (IC₅₀) value of drugs for distinct risk cohorts was estimated using “oncoPredict” package (P<0.05).

Statistical analysis

All analyses were executed in R software (v 4.2.2). Differences between groups were analyzed by Wilcoxon test. P<0.05 was considered statistically significant.

Results

A total of 10 types of cell subpopulations were annotated

To identify single cell signatures of GC samples, scRNA-seq analysis was conducted on the GSE183904 dataset, 151,070 cells and 25,370 genes retained for further analysis (Figure S1A). Furthermore, top 2,000 high variable genes were selected, with the top 10 genes shown in Figure S1B. PCA was performed, and turning point graph flattened at 30 dimensions, so top 30 PCs were chosen in analysis (Figure 1A). Twelve cell clusters were then identified (Figure 1B), followed by annotation of 10 distinct cell subpopulations, such as CD8⁺ natural killer T-like (NKT-like), cancer, memory B, myeloid, endothelial, stromal, lymphoid, mast, memory CD8⁺ T, and neuroendocrine cells (Figure 1C,1D). To diminish doublets of single cell data, “DoubletFinder” package was employed, removing 11,330 cells (7.5%), leaving 139,740 cells for analysis (Figure S1C). Doublet of endothelial cells and stromal cells had the most proportion than other cell subpopulations (Figure 1E). By removing doublets, it was possible to effectively reduce the interference of technical noise.

Figure 1 scRNA-seq analysis in the GSE183904 dataset. (A) Top 30 PCs. (B) The identification of cell clusters. (C) The marker genes of cell subpopulations. (D) The clustering of cell subpopulations. (E) The proportion of doublet in cell subpopulations. NKT-like, natural killer T-like; PCs, principle components; scRNA-seq, single cell RNA sequencing; UMAP, uniform manifold approximation and projection.

CD8⁺ NKT-like cells were identified as key cells in GC samples

Cell subpopulations expression in distinct samples was analyzed, revealing that CD8⁺ NKT-like cells were significantly higher in GC samples (Figure 2A). Additionally, ICDRGs-related score of CD8⁺ NKT-like cells was significantly higher in GC samples compared to controls (Figure 2B,2C), marking these cells as key cells for pseudotime trajectories analysis. Through pseudotime trajectories analysis, entire stage of CD8⁺ NKT-like cells differentiation were appeared in GC samples (Figure 2D). The expression of ICDRGs during differentiation revealed that HSP90AA1 and HMGB1 were highly expressed across all differentiation stage, while FOXP3, CD4, and IL1R1 had more expression in the early stage, and IFNG, CD8A, CD8B, and PRF1 were more expressed in late stage (Figure S1D). Finally, communication of cell subpopulations indicated that cancer cells, stromal cells, and endothelial cells had more and stronger interactions with other cells in GC samples (Figure 2E). Differential expression analysis of CD8⁺ NKT-like cells identified 788 DEGs between GC and control samples (see table online https://cdn.amegroups.cn/static/public/jgo-2025-aw-831-1.xlsx).

Figure 2 CD8⁺ NKT-like cells were identified as key cells in GC. (A) The proportion of cell subpopulations in GC samples and control samples in the GSE183904 dataset. (B) The difference of cell subpopulations between GC and control samples. (C) The AUCell algorithm of cell subpopulations. (D) Pseudotime trajectories analysis of CD8⁺ NKT-like cells. (E) The cell communication of cell subpopulations in GC and control samples. *, P<0.05; ***, P<0.001; ****, P<0.0001. GC, gastric cancer; NKT-like, natural killer T-like; UMAP, uniform manifold approximation and projection.

Identification and analysis of 143 DE-ICDRGs correlated with GC

Based on ICDRGs-related ssGSEA scores, GC samples in the TCGA-GC dataset were categorized into high and low scoring groups (minprop =0.540), with most ICDRGs showing higher expression in high scoring group (Figure 3A). To find ICD-related module genes, WGCNA was performed, and 4 outlier samples were removed, remaining 379 samples for analysis (Figure 3B). A soft threshold value of 9 was selected for WGCNA to ensure scale-free network characteristics (Figure 3C). By setting the least number of gene as 250, a total of 6 modules were obtained (Figure 3D). Among these modules, green module had the highest correlation with ssGSEA score, was selected as the key module (|cor| =0.44, P<0.05), with identification of 1,204 key genes (Figure 3E). The intersection of 788 DEGs and 1,204 key genes yielded 143 DE-ICDRGs (Figure 3F). Enrichment analyses of these DE-ICDRGs identified 644 GO terms, including antigen processing and presentation, endocytic vesicle membrane, and antigen binding, and 42 KEGG pathways, containing antigen processing and presentation, EBV infection, and cell adhesion molecules (Figure 3G,3H).

Figure 3 Identification and analysis of 143 DE-ICDRGs in the TCGA-GC dataset. (A) The distribution of ICDRGs in high and low ssGSEA scoring groups. (B) Clustering analysis of WGCNA. (C) Selection of soft threshold value. (D) Dynamic tree cutting analysis. (E) Correlation of modules with ssGSEA score. (F) Venn diagram of module genes and key genes in CD8⁺ NKT-like cells. (G,H) GO (G) and KEGG (H) enrichment analyses of DE-ICDRGs. BP, biological process; CC, cellular component; DE-ICDRGs, differentially expressed-immunogenic cell death-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; MHC, major histocompatibility complex; NKT-like, natural killer T-like; ssGSEA, single-sample gene set enrichment analysis; TAP, transporter associated with antigen processing; TCGA, The Cancer Genome Atlas; WGCNA, weighted gene co-expression network analysis.

Identification and expression analysis of prognostic genes in GC

Univariate Cox regression identified 7 genes—CD69, CXCR4, GLUL, GLIPR1, RAB8B, TAP1, and PDE4B—as GC survival-related genes (Figure 4A). Further selection of prognostic genes for GC samples at lambda.min =0.0135 revealed CXCR4, GLUL, GLIPR1, RAB8B, and TAP1 (Figure 4B). Through expression analysis, we found that all prognostic genes had significant difference between distinct samples (P<0.05) (Figure 4C). Except CLUL, remaining prognostic genes had higher expression in GC samples.

Figure 4 Identification and analysis of prognostic genes. (A) Selection of 7 genes by univariate Cox regression analysis. (B) LASSO analysis for prognostic genes. (C) The expression of prognostic genes in GC and control samples. *, P<0.05; **, P<0.01; ****, P<0.0001. CI, confidence interval; GC, gastric cancer; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator.

Prognostic model predicted the survival conditions of GC samples

Prognostic model was constructed based on selected genes, and risk scores for GC patients were determined. Using risk score median, GC samples from both TCGA-GC and GSE62254 datasets were classified into high- and low-risk cohorts. As risk score increased, GC samples death probability also rose in both datasets (Figure 5A,5B). Furthermore, GC survival probability in different cohorts was evaluated by KM survival curves, and results revealed that low-risk patients had higher survival probability in both datasets (P<0.05) (Figure 5C,5D). To investigate the efficacy of prognostic model, ROC curves were drafted in both datasets, and area under the curve (AUC) values in diverse years were all greater than 0.6 (Figure 5E,5F), which indicated that prognostic model predicted the survival probability of GC samples.

Figure 5 Construction and validation of prognostic model in the TCGA-GC and GSE62254 datasets. (A,B) Survival curves of GC samples in training (A) and validation (B) sets. (C,D) Kaplan-Meier survival curves of GC samples in training (C) and validation (D) sets. (E,F) ROC curves for prognostic model in training (E) and validation (F) sets. AUC, area under the curve; GC, gastric cancer; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operator characteristic; TCGA, The Cancer Genome Atlas.

The nomogram could predict survival of GC patients

To develop a prognostic model, univariate and multivariate Cox analyses selected risk score, age, and N/M stage as independent prognostic factors in GC samples (Figure 6A,6B). Relied on above factors, a nomogram was created to predict GC’s survival probability (Figure 6C). The calibration curve closed to ideal curve (Figure 6D), and AUC value of ROC curve was greater than 0.69 (Figure 6E), which indicated that nomogram had an ability to predict the survival probability for GC patients.

Figure 6 Construction and validation of the nomogram. (A,B) Univariate (A) and multivariate (B) Cox regression analysis for the selection of independent prognostic factors. (C) The nomogram constructed by independent prognostic factors. (D) Calibration curve for nomogram. (E) ROC curves of nomogram. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; M, metastasis; N, node; OS, overall survival; ROC, receiver operator characteristic; T, tumor.

Immune conditions in GC samples

Disruption of the immune microenvironment is a major cause of poor prognosis in GC (36), so this study explored immune differences between different groups, with 13 types, such as naive B cells and naive CD4 T cells, showing significant difference (Figure 7A,7B). Further correlation analysis indicated that TAP1 had the highest positive and negative correlation with activated memory CD4 T cells (cor =0.47, P<0.05) and resting memory CD4 T cells (cor =−0.34, P<0.05), respectively (Figure 7C). Moreover, TIDE, ImmuneScore, StromalScore, and EstimateScore were higher in high-risk patients (Figure 7D,7E), suggesting that these patients had greater immune escape and poorer immunotherapy outcomes.

Figure 7 Immune infiltration analysis in the TCGA-GC dataset. (A) The immune landscape of immune cells in GC and control samples. (B) The difference of immune cells between GC and control samples. (C) The correlation of immune cells with prognostic genes. (D) The difference of TIDE score between high and low risk cohorts. (E) The difference of StromalScore, ImmuneScore and EstimateScore between high and low risk cohorts. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. GC, gastric cancer; NK, natural killer; TCGA, The Cancer Genome Atlas; TIDE, tumour immune dysfunction and exclusion.

Enriched pathways and drug sensitivity analyses in GC samples

To explore enriched pathways in GC samples, GSEA was performed, revealing 53 KEGG pathways and 4 hallmark pathways. Among top five pathways, DNA replication was increased, while origin unwinding and elongation, pre-IC formation, variant Igh mmset fusion to transcriptional activation, and neuroactive ligand receptor interaction were decreased (Figure 8A). Furthermore, enriched hallmark pathways included kras signaling dn, myogenesis, and pancreas beta cells (Figure 8B), all showing increased trends. Additionally, drugs sensitivity analysis identified 198 drugs with significant different IC₅₀ value between two cohorts. Top 10 with the most markedly different drugs were displayed in Figure 8C. AZD8186, BMS.754807, Doramapimod, JQ1, NU7441 had higher IC₅₀ value in low-risk cohort, while Erlotinib, Gefitinib, Lapatinib, SCH772984, and VX.11e had lower IC₅₀ value in low-risk cohort, indicating that drug choices could be optimized based on patient’s risk group for better treatment efficacy.

Figure 8 Enriched pathways and drug sensitivity analysis in the TCGA-GC dataset. (A,B) Enriched KEGG (A) and hallmark pathways (B) in GC samples. (C) The difference of half maximal inhibitory concentration of drugs between high- and low-risk cohorts. GC, gastric cancer; IC₅₀, half maximal inhibitory concentration; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA, The Cancer Genome Atlas.

Discussion

GC remains a formidable clinical challenge due to its aggressive progression and poor prognosis, underscoring the urgent need for robust prognostic genes and therapeutic targets (37). Programmed cell death (PCD) refers to a series of active and orderly cell death processes regulated by genes, including apoptosis, programmed necrosis, pyroptosis, ferroptosis, and other types (38,39). ICD, as a special form of PCD (40), is a regulated cell death pattern that mediates the interaction between innate and adaptive immunity and is now recognized as a key factor influencing tumor immunogenicity and treatment resistance (41). Our scRNA-seq analysis delineated 10 cellular subpopulations, identifying CD8⁺ NKT-like cells as the key cells in the GC TME. These cells exhibited significantly elevated abundance and ICDRG scores in GC samples, compared with controls. Pseudotime analysis revealed their dynamic functional states during GC progression, characterized by distinct ICDRG expression patterns. The functional impairment of CD8⁺ NKT-like cells within the GC TME is consistent with broader immunosuppressive mechanisms seen in peritoneal metastases. In these metastases, cytotoxic lymphocytes, including CD8⁺ T and NKT-like cells, are significantly depleted, with their levels inversely correlating with tumor burden (42). This lymphocyte attrition creates a permissive niche for tumor progression, further exacerbated by M2-like macrophage expansion. Notably, the loss of cytotoxic function in senescent NKT-like cells—characterized by glucocorticoid receptor downregulation and cytokine dysregulation (43)—may parallel the functional exhaustion observed in our GC cohort, suggesting conserved mechanisms of immune evasion across malignancies. In summary, by combining the previous single-cell research conclusions (19,44,45) with the scRNA-seq analysis results of this study, we believe that CD8⁺ NKT-like cells represent a dynamic state that is regulated by the environment rather than a fixed cell subset with a specific function. The expression level of these cells is associated with the expression of ICDRGs, but the causal relationship between the two cannot be clearly inferred from the current data.

Integrating multi-omics data, we systematically identified five ICD-related prognostic genes: CXCR4, GLUL, GLIPR1, RAB8B, and TAP1 in GC samples. CXCR4, a chemokine receptor implicated in tumor metastasis and immune evasion, was upregulated in GC tissues compared to adjacent non-tumor tissues, aligning with prior studies linking its overexpression to advanced disease stages and reduced survival (46,47). Mechanistically, CXCR4 facilitates GC progression by enhancing glycolytic metabolism and promoting cell migration and invasion (48). Cao et al. demonstrated that CXCR4 mediates Treg-cell-dependent immunosuppression via regulating the Rho-GTPase/NF-κB signaling axis, thereby facilitating tumor immune evasion (49). Meanwhile, our results suggested its role as a driver of immunosuppression and metastatic adaptation, suggesting that targeting CXCR4 could disrupt both tumor-intrinsic and microenvironmental pathways. GLUL, a key enzyme in glutamine metabolism, has been associated with poor prognosis in multiple cancers, including hepatocellular carcinoma and breast cancer (50,51). Notably, Mendelian randomization analysis (52) identified GLUL as a putative risk factor for GC, which aligns with the upregulation of GLUL in high-risk GC patients in our study, suggesting causality and enhancing the biological significance of our finding. GLUL-driven metabolic reprogramming may sustain tumor proliferation under nutrient stress (53), highlighting its dual role as a metabolic checkpoint and prognostic indicator. In addition, relevant studies employing cellular assays have elucidated a novel mechanism through which GLUL contributes to metal-induced carcinogenesis and cancer progression; this mechanism is governed by the MYC/METTL16/YTHDC1/GLUL signaling cascade, highlighting the tumor-promoting function of GLUL in cancer progression (54). GLIPR1, a glioma pathogenesis-related protein, exhibited significant prognostic value in our cohort. Although GLIPR1 has been less well studied in GC, Ye et al. still found that knockdown of this gene enhances CAR-T cell cytotoxicity against GC by alleviating immunosuppressive signals (55). This is consistent with our finding that GLIPR1 overexpression correlates with immune exclusion and poor survival, positioning it as a novel therapeutic target for reversing immune tolerance. RAB8B, a small GTPase involved in vesicular trafficking, has not been extensively studied in GC. However, its oncogenic roles in breast and pancreatic cancers (56,57) suggest that conserved mechanisms may also play a role in GC progression. RAB8B may regulate cytokine secretion or immune cell recruitment, warranting functional validation to elucidate its contribution to ICD-related pathways. TAP1, a critical mediator of antigen processing and presentation, emerged as a key determinant of immune responsiveness. TAP1 deficiency compromises tumor immunogenicity and predicts resistance to ICIs (58). Our data support its synergistic role with oxaliplatin and anti-programmed death-ligand 1 (PD-L1) therapies (59), providing a molecular rationale for combining ICD inducers with immunotherapy to enhance clinical outcomes. In summary, these prognostic genes are associated with various pathways, like chemokine signaling, metabolic reprogramming, immune checkpoint regulation, intracellular trafficking, and antigen presentation. Their collective dysregulations induce a high-risk GC phenotype characterized by aggressive biology and immune suppression.

The risk model constructed using these five genes effectively stratified GC patients into distinct risk groups in both training and validation sets, with high-risk patients exhibiting significantly shorter OS (P<0.001) and increased mortality. To enhance clinical applicability, we integrated the risk score with independent clinical prognostic factors (age, N stage, M stage) into a comprehensive nomogram. This nomogram provided a user-friendly tool for clinicians to estimate individual patient survival probabilities at 3, 4, and 5 years. Similar integrated models have proven valuable in guiding treatment decisions in other cancers (60,61). Our nomogram facilitates personalized risk assessment, potentially aiding in selecting patients for more intensive surveillance or novel therapeutic approaches. A recent study has constructed interpretable pre-trained multi-omics models, which integrate multi-level molecular data and have been applied to multi-task learning in the field of oncology, demonstrating superior robustness, interpretability, and cross-task generalization capabilities (62). Moreover, given that the genetic and transcriptomic structures may exhibit significant heterogeneity among different populations, extending this analysis framework to large-scale national biological sample banks and population genomics studies for verification is of crucial significance for evaluating the universality and clinical application value of the selected ICD-related features in different patient populations (63). Based on the current research progress, we will subsequently conduct studies on the optimization of multi-omics models and the clinical validation of ICD features, providing new ideas and evidence for tumor immunotherapy.

Disruption of the tumor immune microenvironment (TIME) is a well-established driver of poor prognosis in GC, contributing to tumor progression, metastasis, and therapy resistance (36). To characterize the immune landscape associated with our prognostic genes-based risk stratification, we performed comprehensive immune infiltration analysis. Our results revealed significant differences in the abundance of 13 distinct immune cell types between high- and low-risk GC groups. Notably, among the prognostic genes identified, TAP1 exhibited the strongest positive correlation with activated memory CD4⁺ T cells (cor =0.47, P<0.001) and a significant inverse correlation with resting memory CD4⁺ T cells (cor =−0.34, P<0.001). This finding aligns with the critical role of TAP1 in antigen processing and presentation. TAP1 deficiency impairs MHC-II-mediated antigen presentation (58), leading to suboptimal T cell receptor (TCR) engagement. This suboptimal signaling favors the differentiation of CD4⁺ T cells towards an anergic or quiescent memory phenotype (resting) rather than an activated effector state (activated) (64,65). Consequently, the accumulation of resting memory CD4⁺ T cells, which are functionally inert and often express inhibitory receptors like PD-1 and CTLA-4 (66), fails to adequately support DC licensing or CD8⁺ T cell cytotoxicity. This creates an immunosuppressive niche conducive to tumor progression. Collectively, our immune infiltration analysis not only quantifies the differential abundance of specific immune cell populations (13 types) associated with the ICD-based risk signature but also mechanistically links TAP1 dysregulation to a specific functional impairment (reduced antigen presentation) and its immunological consequence (skewed CD4⁺ T cell differentiation), providing a deeper understanding of the immunosuppressive mechanisms operational in high-risk GC. It should be noted that the immune cell infiltration analysis in this study is based on the CIBERSORT algorithm for bulk RNA sequencing, and inherent variability in results exists between different algorithms. Therefore, the aforementioned immune analysis results should be regarded solely as supportive evidence and a tool for hypothesis generation. The reliability of these conclusions requires further validation through experimental methods such as single-cell sequencing and flow cytometry.

In this study, we identified five prognostic genes—CXCR4, GLUL, GLIPR1, RAB8B, TAP1—associated with ICD in GC samples by comprehensive analysis, which provided new targets for the treatment of GC. However, it has to be admitted that our study has some limitations. Firstly, our study focuses on the key cell type (CD8⁺ NKT-like cells) and prognostic gene bioinformatics aspect, but lacks the validation of in vitro and in vivo experiments, so further experimental validation is necessary. Secondly, multigene expression characteristics can serve as an integrated indicator for characterizing complex pathological states, and their cross-disciplinary application value is worthy of in-depth exploration. Additionally, the prognostic genes identified in this study via Cox regression analysis only reflect a statistical correlation with clinical outcomes, rather than clarifying whether this association represents a direct causal effect, an indirect regulatory role, or a concomitant correlation mediated by unmeasured confounding factors. Thus, the causal relationship between this gene set and tumor prognosis awaits further verification in subsequent studies. Whole-genome duplication (WGD) is a widespread event of complete chromosome set replication in various tumors, playing a crucial role in regulating gene expression and driving cell carcinogenesis (67). Given this, it is worthy of in-depth exploration whether the abnormal expression of the prognostic-related genes identified in this study is associated with WGD. Finally, the scRNA-seq data, while providing valuable insights into cellular heterogeneity, represents a static snapshot; longitudinal studies tracking ICD dynamics and immune cell interactions during disease progression or treatment are needed.

Conclusions

This study identified CD8^⁺ NKT-like cells as key cells and five ICD-related prognostic genes—CXCR4, GLUL, GLIPR1, RAB8B, TAP1 in GC samples through integrated scRNA-seq and multi-omics bioinformatics approaches. Based on prognostic genes, a prognostic model was constructed and it predicted well the survival probability of GC samples. Furthermore, we elucidated the immunosuppressive microenvironment associated with high-risk patients—particularly TAP1-mediated CD4⁺ T-cell dysfunction—and identified differential drug sensitivity patterns between risk groups. Collectively, these findings provided a framework for personalized risk assessment and precision immunotherapy in GC samples.

Acknowledgments

We acknowledge TCGA and GEO database for providing their platforms and contributors for uploading their meaningful datasets.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-aw-831/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-aw-831/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-aw-831/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Smyth EC, Nilsson M, Grabsch HI, et al. Gastric cancer. Lancet 2020;396:635-48. [Crossref] [PubMed]
2. Ajani JA, D'Amico TA, Bentrem DJ, et al. Gastric Cancer, Version 2.2025, NCCN Clinical Practice Guidelines In Oncology. J Natl Compr Canc Netw 2025;23:169-91. [Crossref] [PubMed]
3. Huang X, Liu Y, Qian C, et al. CHSY3 promotes proliferation and migration in gastric cancer and is associated with immune infiltration. J Transl Med 2023;21:474. [Crossref] [PubMed]
4. Machlowska J, Baj J, Sitarz M, et al. Gastric Cancer: Epidemiology, Risk Factors, Classification, Genomic Characteristics and Treatment Strategies. Int J Mol Sci 2020;21:4012. [Crossref] [PubMed]
5. Zhao W, Lin J, Cheng S, et al. Comprehensive analysis of COMMD10 as a novel prognostic biomarker for gastric cancer. PeerJ 2023;11:e14645. [Crossref] [PubMed]
6. Bejarano L, Jordāo MJC, Joyce JA. Therapeutic Targeting of the Tumor Microenvironment. Cancer Discov 2021;11:933-59. [Crossref] [PubMed]
7. Hinshaw DC, Shevde LA. The Tumor Microenvironment Innately Modulates Cancer Progression. Cancer Res 2019;79:4557-66. [Crossref] [PubMed]
8. Bader JE, Voss K, Rathmell JC. Targeting Metabolism to Improve the Tumor Microenvironment for Cancer Immunotherapy. Mol Cell 2020;78:1019-33. [Crossref] [PubMed]
9. Bai R, Cui J. Development of Immunotherapy Strategies Targeting Tumor Microenvironment Is Fiercely Ongoing. Front Immunol 2022;13:890166. [Crossref] [PubMed]
10. Wei G, Zhang H, Zhao H, et al. Emerging immune checkpoints in the tumor microenvironment: Implications for cancer immunotherapy. Cancer Lett 2021;511:68-76. [Crossref] [PubMed]
11. Jin MZ, Wang XP. Immunogenic Cell Death-Based Cancer Vaccines. Front Immunol 2021;12:697964. [Crossref] [PubMed]
12. Chiaravalli M, Spring A, Agostini A, et al. Immunogenic Cell Death: An Emerging Target in Gastrointestinal Cancers. Cells 2022;11:3033. [Crossref] [PubMed]
13. Yu S, Xiao H, Ma L, et al. Reinforcing the immunogenic cell death to enhance cancer immunotherapy efficacy. Biochim Biophys Acta Rev Cancer 2023;1878:188946. [Crossref] [PubMed]
14. Arimoto KI, Miyauchi S, Liu M, et al. Emerging role of immunogenic cell death in cancer immunotherapy. Front Immunol 2024;15:1390263. [Crossref] [PubMed]
15. Zhao S, Sun D, Yu H, et al. Oxaliplatin accelerates immunogenic cell death by activating the cGAS/STING/TBK1/IRF5 pathway in gastric cancer. FEBS J 2025;292:3814-28. [Crossref] [PubMed]
16. Zhu X, Zheng W, Wang X, et al. Enhanced Photodynamic Therapy Synergizing with Inhibition of Tumor Neutrophil Ferroptosis Boosts Anti-PD-1 Therapy of Gastric Cancer. Adv Sci (Weinh) 2024;11:e2307870. [Crossref] [PubMed]
17. Wang L, Chen J, Zuo Q, et al. Calreticulin enhances gastric cancer metastasis by dimethylating H3K9 in the E-cadherin promoter region mediating by G9a. Oncogenesis 2022;11:29. [Crossref] [PubMed]
18. Kong Q, Zhang Z, Liang Z. Upregulating miR-637 aggravates endoplasmic reticulum stress-induced apoptosis in gastric cancer cells by suppressing Calreticulin. Anim Cells Syst (Seoul) 2020;24:267-74. [Crossref] [PubMed]
19. Kim S, Leem G, Choi J, et al. Integrative analysis of spatial and single-cell transcriptome data from human pancreatic cancer reveals an intermediate cancer cell population associated with poor prognosis. Genome Med 2024;16:20. [Crossref] [PubMed]
20. Ji G, Yang Q, Wang S, et al. Single-cell profiling of response to neoadjuvant chemo-immunotherapy in surgically resectable esophageal squamous cell carcinoma. Genome Med 2024;16:49. [Crossref] [PubMed]
21. Fang K, Ohihoin AG, Liu T, et al. Integrated single-cell analysis reveals distinct epigenetic-regulated cancer cell states and a heterogeneity-guided core signature in tamoxifen-resistant breast cancer. Genome Med 2024;16:134. [Crossref] [PubMed]
22. Luo Q, Teschendorff AE. Cell-type-specific subtyping of epigenomes improves prognostic stratification of cancer. Genome Med 2025;17:34. [Crossref] [PubMed]
23. Li D, Zhou Z, Hui Y, et al. EHBMT, a method for visualizing tumor evolution, identifies a surge in gastric cancer with hybrid epithelial-mesenchymal phenotypes due to an inflammatory microenvironment. Exp Mol Med 2025;57:2440-57. [Crossref] [PubMed]
24. Dai SL, Pan JQ, Su ZR. Multi-omics features of immunogenic cell death in gastric cancer identified by combining single-cell sequencing analysis and machine learning. Sci Rep 2024;14:21751. [Crossref] [PubMed]
25. Hao Y, Stuart T, Kowalski MH, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42:293-304. [Crossref] [PubMed]
26. Kumar V, Ramnarayanan K, Sundar R, et al. Single-Cell Atlas of Lineage States, Tumor Microenvironment, and Subtype-Specific Expression Programs in Gastric Cancer. Cancer Discov 2022;12:670-91. [Crossref] [PubMed]
27. Zhu S, Cheng Q, Zou M, et al. Combining bulk and scRNA-seq to explore the molecular mechanisms governing the distinct efferocytosis activities of a macrophage subpopulation in PDAC. J Cell Mol Med 2024;28:e18266. [Crossref] [PubMed]
28. Cao L, Lu X, Wang X, et al. From single-cell to spatial transcriptomics: decoding the glioma stem cell niche and its clinical implications. Front Immunol 2024;15:1475235. [Crossref] [PubMed]
29. Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods 2017;14:979-82. [Crossref] [PubMed]
30. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
31. Simon N, Friedman J, Hastie T, et al. Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J Stat Softw 2011;39:1-13. [Crossref] [PubMed]
32. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med 2013;32:5381-97. [Crossref] [PubMed]
33. Zhang WY, Chen ZH, An XX, et al. Analysis and validation of diagnostic biomarkers and immune cell infiltration characteristics in pediatric sepsis by integrating bioinformatics and machine learning. World J Pediatr 2023;19:1094-103. [Crossref] [PubMed]
34. Cao T, Li Q, Huang Y, et al. plotnineSeqSuite: a Python package for visualizing sequence data using ggplot2 style. BMC Genomics 2023;24:585. [Crossref] [PubMed]
35. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
36. Hong Y, Luo C, Liu Y, et al. Development of a PANoptosis-related LncRNAs for prognosis predicting and immune infiltration characterization of gastric Cancer. Sci Rep 2025;15:7699. [Crossref] [PubMed]
37. Joshi SS, Badgwell BD. Current treatment and recent progress in gastric cancer. CA Cancer J Clin 2021;71:264-79. [Crossref] [PubMed]
38. van Doornik MC, Cats BP, Barth PG, et al. Intra-uterine tachycardia associated with multicystic encephalomalacia (MCE). Eur J Obstet Gynecol Reprod Biol 1985;20:191-5. [Crossref] [PubMed]
39. Yang J, Hwang CS. Nα-terminal acetylation meets ferroptosis via N-degron pathway. Mol Cells 2024;47:100160. [Crossref] [PubMed]
40. Chen SP, Tian J. Bibliometric analysis of programmed cell death and immunogenic cell death in hepatocellular carcinoma immunotherapy: global trends and future directions. Discov Oncol 2025;16:1208. [Crossref] [PubMed]
41. Galluzzi L, Vitale I, Warren S, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J Immunother Cancer 2020;8:e000337. [Crossref] [PubMed]
42. Takahashi K, Kurashina K, Yamaguchi H, et al. Altered intraperitoneal immune microenvironment in patients with peritoneal metastases from gastric cancer. Front Immunol 2022;13:969468. [Crossref] [PubMed]
43. Hodge G, Jersmann H, Tran HB, et al. Lymphocyte senescence in COPD is associated with loss of glucocorticoid receptor expression by pro-inflammatory/cytotoxic lymphocytes. Respir Res 2015;16:2. [Crossref] [PubMed]
44. Choi Y, Shim M, Kim SH, et al. Single-cell transcriptomics of the myeloid milieu reveals an angiogenic niche in triple-negative breast cancer. Exp Mol Med 2025;57:2487-504. [Crossref] [PubMed]
45. Deng Y, Guo C, Liu X, et al. Single-cell transcriptomic profiling reveals liver fibrosis in colorectal cancer liver metastasis. Exp Mol Med 2025;57:2517-32. [Crossref] [PubMed]
46. Chen G, Zhou Z, Jin J, et al. CXCR4 is a prognostic marker that inhibits the invasion and migration of gastric cancer by regulating VEGF expression. Oncol Lett 2021;22:587. [Crossref] [PubMed]
47. Yu S, Wu T, Wang J, et al. Combined Evaluation of Expression of CXCR4 and Nrf2 as Prognostic Factor for Patients with Gastric Carcinoma. Anticancer Agents Med Chem 2018;18:388-93. [Crossref] [PubMed]
48. Xu L, Liu J, An Y, et al. Glycolysis-related genes predict prognosis and indicate immune microenvironment features in gastric cancer. BMC Cancer 2024;24:979. [Crossref] [PubMed]
49. Cao C, Xu M, Peng T, et al. Blocking CXCR4(+) CD4(+) T cells reprograms T(reg)-mediated immunosuppression via modulating the Rho-GTPase/NF-κB signaling axis. Genome Med 2025;17:85. [Crossref] [PubMed]
50. Liu P, Lu D, Al-Ameri A, et al. Glutamine synthetase promotes tumor invasion in hepatocellular carcinoma through mediating epithelial-mesenchymal transition. Hepatol Res 2020;50:246-57. [Crossref] [PubMed]
51. Wang Y, Fan S, Lu J, et al. GLUL Promotes Cell Proliferation in Breast Cancer. J Cell Biochem 2017;118:2018-25. [Crossref] [PubMed]
52. Fan W, Wu Z, Xu S, et al. Mendelian randomization analysis of plasma proteins reveals potential novel tumor markers for gastric cancer. Sci Rep 2025;15:3537. [Crossref] [PubMed]
53. Bott AJ, Maimouni S, Zong WX. The Pleiotropic Effects of Glutamine Metabolism in Cancer. Cancers (Basel) 2019;11:770. [Crossref] [PubMed]
54. Xie YX, Wang L, Zhou ZH, et al. m(6)A RNA methyltransferase METTL16 induces Cr(VI) carcinogenesis and lung cancer development through glutamine biosynthesis and GLUL expression. J Hazard Mater 2024;480:136093. [Crossref] [PubMed]
55. Ye J, Liu Q, He Y, et al. Combined therapy of CAR-IL-15/IL-15Rα-T cells and GLIPR1 knockdown in cancer cells enhanced anti-tumor effect against gastric cancer. J Transl Med 2024;22:171. [Crossref] [PubMed]
56. Quinting T, Heymann AK, Bicker A, et al. Myoglobin Protects Breast Cancer Cells Due to Its ROS and NO Scavenging Properties. Front Endocrinol (Lausanne) 2021;12:732190. [Crossref] [PubMed]
57. Anand S, Khan MA, Khushman M, et al. Comprehensive Analysis of Expression, Clinicopathological Association and Potential Prognostic Significance of RABs in Pancreatic Cancer. Int J Mol Sci 2020;21:5580. [Crossref] [PubMed]
58. Tu Z, Li K, Ji Q, et al. Pan-cancer analysis: predictive role of TAP1 in cancer prognosis and response to immunotherapy. BMC Cancer 2023;23:133. [Crossref] [PubMed]
59. Zhao Y, Liu Z, Deng K, et al. Identification of TAP1 as a T-cell related therapeutic target in gastric cancer by mediating oxalipliatin-related synergistic enhancement of immunotherapy. Int Immunopharmacol 2024;132:111998. [Crossref] [PubMed]
60. Wang S, Zhang S, Li X, et al. HIGD1B, as a novel prognostic biomarker, is involved in regulating the tumor microenvironment and immune cell infiltration; its overexpression leads to poor prognosis in gastric cancer patients. Front Immunol 2024;15:1415148. [Crossref] [PubMed]
61. Wang S, Wang R, Hu D, et al. Machine learning reveals diverse cell death patterns in lung adenocarcinoma prognosis and therapy. NPJ Precis Oncol 2024;8:49. [Crossref] [PubMed]
62. Wang FA, Zhuang Z, Gao F, et al. TMO-Net: an explainable pretrained multi-omics model for multi-task learning in oncology. Genome Biol 2024;25:149. [Crossref] [PubMed]
63. Lee H, Kim W, Kwon N, et al. Lessons from national biobank projects utilizing whole-genome sequencing for population-scale genomics. Genomics Inform 2025;23:8. [Crossref] [PubMed]
64. Logunova NN, Kriukova VV, Shelyakin PV, et al. MHC-II alleles shape the CDR3 repertoires of conventional and regulatory naïve CD4(+) T cells. Proc Natl Acad Sci U S A 2020;117:13659-69. [Crossref] [PubMed]
65. Ishina IA, Zakharova MY, Kurbatskaia IN, et al. MHC Class II Presentation in Autoimmunity. Cells 2023;12:314. [Crossref] [PubMed]
66. Axelrod ML, Cook RS, Johnson DB, et al. Biological Consequences of MHC-II Expression by Tumor Cells in Cancer. Clin Cancer Res 2019;25:2392-402. [Crossref] [PubMed]
67. Chang E, An JY. Whole-genome doubling is a double-edged sword: the heterogeneous role of whole-genome doubling in various cancer types. BMB Rep 2024;57:125-34. [Crossref] [PubMed]