- 1College of First Clinical Medicine, Shandong University of Traditional Chinese Medicine, Jinan, China
- 2China Institute of Sport and Health Science, Beijing Sport University, Beijing, China
- 3Department of Rheumatology and Immunology, Tongren Hospital, School of Medicine, Shanghai Jiao Tong University, Shanghai, China
- 4Department of Oncology, Affiliated Hospital of Shandong University of Traditional Chinese Medicine, Jinan, China
- 5Department of Gynecology, Longhua Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China
Background: Globally, gastric cancer (GC) stands as the fifth most prevalent form of malignant neoplasm and represents a significant contributor to mortality associated with oncological conditions. Despite advancements in therapeutic strategies for GC, the outcomes for patients with advanced stages of the disease continue to be unfavorable, largely due to tumor heterogeneity and the challenges posed by resistance to therapeutic agents. Metabolic reprogramming is pivotal in driving the advancement of GC, contributing to the development of resistance to pharmacological treatments and facilitating the cancer’s ability to evade immune surveillance. Developing multi-target comprehensive treatment strategies by integrating tumor microenvironment (TME) modulation holds promise for significantly improving therapeutic efficacy.
Methods: The study analyzed GC and identified key cell subtypes by integrating data derived from single-cell RNA-sequencing (scRNA-seq) alongside spatial transcriptomics information. Cell type identification was accomplished using the tool of Seurat, and the spatial distribution of cell types was revealed through the Robust Cell Type Decomposition technique. CellChat was used to analyze the interactions between key cell subtypes and other cells, and the “StLearn” package was employed to investigate spatial cell communication in depth. Additionally, the functional role of the key molecule ELK4 was validated through in vitro experiments.
Results: This research utilized scRNA-seq combined with spatial transcriptomics to comprehensively analyze GC, identifying the C1 NDUFAB1+ subtype, which exhibited high proliferative activity, metabolic reprogramming capabilities, and immune evasion properties. It was found that the C1 NDUFAB1+ subtype closely interacted with fibroblasts and pericytes via the PARs signaling pathway. Additionally, in vitro experiments confirmed that knockdown of ELK4 substantially curbed tumor cell proliferation, migration, and invasion.
Conclusion: This study revealed the main significance of the C1 NDUFAB1+ subtype in GC, elucidating its core mechanisms in tumor progression, metabolic reprogramming, and immune evasion. ELK4 was identified as a key regulatory factor that markedly enhanced the proliferation, migratory capacity, and invasive potential of tumor cells, while changes in the TME were a driving force behind immune suppression and drug resistance. The findings underscored the importance of developing specific therapeutic targets, targeting metabolic reprogramming, and overcoming immune evasion, providing new theoretical foundations.
Introduction
Gastric cancer (GC) is positioned as the fifth leading frequent type of cancer globally and holds a significant position among cancer-related causes of death, representing a significant risk to public health (1, 2). Referring to the global cancer data from 2022, the incidence of GC surpassed 968,000, with nearly 660,000 deaths, representing roughly 6.8% of total cancer mortality (3). The incidence of GC exhibits notable regional variations, with East Asia being a high-incidence region, while the African continent reports relatively lower rates (1, 3). Based on histological characteristics, clinical behavior, and molecular features, GC can be classified into several distinct types, with intestinal type and diffuse type being the two primary histological categories (4, 5). The intestinal type is frequently linked to intestinal metaplasia and Helicobacter pylori infection, showing higher prevalence in regions with elevated GC rates. On the other hand, the diffuse type is associated with worse clinical results and occurs frequently found in younger individuals and female patients (6, 7).
GC treatment primarily includes surgical resection, radiotherapy, chemotherapy, immunotherapy, and targeted treatment (8). Surgical resection remains the primary treatment for early-stage GC, but is associated with high postoperative recurrence rates and significant risk of distant metastasis (9). For advanced-stage patients, while chemotherapy can prolong survival (10), it often leads to severe adverse effects including myelosuppression and neurotoxicity (11). In recent years, immune checkpoint blockers (like PD-1/PD-L1 blockers) have exhibited certain efficacy in some GC patients (8, 12), but the response rate remains low, and some patients develop resistance after treatment. Clinical data showed that the objective response rate of single-agent PD-1/PD-L1 inhibitors in advanced GC patients was approximately 9%-14% (13, 14). More importantly, some advanced GC patients developed resistance following immune checkpoint inhibitor treatment (15). Drug resistance may arise from various molecular mechanisms, such as abnormal activation of the epidermal growth factor receptor (EGFR) signaling cascade (16–18), dysregulation of the PI3K/AKT/mTOR pathway (19–21), and the role of immunosuppressive cells in the tumor microenvironment (TME) (22). Additionally, research has demonstrated that alterations in cellular metabolism greatly aid in tumor advancement and resistance to therapeutic agents, and immune evasion of GC (23). GC cells undergo metabolic reprogramming, such as enhanced glycolysis, glutamine metabolism, and fatty acid synthesis, to fulfill the necessary energy for rapid cell division and maintain survival advantages (24, 25). Therefore, combining research on metabolic reprogramming with TME modulation strategies to develop multi-target, multi-level comprehensive treatment approaches hold promise for significantly boosting the efficacy of GC therapy and providing new breakthroughs to overcome drug resistance and tumor heterogeneity.
The recent innovations of single-cell RNA-sequencing (scRNA-seq) platform has offered superior precision in the study of GC (26–29). This technology helps identify cell subtypes and metabolic features associated with drug resistance (30–32), thereby guiding the combination of targeted drugs and immune checkpoint inhibitors. It can also be used to monitor dynamic changes in tumor cells and the TME during treatment, providing real-time evidence for adjusting treatment plans. Furthermore, the introduction of spatial transcriptomics (ST) has further expanded the depth of research, with its advantage lying in the ability to resolve the spatial distribution and interactions of different cell types within the TME, offering a more comprehensive perspective on tumor heterogeneity and TME complexity.
In this study, we incorporated ST following the implementation of scRNA-seq to thoroughly investigate tumor heterogeneity and spatial distribution characteristics within the GC TME. Through scRNA-seq, we identified four tumor cell subtypes, among which the C1 NDUFAB1+ subtype exhibited significant metabolic activity and stem cell-like properties. Its metabolic features were closely associated with metabolic reprogramming, suggesting that this subtype may play a key driving role in tumor progression. Additionally, we discovered strong interactions between C1 NDUFAB1+ subtype and fibroblasts as well as pericytes throughout the TME, primarily mediated through the PARs signaling pathway, particularly via the PRSS3-F2R ligand-receptor pair. ST further confirmed the spatial enrichment of C1 NDUFAB1+ subtype in specific regions of GC tissues, underlining their major contribution to tumor development. To elucidate further the molecular regulatory mechanisms of C1 NDUFAB1+ subtype, we identified ELK4, a key transcription factor (TF) highly expressed in this subtype. ELK4 has been recognized as a promising target for treating multiple malignancies (33, 34), as inhibiting its expression or activity can constrain the proliferation, migration, and immune evasion of tumor cells, while enhancing sensitivity to chemotherapy and immunotherapy. To validate the regulatory capacity of ELK4 in GC, we executed in vitro tests with two GC cell lines (NCI-N87 and AGS). The findings indicated that silencing ELK4 significantly impaired the proliferation, migration, and invasion capabilities of GC cells. Furthermore, based on genes associated with C1 NDUFAB1+ subtype, we developed and verified a predictive risk scoring model, revealing that individuals in the high-risk group exhibited poorer clinical outcomes, and their TME exhibited significant immunosuppressive characteristics, further supporting the significant impact of C1 NDUFAB1+ subtype in tumor progression and immune evasion.
Our study comprehensively revealed the key role of C1 NDUFAB1+ subtype in metabolic reprogramming, tumor progression, and immune microenvironment regulation in GC. We identified the C1 NDUFAB1+ subtype with distinct metabolic characteristics using scRNA-seq combined with ST, and investigated its unique spatial distribution patterns and interaction networks with stromal cells in the TME. Furthermore, we successfully established a risk scoring model based on the C1 NDUFAB1+ subtype. These discoveries open new avenues for the development of combined therapies targeting immune evasion and metabolic regulation, while also offering new theoretical foundations and potential targets for the precision treatment of GC.
Materials and methods
Acquisition of GC data
In this study aimed at scRNA-seq analysis of GC, we accessed the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the accession number GSE183904. Subsequently, we selected 20 GC tissue samples for analysis. Additionally, we obtained bulk RNA-seq data from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). As the datasets utilized in this research were sourced from publicly available repositories, no ethical review was required.
Single-cell sequencing data handling and interpretation
During the data import phase, we used R software (v4.2.2) and the Seurat package (v4.3.0) to perform initial handling of the single-cell sequencing data (35, 36). To ensure data reliability, we excluded low-quality cells by purifying, filtering, and removing doublets (37–39). The filtering criteria included (1): nFeature range (300–5,000) (2); nCount range (500–50,000) (3); mitochondrial gene expression proportion (≤25%) (4); red blood cell gene expression proportion (≤5%). After strict filtering, a total of 92,566 high-quality cells were preserved for subsequent study.
Data normalization was a critical step in the analysis. We used the “NormalizeData” function (40–42) to normalize the samples. To capture heterogeneity among cells, we employed the “FindVariableFeatures” function (43–45) to identify the top 2,000 most variable genes (46–48). Subsequently, we standardized the gene expression data using the “ScaleData” function (49–51) and performed Principal Component Analysis (PCA) (41, 52, 53). Batch effects were corrected using the Harmony R package (v0.1.1). Ultimately, dimensionality reduction and visualization were achieved by applying Uniform Manifold Approximation and Projection (UMAP) based on the top 30 principal components, further revealing differences in gene expression among cells (54–56).
Identification and annotation of cell types
We utilized the “FindClusters” and “FindNeighbors” functions (57) to perform clustering analysis on the cells, which initially determined the distribution of cell subtypes. To further dissect the heterogeneity of GC cells, we employed the “FindAllMarkers” function to detect genes with differential expression patterns within each subtype and annotated the cell subtypes by calculating the average expression of marker genes. The selection of marker genes integrated known information from the CellMarker database (http://xteam.xbio.top/CellMarker/) and relevant data from the literature, with manual curation ensuring the accuracy of the annotations.
Copy number variation assessment
We employed the inferCNV R program (v1.6.0) (58) to analyze scRNA-seq data, distinguishing malignant tumor cells from normal cells by calculating CNV levels. Using endothelial cells (ECs) as a reference, we effectively identified chromosomal CNVs in tumor cells. By comparing gene expression intensities between tumor cells and normal cells, inferCNV inferred amplifications or deletions in chromosomal regions.
Analysis of spatial transcriptome
We collected ST 1 and ST 2 slides from the ST dataset (GSE251950) and employed the Robust Cell Type Decomposition (RCTD) technique to map cell types identified in the scRNA-seq dataset to the spatial ST data, aiming to reveal the distribution of different cell types across spatial regions. First, we used the “FindAllMarkers” function to identify marker genes for each cell type, filtering for markers with a positive log2-transformed fold change. Subsequently, following the standard RCTD analysis pipeline, we analyzed the reference scRNA-seq data and Visium ST data in a complete doublet mode. To further determine the spatial distribution of cell types, we localized the cell types identified in the scRNA-seq data to spatial regions significantly enriched by Multimodal Intersection Analysis (MIA). Additionally, we utilized the “StLearn” package in Python to explore interactions in detail within spatial regions.
Biological function enrichment analysis
Using the “ClusterProfiler” package, we executed functional profiling on differentially expressed genes (DEGs) in tumor cell subtypes through Gene Ontology (GO) (59–62) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (63–65). Besides, Gene Set Enrichment Analysis (GSEA) was employed to evaluate the overall expression patterns of gene sets, further revealing the functional characteristics of tumor cell subtypes. To more comprehensively analyze the variability in gene expression data, we applied Gene Set Variation Analysis (GSVA), computing the enrichment scores of specific gene sets across individual samples.
Differentiation trajectory and stemness analysis
Based on the Slingshot (v2.6.0) trajectory inference framework, we constructed a developmental trajectory map of tumor cell subtypes and revealed the temporal characteristics of their differentiation pathways through pseudotime analysis. Smooth trajectory curves were derived using the “getLineages” and “getCurves” tools, and changes in cell expression levels were visualized in reduced-dimensional space. The initial and final points of the differentiation trajectories were determined based on the biological characteristics of GC tumor cell subtypes and the expression patterns of their marker genes. Additionally, we employed CytoTRACE to evaluate the stemness levels of different tumor cell subtypes, with scores ranging from 0 to 1, where higher scores indicated stronger cell stemness.
Cell communication
We leveraged the “CellChat” (v1.6.1) (66) for cell type interaction analysis in the GC microenvironment, constructing a cell communication network. The “netVisual_diffInteraction” function was used to visualize the changes in communication intensity between cells, while the “identifyCommunicationPatterns” function was employed to identify complex communication patterns. By integrating the CellChatDB database, we identified key signaling pathways and ligand-receptor pairs. A P-value cutoff of 0.05 was established in the analysis to assess the statistical significance of cell-cell interactions.
Transcriptional regulatory network analysis
We used the pySCENIC package (v0.10.0) (67) to reconstruct the gene regulatory network of tumor cells in GC, identifying stable cell states and assessing transcriptional activity. AUCell (68) was employed to compute the activity of regulons, with a score threshold set at 0.2, and significance was evaluated through 100 random permutations. Key TFs were selected based on their regulon specificity scores and AUCell activity.
Construction and validation of a prognostic model for GC
First, we extracted pertinent data from the TCGA database and identified eight genes with significant prognostic relevance through univariate Cox regression analysis (69–72). To mitigate multicollinearity among the genes, LASSO (73, 74) and multivariate Cox regression analysis was subsequently employed, ultimately identifying seven key prognostic genes and calculating their coefficient values. A risk scoring model was developed using the following formula: Risk Score = , with X denoting the coefficient and Y denoting the gene expression. Using the optimal threshold identified by the “surv_cutpoint” algorithm, patients were separated into high- and low-risk groups. To validate the predictive ability of the model, we utilized Kaplan-Meier and Receiver Operating Characteristic (ROC) curves to assess the model’s predictive accuracy (75–79). Additionally, we constructed a nomogram incorporating age, gender, race, stage, risk score, and TNM classification (80). The model was further validated through Area Under the Curve (AUC) values (81–83) based on both the C-index and true positive rate.
Analysis of the immune landscape and therapeutic response
The CIBERSORT and ESTIMATE (84, 85) were used to assess immune infiltration. Subsequently, we utilized the Tumor Immune Dysfunction and Exclusion (TIDE) program to evaluate the response to immunotherapy in the high- and low-risk groups. Afterward, we employed the “pRRophetic” package (v0.5) to estimate the half-maximal inhibitory concentration (IC50) values for drug.
Cell culture
The NCI-N87 GC cell line was cultured in 1640 basal medium enriched with 10% fetal bovine serum (FBS), 1% penicillin-streptomycin (P/S), 1% GlutaMAX-1 glutamine, and 1% sodium pyruvate, while the AGS GC cell line was cultured in F12K medium supplemented with 10% premium FBS and 1% penicillin-streptomycin (P/S). Both GC cell lines were maintained under standard conditions (37°C, 5% CO2, 95% humidity).
Cell transfection
In the cell transfection experiment, we utilized small interfering RNAs (siRNAs) provided by GenePharma (Suzhou, China) to inhibit ELK4 expression. Cells were grown in 6-well plates until achieving 50% confluency, followed by transfection using Lipofectamine 3000 RNAiMAX (Invitrogen, USA). Two specific siRNA sequences, siELK4-1 (CAUUCAACAUGAUUGCAUU) and siELK4-2 (CUCAGAUACUAUUAUGUAA), were designed for the experiment, along with a negative control siRNA (si-NC) as a control.
Cell viability assay
We employed the CCK-8 method to evaluate the viability of transfected NCI-N87 and AGS cells. Cells were plated in 96-well plates at a concentration of 5×10³ cells per well and incubated for 24 h to allow proper adherence. Subsequently, 10 μL of CCK-8 reagent (A311-01, Vazyme) was added to each well, gently mixed, and incubated at 37°C in the dark for 2 h. Following incubation, absorbance readings were taken at 450 nm using a microplate spectrophotometer (A33978, Thermo). The experiment was conducted daily from day 1 to day 4 post-transfection, and the absorbance values were recorded and averaged. A curve depicting changes in cell viability over time was plotted to visually reflect the dynamic changes in cell viability levels.
Quantitative real-time polymerase chain reaction
We first extracted RNA from the cells using TRIzol reagent, followed by reverse transcription. Then, the amplification process was monitored in real-time using fluorescent dyes. The primer sequences used in the experiment were: F: GGATTCGCAAGAACAAGCCT, R: TCAATCCTGCCCACTGTCAT.
5-ethynyl-2’-deoxyuridine proliferation assay
We seeded the transfected NCI-N87 and AGS cells in 6-well plates at a density of 5×10³ cells per well and cultured them for 24 h to allow proper adherence. Subsequently, 2× EDU working solution was added to the medium, and the cells were incubated at 37°C for 2 h. Post-incubation, the supernatant was aspirated, and the adherent cells were subjected to dual PBS washes for the removal of non-incorporated reagent. Following this step, the cellular specimens underwent fixation in 4% paraformaldehyde for a duration of 0.5 h. Afterward, a permeabilization procedure was performed with a mixture of 2 mg/mL glycine and 0.5% Triton X-100 to optimize staining performance.
To conclude the procedure, the cells were treated with a 1:1 mixture of 1ml 1X Apollo and Hoechst 33342 solution for 0.5 h to label proliferating cells and nuclei. After staining, images were observed and captured, and the frequency of EDU-incorporated cells was quantified and statistically evaluated.
Wound healing assay
Transfected cellular populations were initially seeded into 6-well tissue culture plates and allowed them to grow until they achieved nearly 95% confluency. Uniform linear wounds were generated using a sterile 200 μL pipette tip on the cell monolayer to simulate wounds. Post-scratching, the wells underwent additional washing steps, and serum-free medium was applied to prevent serum interference with cell migration. Following this, images of the wound areas were then acquired at the initial time point (0 h) and after 48 h to track the migratory behavior of the cells over time. After the experiment, the scratch widths were quantitatively analyzed using Image-J software to calculate the cell migration rate.
Transwell assay
After a 1-day serum starvation period, Matrigel (BD Biosciences, USA) was thoroughly mixed with the cell suspension and applied to the upper compartment of the transwell, with serum-added medium being introduced into the lower chamber. Upon completion of a 2-day incubation in culture dishes, cellular specimens were fixed using 4% paraformaldehyde. Finally, crystal violet staining was applied to quantitatively assess migratory and invasive capabilities.
Statistical analysis
The research employed R and Python computational tools for statistical computations and data interpretation, with the Wilcoxon rank-sum test and Spearman correlation methodology being the principal approaches for evaluating intergroup statistical significance. For the significance determination criteria, a two-tailed test was used to calculate P-values, with a statistical threshold set at 0.05. The statistical significance levels were denoted using asterisks, with * representing P-values below 0.05, ** indicating P < 0.01, *** signifying P < 0.001, and **** marking P-values less than 0.0001. This hierarchical system enhanced the interpretability of the results by quantifying the strength of differences.
Results
Single-cell profiling uncovered heterogeneity and subtype-specific molecular signatures in GC
Based on GC scRNA-seq data, we thoroughly analyzed the cellular composition and functional characteristics of its microenvironment. Our analysis workflow was shown in Figure 1. We first performed quality control and batch effect correction on the collected GC tissues. Then, after dimensionality reduction and clustering of the high-quality filtered cells, we identified 10 major cell types, including epithelial cells (EPCs), ECs, fibroblasts, myeloid cells, pericytes, mast cells (MCs), T cells and NK cells, B cells, plasma cells, and proliferating cells. The distribution profiles of 20 individual samples were additionally presented in Figure 2A, stratified by cell cycle phases (G1, G2/M, S) and group types (intestinal, diffuse). The heterogeneity of EPCs in the TME is a complex and significant phenomenon, profoundly influencing tumor development and treatment response. GC, a cancerous growth arising from the stomach mucosal epithelium, is closely associated with EPCs in its occurrence and progression (86). For the discrimination of cancerous cells from normal cellular populations and examination of tumor heterogeneity, we used ECs as a comparative baseline and identified tumor cells within EPCs through inferCNV analysis (Supplementary Figure 1). Following this, we performed a more detailed examination regarding the tumor cells and classified four tumor cell subtypes based on the expression levels of marker genes (C0 MUC5AC+ subtype, C1 NDUFAB1+ subtype, C2 SRGN+ subtype, C3 HEPACAM2+ subtype). In Figure 2B, we demonstrated the clustering of these four tumor cell subtypes and the expression distribution of CNVscore, Cell-Stemness-AUC, nFeature-RNA, and nCount-RNA across all tumor cells. Additionally, we displayed the expression of the top 10 marker genes, the distribution of named genes, as well as the significantly upregulated and downregulated genes in each subtype (Figures 2C–E). Interestingly, we observed that C1 NDUFAB1+ subtype exhibited higher expression levels of nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score in Figure 2F. Therefore, we speculated that C1 NDUFAB1+ subtype might be in a more active cell cycle state and could potentially possess stronger malignant characteristics. In Figure 2G, we depicted the proportion of different samples for each subtype in the two groups, finding that the percentage of intestinal samples in C1 NDUFAB1+ subtype was higher than that of diffuse samples. As expected, the cell cycle phases of C1 NDUFAB1+ subtype showed a preference for the G2/M and S phases (Figure 2H).

Figure 1. GC scRNA-seq: analysis workflow. The analysis process of GC scRNA-seq data covered data normalization, identification of key cell subtypes, and functional interpretation.

Figure 2. Detailed single-cell profiling of GC. (A) UMAP plot in the upper left corner displayed the distribution of 10 different cell types, while the UMAP plots in the lower right corner displayed the distribution of 20 different samples, cell cycle phases, and groups, respectively. (B) The circular plot showed the clustering of four tumor cell subtypes identified in GC, outlined by contour curves. The outer, middle, and inner axes represented the logarithmic scale of the clusters for each tumor cell subtype, along with group and phase. UMAP plots were arranged in a clockwise direction starting from the top left corner, displaying the expression distributions of CNVscore, Cell-Stemness-AUC, nFeature-RNA, and nCount-RNA across all tumor cells. (C) Bubble plot displayed the differential expression of the top 10 marker genes across the four tumor cell subtypes and two groups. The pie charts displayed the proportions of G1, G2/M, and S phases, while the violin plots presented the expression levels of G2/M.Score, S.Score, and nCount-RNA. The size of the bubbles indicated the percentage of gene expression, and the color intensity represented the level of gene expression. (D) UMAP plots showcased the distribution of named genes for each tumor cell subtype. (E) Volcano plots highlighted the differentially upregulated and downregulated genes in each tumor cell subtype. (F) Violin plots illustrated the expression levels of nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score across different tumor cell subtypes. (G) Box plots described the proportion of different samples in each subtype across the two groups. (H) Heatmap assessed phase preference for each subtype using the Ro/e score.
Spatial multi-omics features and co-localization of key subtypes in GC
To investigate the spatial multi-omics characteristics of GC and explore the spatial expression patterns of key molecular events during tumor progression, we conducted ST analysis on two collected GC tissue sections. We utilized the RCTD deconvolution method to display the first cell types inferred at selected points on ST 1 slide in Figure 3A. Figure 3B illustrated the ST landscape of C1 NDUFAB1+ subtype, which was consistent with the first cell types inferred above. To validate the accuracy of RCTD, we further analyzed the data using the MIA method (Figures 3C, D). Analysis revealed that the C7 cluster showed maximal enrichment within C1 NDUFAB1+ subtype, and the spatial distribution characteristics of the C7 cluster on ST 1 slide were similar to those of C1 NDUFAB1+ subtype. In addition, we performed RCTD analysis on ST 2 slide, which revealed high expression in regions spatially corresponding to C1 NDUFAB1+ subtype (Supplementary Figures 2A, B). Furthermore, Supplementary Figure 2C demonstrated that nCount-Spatial, nFeature-Spatial, G2/M.Score, and S.Score exhibited spatial expression patterns similar to those of C1 NDUFAB1+ subtype.

Figure 3. Spatial distribution features and enrichment analysis of tumor cell subtypes in GC. (A) The ST feature map demonstrated the first cell types inferred at selected points on ST 1 slide. Each point represented the cell type with the highest probability within that location. (B) The ST feature map revealed the spatial expression pattern of the C1 NDUFAB1+ subtype. The intensity of the color indicated the relative strength of expression. (C) MIA analysis assessed the enrichment degree of spatial clusters associated with various cell types on ST 1 slide. (D) The results of spatial spot clustering on ST 1 slide were visualized. (E) The word cloud diagrams depicted the main biological processes of each tumor cell subtype. (F) Heatmap illustrated the top 5 enriched GOBP and KEGG terms in tumor cell subtype. (G) Heatmaps respectively showed the top 20 enriched metabolism-related pathways across all tumor cell subtypes and C1 NDUFAB1+ subtype.
Biological functions and metabolic analysis of tumor cell subtypes in GC
These characteristics of C1 NDUFAB1+ subtype suggested that they might be a key driver subtype in tumor progression. Consequently, we proceeded to investigate their biological functions. Figures 3E, F illustrated the main biological processes enriched by DEGs in different tumor cell subtypes. We found that C0 MUC5AC+ subtype was primarily associated with junction, adhesion, and splicing, and were mainly enriched in pathways such as cell-cell junction organization, RNA splicing, miRNA metabolic process, regulation of miRNA metabolic process, and DNA-templated transcription elongation. C1 NDUFAB1+ subtype was mainly associated with mitochondrial, localization, and triphosphate, and were predominantly concentrated in pathways such as cytoplasmic translation, oxidative phosphorylation, aerobic respiration, ATP synthesis coupled electron transport, and mitochondrial ATP synthesis coupled electron transport. C2 SRGN+ subtype was primarily associated with leukocyte, production, and immune, and were mainly enriched in pathways such as leukocyte cell-cell adhesion, regulation of T cell activation, regulation of leukocyte cell-cell adhesion, regulation of cell-cell adhesion, and positive regulation of leukocyte cell-cell adhesion. C3 HEPACAM2+ subtype was mainly associated with peptide, secretion, and splicing, and demonstrating major representation in functional pathways such as protein folding, regulation of mRNA splicing and via spliceosome, response to unfolded protein, regulation of RNA splicing, and response to glucose. Meanwhile, in Figure 3F, we illustrated the KEGG terms enriched in each tumor cell subtype. C0 MUC5AC+ subtype was primarily associated with pathogenic Escherichia coli infection. C1 NDUFAB1+ subtype showed significant enrichment in ribosome, Parkinson disease, Huntington disease, prion disease, and oxidative phosphorylation. C2 SRGN+ subtype was predominantly linked to Th17 cell differentiation, Th1 and Th2 cell differentiation, type I diabetes mellitus, graft-versus-host disease, and allograft rejection. Finally, C3 HEPACAM2+ subtype was mainly related to protein processing in endoplasmic reticulum, vibrio cholerae infection, phagosome and protein export. Subsequently, in Figure 3G, we displayed the top 20 enriched metabolism-related pathways across different tumor cell subtypes. Notably, compared to other subtypes, the C1 NDUFAB1+ subtype showed predominant enrichment in oxidative phosphorylation, glutathione metabolism, glycolysis/gluconeogenesis, sulfur metabolism, and citrate cycle (TCA cycle). The above results indicated that C1 NDUFAB1+ subtype exhibited high metabolic activity, a characteristic that likely supported their rapid proliferation and survival, making them a key driver subtype in tumor progression. Meanwhile, C0 MUC5AC+ subtype, C2 SRGN+ subtype, and C3 HEPACAM2+ subtype played important roles in cell-cell connections, immune responses, and stress responses, respectively, highlighting the heterogeneity and complexity of the tumor.
Differentiation trajectories and stemness analysis
The differentiation process of tumor cells is one of the pivotal elements in tumorigenesis and progression, and stemness, as an important characteristic of tumor cells, directly influences their self-renewal, proliferation, and differentiation capabilities (87). Understanding the differentiation trajectories and stemness features of different tumor cell subtypes is crucial for revealing tumor heterogeneity, prognosis, and potential therapeutic targets. Therefore, we systematically explored the differentiation lineages and related characteristics of different tumor cell subtypes using a combination of analytical methods. First, we ranked the stemness of different tumor cell subtypes using the CytoTRACE method, and the observations supported that C1 NDUFAB1+ subtype possessed increased stemness properties, while C0 MUC5AC+ subtype exhibited lower stemness (Figures 4A, B). This finding suggested that C1 NDUFAB1+ subtype might exhibit enhanced capacity for self-renewal and differentiation. To further investigate their stemness features, we displayed the variations in the expression of stemness genes among distinct subtypes (Figure 4C) and visualized the distribution of significantly expressed stemness genes (CTNNB1, ABCG2, MYC, LGR5) in C1 NDUFAB1+ subtype (Figure 4D). The high expression patterns of these genes further supported the high stemness characteristics of C1 NDUFAB1+ subtype. To comprehensively understand the differentiation dynamics of tumor cells, we used Slingshot to evaluate the differentiation lineages of four subtypes (Figure 4E). The results showed that both lineage 1 and lineage 2 passed through C0 MUC5AC+ subtype in the early differentiation stage and diverged to different endpoints after passing through C1 NDUFAB1+ subtype. This observation indicated that C1 NDUFAB1+ subtype might serve as critical factors in cellular differentiation mechanisms. To further validate this hypothesis, we displayed the distribution of lineage 1 and lineage 2 fitted by Slingshot, clearly showing their progression along the inferred pseudotime order (Figure 4F). Additionally, the pseudotime trajectories of named genes in the four tumor cell subtypes within lineage 1 and lineage 2 were highly consistent with the Slingshot results (Figure 4G), further supporting our findings.

Figure 4. Developmental and differentiation features of tumor cell subtypes in GC. (A) The two-dimensional graphs depicted the CytoTRACE results of the predicted order (left) and distribution characteristics (right) for the four tumor cell subtypes. (B) Box plot described the stemness ranking of the tumor cell subtypes predicted by CytoTRACE. (C) Bubble plot displayed the differential expression of the top stemness genes for each tumor cell subtype. (D) UMAP plots visualized the distribution of significantly expressed stemness genes in C1 NDUFAB1+ subtype. (E) UMAP plot utilized Slingshot analysis to display two cellular differentiation lineages of the four tumor cell subtypes over time. Solid lines and arrows represented the differentiation trajectories and the direction of cell differentiation from naive to mature, respectively. (F) UMAP plots illustrated the distribution of lineage 1 and lineage 2 fitted by Slingshot, showing their progression along the inferred pseudotime order. (G) Dynamic trend plots depicted the pseudotime trajectories of named genes across four tumor cell subtypes in lineage 1 and lineage 2. (H) Heatmap demonstrated the key biological processes related to the two lineages in tumor cell differentiation, as presented in the GO enrichment analysis results. The ridge plots displayed the pseudotime density changes of the four tumor cell subtypes. The trajectory plots showed the pseudotime score changes of S.Score and G2/M.Score across the two lineages.
Based on these observations, we speculated that C1 NDUFAB1+ subtype might play a critical role in the differentiation process of GC cells. Specifically, C0 MUC5AC+ subtype exhibited lower stemness and higher differentiation levels, while C1 NDUFAB1+ subtype demonstrated higher stemness, proliferative capacity, and differentiation potential. To more deeply reveal the biological significance of these different lineages, we explored key biological processes associated with tumor cell differentiation (Figure 4H). The evidence suggested that C1 cluster was mainly tied to immune and humoral responses, C2 cluster was related to digestive and structural functions, and C3 cluster was closely linked to Wnt signaling pathway and morphogenesis. Additionally, C4 cluster was significantly associated with visual, brain, and pancreatic functions. These enrichment results not only revealed the potential functional differences among different tumor cell subtypes but also provided new insights into their roles in tumorigenesis and development.
Cell-cell communication and signaling pathway regulatory network analysis
To better elucidate the involvement of C1 NDUFAB1+ subtype in the TME, we investigated the complex interactions between tumor cells and other cell types. Using CellChat, we executed a detailed investigation of the communication networks between C1 NDUFAB1+ subtype as both source and target with other cell types (Figure 5A). We found that the interactions between C1 NDUFAB1+ subtype and fibroblasts, as well as pericytes, were stronger compared to other cell types. We examined the main driving patterns of all cell types and the interaction proteins under these patterns (Figures 5B, C). The investigations indicated that the outgoing and incoming signals of C1 NDUFAB1+ subtype were primarily driven by Pattern 3, involving significant contributions from CDH, EDN, OCLN, DESMOSOME, and EPHA. In Figure 5D, we further presented the relative signal strength of various signaling pathways across different cell types and all cell types under different signal patterns. Through signal network analysis, we discovered that C1 NDUFAB1+ subtype primarily communicated with fibroblasts and pericytes via the PARs signaling pathway (Figure 5E). Additionally, in Figure 5F, we observed that C1 NDUFAB1+ subtype mainly took on the responsibilities of sender, mediator, and influencer, while fibroblasts and pericytes acted as receivers. Subsequently, we used a hierarchical graph to reveal that C1 NDUFAB1+ subtype might regulate themselves through autocrine mechanisms, while influencing fibroblasts and pericytes through paracrine mechanisms (Figure 5G). Furthermore, we visualized the expression of key ligand-receptor pairs in the PARs signaling pathway across four tumor cell subtypes and nine other cell types (Figures 5H, I). We found that PRSS3 and F2R were highly expressed in C1 NDUFAB1+ subtype, fibroblasts, and pericytes. We then further analyzed the communication network within the PRSS3-F2R ligand-receptor pair, confirming that the interactions between C1 NDUFAB1+ subtype and fibroblasts and pericytes could be mediated by PRSS3-F2R within the PARs signaling pathway (Figures 5J, K).

Figure 5. Cell communication in single-cell transcriptomics. (A) Circle diagrams represented the interactions between C1 NDUFAB1+ subtype as both the source (upper) and the target (lower) with other cell types in terms of the number (left) and intensity (right). (B) Heatmaps showed the outcoming (left) and ingoing (right) communication patterns of the four tumor cell subtypes and nine other cell types across three cell communication patterns, along with the contributions of different proteins in each communication pattern. (C) Sankey diagrams predicted the outgoing communication patterns (upper) of the four tumor cell subtypes and nine other cell types as secreting cells, and the incoming communication patterns (lower) as target cells, along with the signaling pathways under the three cell communication patterns. (D) Heatmaps displayed the relative intensity of various signaling pathways in the four tumor cell subtypes and nine other cell types under the outgoing and incoming signaling patterns, while the bar charts illustrated the relative signal intensity of different cell types. (E) Heatmap showed the communication probability of the four tumor cell subtypes and nine other cell types under the PARs signaling network. (F) Heatmap depicted the centrality scores within the PARs signaling pathway network. (G) Hierarchical diagram portrayed the interactions between four tumor cell subtypes and nine other cell types in the PARs signaling pathways network. (H, I) Bubble plot and violin plot compared the expression of significant ligands and receptors in the PARs signaling pathways across four tumor cell subtypes and nine other cell types. (J, K) Chord plot and circle diagram displayed the communication network in the PRSS3-F2R ligand-receptor pair.
Integrating single-cell and ST to analyze intercellular spatial interactions
To further understand the spatial distribution and spatial interactions of different cell types within the TME, we integrated scRNA-seq data with ST data using the RCTD deconvolution technique, delineating the spatial architecture of different cell types on ST 2 slide (Figure 6A). Subsequently, we analyzed ST 2 slide using the “Stlearn” package in Python, displaying the top 10 significant ligand-receptor pairs within the spots in Figure 6B. We illustrated the spatial interaction strength and statistical value of the THBS2-ITGB1 interaction pair and performed spatial enrichment analysis. By combining the results with Figure 6A, we observed that the THBS2-ITGB1 interactions were predominantly concentrated at the boundaries of tumor cells. This suggested that the ligand-receptor interactions between THBS2 and ITGB1 might play a regulatory role in spatial communication between tumor cells and other cell types (Figure 6C). In Figure 6D, we visualized the strength of cell-cell interactions mediated by the THBS2-ITGB1 ligand-receptor pair, revealing higher communication intensity between C1 NDUFAB1+ subtype and fibroblasts as well as pericytes. The results in Figures 6E, F further confirmed this observation and unveiled the complexity of this spatial cell-cell communication pattern. Based on these findings, we speculated that C1 NDUFAB1+ subtype might have reshaped the network of intercellular interactions within the TME through abnormal communication crosstalk with fibroblasts and pericytes.

Figure 6. Spatial cell communication. (A) RCTD analysis predicted the cell types at each spatial spot on ST 2 slide. (B) The ranking plot displayed the top 10 significant ligand-receptor pairs within the spots. (C) The interaction strength of the THBS2-ITGB1 ligand-receptor pair was represented across all spots (left), in significant spots (middle), and its statistical value was shown for each spot (right). (D) The interaction heatmap visualized the intensity of intercellular interactions mediated by the THBS2-ITGB1 ligand-receptor pair. (E, F) Chord plot and network diagram depicted the spatial interactions among different cell types in the THBS2-ITGB1 ligand-receptor pair.
Investigation of C1 NDUFAB1+ subtype based on TF regulatory networks
Investigating TFs is crucial for exploring tumor heterogeneity and the TME, as TFs play a pivotal role in regulating gene expression and cellular behavior. Figure 7A displayed the clustering analysis of tumor cells based on their gene expression levels, followed by the clustering patterns of different tumor cell subtypes according to regulon activity scores (Figure 7B). Subsequently, we identified three regulatory modules through hierarchical clustering of TFs within the tumor cell subtypes (Figure 7C). We visualized the expression distribution of TFs across each module and the expression levels of different tumor cell subtypes within each module (Figures 7D, E), revealing that C1 NDUFAB1+ subtype displayed elevated expression within the M3 module relative to other subtypes. Meanwhile, within the M3 module, C1 NDUFAB1+ subtype showed higher regulon activity scores relative to other tumor cell subtypes (Figure 7F). Figure 7G illustrated the top 5 TFs in C1 NDUFAB1+ subtype, namely ZNF615, TFDP1, E2F1, ETV4, and ELK4. In Figure 7H, we observed the top 5 TFs with higher specificity scores in C1 NDUFAB1+ subtype. We further analyzed the expression levels of the top 5 TFs in C1 NDUFAB1+ subtype across different subtypes (Figure 7I) and visualized their expression distribution within the M3 module using UMAP plots (Figure 7J). We noted that ELK4 exhibited elevated expression levels relative to other subtypes, and its expression within the M3 module was more pronounced relative to the other four TFs. Therefore, we conducted in vitro experiments to further investigate the effects of ELK4 on GC cells.

Figure 7. Characterization of TF activity and regulatory modules in tumor cell subtypes. (A) UMAP plot revealed the clustering patterns of tumor cells, which were determined by gene expression levels. (B) UMAP plots colored and visualized distinct clustering patterns among tumor cells based on the regulon activity scores. (C) Heatmap demonstrated the identification of three regulatory modules through hierarchical clustering of TFs within tumor cell subtypes. (D) UMAP plots presented the distribution of expression levels of TFs in each module. (E) Violin plots illustrated the expression levels of four tumor cell subtypes within each module. (F) Scatter plots demonstrated the ranking of regulon activity scores of TFs across four tumor cell subtypes within each module. (G) Heatmap exhibited the expression of the top 5 TFs in each tumor cell subtype. (H) Scatter plots demonstrated the ranking of specificity scores for the top 5 TFs in each tumor cell subtype. (I) Violin plots illustrated the expression levels of the top 5 TFs in C1 NDUFAB1+ subtype across each tumor cell subtype. (J) UMAP plots visualized the expression distribution of the top 5 TFs in C1 NDUFAB1+ subtype within the M3 module.
In vitro experiments validated the regulatory role of ELK4 throughout the advancement of GC
To gain deeper insights into the involvement of ELK4 in GC advancement, in vitro assay was carried out with NCI-N87 and AGS cell populations. Originally, we performed knockdown of ELK4 and then categorized into three distinct groups: si-NC, siELK4-1, and siELK4-2. Subsequently, we observed that the relative expression levels of both mRNA and protein were markedly decreased in the siELK4–1 and siELK4–2 groups when measured against the si-NC group (Figure 8A). Additionally, using the CCK-8 assay, we observed that the mean optical density (OD) values were also significantly decreased in the siELK4–1 and siELK4–2 groups (Figure 8B). This demonstrated that the cell viability of the NCI-N87 and AGS GC cell lines was significantly reduced after ELK4 knockdown. Next, colony formation and EDU staining results revealed that the colony numbers and the cell proliferation rates were markedly decreased in both GC cell populations after ELK4 knockdown (Figures 8C–E). Finally, the cell wound healing and transwell assays showed that the wound healing ability, migration, and invasion capabilities of the two GC cell lines were inhibited following ELK4 knockdown (Figures 8F–H). Through these studies, we uncovered ELK4’s crucial involvement in GC development, providing important experimental evidence for further exploration of ELK4 as a promising target for GC treatment.

Figure 8. In vitro experimental validation. (A) Bar plots displayed the relative expression levels of ELK4 mRNA and ELK4 protein in the NCI-N87 and AGS GC cell lines for si-NC, siELK4-1, and siELK4-2. (B) The CCK-8 assay demonstrated the OD values of si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (C) The colony formation assays compared the results of cell colony formation among si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (D) Bar plots visually showed the colony numbers and cell proliferation rates of si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (E) The EDU staining assays exhibited the results of cell proliferation among si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (F) The cell wound healing assays evaluated the ability of cell wound healing at 0 h and 48 h among si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (G) Transwell assays evaluated the cell migration and invasion abilities among si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. (H) Bar plots compared the results of cell wound healing, migration, and invasion among si-NC, siELK4-1, and siELK4–2 in the two GC cell lines. **P < 0.01, ***P < 0.001.
Construction and validation of a risk scoring model for tumor prognosis prediction
We designed a risk scoring model to probe the molecular mechanisms of tumor prognosis and comprehensively assess patients’ prognostic risks. Initially, eight genes with notable prognostic relevance were ascertained (Figure 9A). To mitigate multicollinearity, narrowing down the selection to seven key prognostic genes (Figure 9B). In-depth profiling of these genes and their coefficients indicated that NHLH2, ATF7, ERG, CREM, and NR3C1 were associated with poor prognosis (Figures 9C, D). Using the optimal threshold for the NDUFAB1+ tumor cells risk score (NTRS), patients were stratified into high and low NTRS groups, followed by DEGs analysis. The findings indicated notable disparities in the expression of the seven genes between the two groups, with the high NTRS group exhibiting poorer prognosis (Figures 9E, F). Kaplan-Meier survival curves further validated these findings, demonstrating lower OS rates in the high NTRS group (Figure 9G). To evaluate the model’s predictive performance, ROC curves were plotted. The AUC values for 1-year, 3-year, and 5-year survival predictions all exceeded 0.6, indicating the model’s robust predictive capability (Figure 9H). PCA demonstrated notable variations in gene distribution between the two groups, with PC1 and PC2 explaining 10.32% and 8.55% of the total variance, respectively (Figure 9I). Additionally, a negative correlation between risk scores and OS was observed, further supporting the model’s reliability (Figure 9J). Further analysis of the relationships among the seven prognostic genes, risk scores, and OS showed that NHLH2, ATF7, ERG, CREM, and NR3C1 exhibited a positive correlation with risk scores but an inverse correlation with OS, suggesting their potential roles as drivers of poor prognosis. In contrast, SOX9 and E2F2 exhibited the opposite trends, indicating potential protective effects on prognosis (Figure 9K). To explore the influence of various risk factors on prognosis, the distribution of clinical features was compared between the high and low NTRS groups (Figure 9L). A nomogram incorporating age, gender, race, stage, risk score, and TNM classification was constructed, demonstrating that the NTRS risk group exhibited the most significant survival differences (Figure 9M). Finally, the model’s accuracy was validated using the C-index and true positive rate, with AUC values for 1-year, 3-year, and 5-year predictions all exceeding 0.6, further confirming the model’s robustness (Figures 9N, O). These findings provide critical theoretical and practical insights for optimizing clinical decision-making, accurately predicting patient prognosis, and improving survival outcomes.

Figure 9. Creation and validation of a prognostic risk scoring model. (A) Forest plot demonstrated the results of univariate Cox regression analysis, where HR < 1 represented protective factors and HR > 1 indicated risk factors. (B) LASSO regression analysis identified prognostic genes by determining the optimal parameters through cross-validation (upper) and generating the coefficient profile based on the optimal lambda (lower). (C) Forest plot illustrated the results of multivariate Cox regression analysis. (D) Bar plot exhibited the coefficients of the genes in the constructed model. (E) The risk curve and scatter plot contrasted the risk scores (upper) and survival outcomes (lower) over time in the high and low NTRS groups, respectively. (F) Heatmap compared the differences in expression levels of prognostic genes in the constructed model between the high and low NTRS groups. (G) Kaplan-Meier curves compared the OS rates between the high and low NTRS groups. (H) ROC curves analysis for survival prediction showed sensitivity and specificity for 1-year, 3-year, and 5-year survival. (I) Scatter plot utilized PCA to demonstrate the differences in gene distribution between the high and low NTRS groups. (J) Scatter plot combined with a linear regression line showed the relationship between risk scores and OS. (K) Heatmap and scatter plots for correlation analysis, showing the correlation coefficients among prognostic genes, OS, and risk. (L) Pie chart compared the proportional distribution of clinical characteristics, including status, age, gender, race, stage, and TNM classification, between the high and low NTRS groups. (M) Nomogram predicted 1-year, 3-year, and 5-year OS using various clinical characteristics, including age, gender, race, stage, risk score, and TNM classification. *P < 0.05, **P < 0.01, ***P < 0.001. (N) Violin plot demonstrated the C-index for 1-year, 3-year, and 5-year predictions, measured by AUC through cross-validation. (O) ROC curves illustrated the predictive performance of the model through AUC values for 1-year, 3-year, and 5-year predictions.
Investigating the effect of NTRS on the immune landscape and prognosis in tumor cells
The initial step in examining NTRS’s complex involvement in the immune microenvironment involved determining the relative abundance of distinct immune cell populations (Figure 10A). Figure 10B manifested pronounced contrasts regarding the estimated proportions of nine immune cell types between both groups. In the high NTRS group, MCs resting, Macrophages M2, B cells naive, and Monocytes were more abundant. In contrast, T cells CD4 memory activated, T cells follicular helper, MCs activated, NK cells resting, and Macrophages M0 were more prevalent in the low NTRS group. Comprehensive analyses were performed to evaluate the interconnections among diverse immune cell types and risk scores, prognostic genes, and OS (Figures 10C, D). MCs resting, Macrophages M2, Monocytes, and B cells naive showed significant positive correlations with risk scores. Conversely, T cells CD4 memory activated and T cells follicular helper exhibited significant negative correlations with risk scores. These findings suggested that the observed differences in immune cell composition across study groups might be closely related to disease risk and prognosis. Next, we compared the differences in signature and TIDE scores between the high and low NTRS groups (Figures 10E, F). The results indicated that the stromal, immune, ESTIMATE, and TIDE scores were significantly higher in the high NTRS group, suggesting a more immunosuppressive TME in this group. Additionally, analysis outcomes of immunological checkpoint-associated genes (Figure 10G) uncovered that the majority immune checkpoints were positively correlated with risk scores. Genes such as ATF7, ERG, CREM, NR3C1, and E2F2 were positively correlated with multiple immune checkpoints, while NHLH2 and SOX9 showed negative correlations. In the high NTRS group, most immune checkpoint-related genes were expressed at higher levels, with only TNFRS14 and LGALS9 showing higher expression in the low NTRS group (Figure 10H).

Figure 10. Analysis of immune profiling, enrichment, and drug sensitivity in the prognostic model. (A) Box plot displayed the estimated proportions of 22 immune cell types. (B) Box plot showcased significant differences in the estimated proportions of nine immune cell types between the high and low NTRS groups. (C) Lollipop plot revealed the results of the correlation analysis between different immune cell types and the risk score. (D) Heatmap manifested the correlation between different immune cell types and prognostic genes, OS, and risk scores. (E) Box plot compared the differences in signature scores (stromal, immune, and ESTIMATE) between the high and low NTRS groups. (F) Violin plot exhibited the comparison of TIDE score between the high and low NTRS groups. (G) Bubble plot presented the results of Spearman correlation analysis between different immune checkpoint-related genes and prognostic genes, OS, and risk scores. (H) Box plot showed significant differences in the expression of various immune checkpoint-related genes between the high and low NTRS groups. (I) Volcano plot exhibited the DEGs between the high and low NTRS groups. (J) Heatmap revealed significant differences in gene expression between the high and low NTRS groups. (K) Bar plot displayed the top 20 enriched pathways of DEGs in the KEGG enrichment analysis. (L) GSEA enrichment analysis of DEGs revealed the enrichment results of major biological processes. (M) GSVA enrichment analysis showed the main biological pathways and gene sets enriched in the high and low NTRS groups. (N) Violin plots combined with box plots illustrated the drug sensitivity differences between the high and low NTRS groups by comparing the IC50 values assessed using different therapeutic drugs. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.
To comprehensively reveal the discrepancies between both groups, we probed the upregulated and downregulated DEGs and distinct gene expression patterns (Figures 10I, J). DEGs were found to be predominantly enriched in pathways such as neuroactive ligand-receptor interaction, calcium signaling pathway, adrenergic signaling in cardiomyocytes, cAMP signaling pathway, and cytoskeleton in muscle cells (Figure 10K). GSEA enrichment analysis (Figure 10L) further demonstrated that upregulated genes were primarily engaged in biological processes such as membrane depolarization during action potential, peptide cross linking, smooth muscle contraction, and sodium ion export across plasma membrane. Downregulated genes were chiefly implicated in kinetochore organization, cell cycle DNA replication, and centromere complex assembly. GSVA analysis highlighted significant differences in major biological pathways and gene set enrichment between both groups (Figure 10M). Finally, we evaluated the sensitivity to different therapeutic drugs (Figure 10N). The high NTRS group showed greater sensitivity to AZD.2281, Imatinib, MG.132, MK.2206, Pazopanib, PF.02341066, and Sunitinib. In contrast, the low NTRS group exhibited higher sensitivity to Metformin, Methotrexate, and Paclitaxel. These findings not only revealed the complex regulatory network of NTRS within the immune landscape and its possible influence on disease prognosis but also laid a crucial foundation for developing precise treatment strategies tailored to different risk score groups.
Discussion
GC maintains a notable position in the epidemiology of digestive system cancers (88). Owing to nonspecific initial clinical manifestations, most cases are diagnosed during progressive disease phases, leading to unfavorable outcomes and a 5-year OS rate below 5% (89). Although the treatment options for GC are diverse, such as surgical removal, chemotherapy, radiation therapy, targeted treatments, and immunotherapeutic approaches, these methods still face significant limitations in practical applications. Firstly, tumor heterogeneity leads to substantial differences in treatment responses among distinct patients or even within different regions of the same tumor, greatly limiting the universality and effectiveness of treatments (90). Secondly, drug resistance is a particularly prominent issue (91), with many patients developing resistance after initial treatment, leading to disease progression (92). Additionally, the complexity of the TME further exacerbates treatment challenges, as the immunosuppressive cells within, such as tumor-associated macrophages and regulatory T cells, promote immune evasion and drug resistance in tumor cells by secreting immunosuppressive factors and metabolites (93–95). Meanwhile, GC cells adapt to microenvironmental changes through metabolic reprogramming, such as enhanced glycolysis, glutamine metabolism, and fatty acid synthesis, to maintain their survival advantages, which also provides a potential mechanism for treatment resistance (96). These factors not only limit the effectiveness of current treatments but also highlight the urgent need to develop new therapeutic strategies.
ScRNA-seq technology not only reveals the molecular characteristics of GC but also deciphers makeup and operational dynamics of various cellular entities in the TME, offering essential perspectives on tumor heterogeneity, drug resistance mechanisms, and immune evasion (97–99). In this research, we pinpointed 10 key cell types that were pivotal in the advancement of GC. Notably, the heterogeneity of EPCs, particularly the identification of four distinct tumor cell subtypes (C0 MUC5AC+, C1 NDUFAB1+, C2 SRGN+, and C3 HEPACAM2+), highlighted the complexity of GC biology. Among these tumor cells, C1 NDUFAB1+ subtype demonstrated elevated expression of genes associated with the cell cycle and indicators of cellular activity, suggesting a more aggressive phenotype. This subtype was more prevalent in intestinal-type GC samples, aligning with its potential role in driving tumor progression. We also observed that C1 NDUFAB1+ subtype showed a preference for the G2/M and S phases, further supporting their proliferative and potentially more malignant characteristics. ST analysis further confirmed these findings, revealing the spatial enrichment of C1 NDUFAB1+ subtype in specific regions of GC tissues. These findings not only supported the proliferative activity and malignant features of C1 NDUFAB1+ subtype but also revealed their spatial heterogeneity within the TME. The specific spatial localization patterns of C1 NDUFAB1+ subtype in tumor tissues were likely closely related to their functional role in tumor progression.
The investigation additionally uncovered notable distinctions in the biological process and metabolic profiles among each subtype. The results showed that C1 NDUFAB1+ subtype was primarily associated with mitochondrial function, localization, and triphosphate metabolism, and was significantly enriched in biological processes including cytoplasmic translation, oxidative phosphorylation, aerobic respiration, ATP synthesis coupled electron transport, and mitochondrial ATP synthesis coupled electron transport. These findings indicated that C1 NDUFAB1+ subtype exhibited high metabolic activity, which likely provided energy support for their rapid proliferation and survival, further solidifying their role as a critical driving subtype in tumor progression. Additionally, C1 NDUFAB1+ subtype also demonstrated significant activity in metabolic pathways encompassing oxidative phosphorylation, glutathione metabolism, glycolysis/gluconeogenesis, sulfur metabolism, and citrate cycle (TCA cycle), further emphasizing their critical role in tumor metabolic reprogramming (100, 101). In contrast, other tumor cell subtypes exhibited distinct functional characteristics. C0 MUC5AC+ subtype was mainly associated with cell junctions, RNA splicing, and miRNA metabolism, potentially playing a role in maintaining intercellular communication and transcriptional regulation (102). C2 SRGN+ subtype was significantly enriched in pathways related to leukocyte-mediated immune responses and the regulation of cell-cell adhesion, suggesting their potential role in the tumor immune microenvironment (103, 104). Meanwhile, C3 HEPACAM2+ subtype was primarily linked to protein folding and endoplasmic reticulum stress responses, possibly functioning in cellular stress responses and the maintenance of protein homeostasis (105). These findings highlighted the functional heterogeneity of GC tumor cell subtypes and revealed the diverse roles of different subtypes in tumor progression.
The differentiation trajectories and stemness analysis of tumor cell subtypes in GC provided profound insights into the cellular dynamics and functional heterogeneity within the TME. Our findings revealed that C1 NDUFAB1+ subtype displayed greater stem-like characteristics relative to other subtypes, suggesting that these cells might possess stronger self-renewal and differentiation potential (106). Stemness genes were essential for supporting the self-renewal, preservation, and differentiation potential of cancer stem cells (107, 108). The high stemness of C1 NDUFAB1+ subtype was likely supported by the elevated expression of key stemness genes. CTNNB1, a pivotal molecule in the Wnt signaling pathway, was involved in cell proliferation and stemness maintenance (109). MYC, a proto-oncogene, regulated the cell cycle and metabolism (110, 111). ABCG2, a multidrug resistance protein, was potentially associated with treatment resistance in tumor cells (112, 113). LGR5, a target gene of the Wnt signaling pathway, was commonly used to mark tumor stem cells (114). The high stemness and metabolic activity of C1 NDUFAB1+ subtype suggested that they might play a critical role in driving tumor advancement and resistance to therapy (115). The differentiation trajectories inferred through Slingshot analysis further highlighted the influential role of C1 NDUFAB1+ subtype in the developmental dynamics of GC cells. Both lineage 1 and lineage 2 originated from C0 MUC5AC+ subtype, which exhibited lower stemness and higher differentiation levels, and diverged after passing through C1 NDUFAB1+ subtype. This observation suggested that C1 NDUFAB1+ subtype might act as a transitional or progenitor-like population, driving the differentiation of tumor cells into distinct functional states. The pseudotime trajectories of named genes within these lineages further validated this hypothesis, showing a clear progression along the inferred differentiation pathways. Functional enrichment analysis of the differentiation lineages revealed cluster-specific biological processes. For instance, the C1 cluster was linked to immune and humoral responses, suggesting a potential role in modulating the tumor immune microenvironment. The C2 cluster showed a connection to digestive and structural functions, possibly reflecting its involvement in maintaining tissue architecture and function. The C3 cluster showed a strong connection to the Wnt signaling pathway (116) and morphogenesis, indicating its potential effect on regulating cell fate and tissue patterning. Finally, the C4 cluster was related to visual, brain, and pancreatic functions, hinting at a broader role in systemic tumor biology. These findings underscored the functional diversity and complexity of tumor cell subtypes in GC. The high stemness and differentiation potential of C1 NDUFAB1+ subtype, coupled with their central role in differentiation trajectories, suggested that this subtype might function as an important driver of tumor progression and heterogeneity.
Next, we carried out a detailed exploration of the cell-cell communication network of C1 NDUFAB1+ subtype within the TME using CellChat. The results revealed that C1 NDUFAB1+ subtype communicated with fibroblasts and pericytes through the PARs signaling pathway. Proteinase activated-receptors (PARs), a class of G protein-coupled receptors, have been demonstrated to play important roles in various cancers (117–119). Particularly in GC, the expression of PARs was closely associated with tumor development and progression (120, 121). The role of the PARs signaling pathway in the TME has been extensively studied, especially in the interactions between tumor cells and stromal cells (122, 123). We observed that PRSS3 and F2R were highly expressed in C1 NDUFAB1+ subtype, fibroblasts, and pericytes, suggesting that the PRSS3-F2R ligand-receptor pair might contribute significantly to mediate communication between these cells. This discovery provided new insights into the molecular mechanisms of cell-cell communication in the TME. To achieve a deeper insight into the spatial distribution and interactions of different cell types within the TME, we integrated scRNA-seq data with ST data. We found that the THBS2-ITGB1 ligand-receptor pair exhibited particularly significant interactions at the boundaries of tumor cells. This phenomenon suggested that THBS2-ITGB1 might act as an important regulator in spatial communication between tumor cells and other cell types. By visualizing the strength of cell-cell interactions facilitated by THBS2-ITGB1, we observed higher communication intensity between C1 NDUFAB1+ subtype and fibroblasts, as well as pericytes. This result further supported the hypothesis that C1 NDUFAB1+ subtype reshaped the network of intercellular interactions within the TME through abnormal signaling crosstalk. Such aberrant communication patterns may have provided favorable conditions for tumor multiplication, infiltration, and metastasis. Through integrating spatial transcriptomic data with histomorphological information, we not only revealed the spatial distribution characteristics of the C1 NDUFAB1+ subtype but also provided a more comprehensive interpretation of its biological significance in tumor progression.
In this study, we conducted a meticulous exploration of C1 NDUFAB1+ subtype through TF regulatory networks. TFs orchestrate the regulation of gene expression and cellular behavior (124, 125), making their investigation crucial for understanding tumor heterogeneity and the TME. Further research revealed that C1 NDUFAB1+ subtype exhibited higher expression levels and regulon activity scores in the M3 module, suggesting that this module might have been a significant player in the functional regulation of C1 NDUFAB1+ subtype. Within the M3 module, C1 NDUFAB1+ subtype displayed higher expression of specific TFs, particularly the top 5 TFs: ZNF615, TFDP1, E2F1, ETV4, and ELK4. Among these, the expression level of ELK4 was particularly prominent in C1 NDUFAB1+ subtype, significantly higher than in other tumor cell subtypes. This finding suggested that ELK4 potentially served as a critical factor in the biological behavior of C1 NDUFAB1+ subtype. Prior studies demonstrated that ELK4, an ETS-family TF, promoted GC progression by activating oncogenic lncRNA SNHG22, and its knockdown suppressed GC cell proliferation and invasion (126). Additionally, ELK4 was shown to drive malignant phenotypes in GC by regulating the KDM5A-PJA2-KSR1 axis (127). To validate this hypothesis, we further investigated the effects of ELK4 on the proliferation, migration, and invasion capabilities of GC cells through in vitro functional experiments.
In vitro studies indicated that silencing ELK4 markedly hindered the growth of NCI-N87 and AGS cell strains. Using CCK-8 assays, colony formation assays, and EDU staining, we observed that the viability and proliferation rates of GC cells considerably diminished after ELK4 knockdown. Additionally, wound healing and transwell assays further confirmed that the migration and invasion capabilities of GC cells were markedly suppressed following ELK4 knockdown. These results collectively implied that ELK4 served as a critical factor in modulating the oncogenic properties of GC cells.
By integrating multi-omics data and clinical prognostic information, we designed and confirmed a predictive risk scoring model using C1 NDUFAB1+ subtype. The results demonstrated that NHLH2, ATF7, ERG, CREM, and NR3C1 were associated with poor prognosis, while SOX9 and E2F2 exhibited potential protective effects. Using the optimal threshold for NTRS, patients were grouped into high and low NTRS cohorts, where the high NTRS cohort showed significantly poorer survival rates. We thoroughly investigated the impact of C1 NDUFAB1+ subtype on the immune microenvironment, particularly their roles in immune escape and metabolic reprogramming. In the high NTRS group, the proportions of MCs resting, Macrophages M2, B cells naive, and Monocytes were higher, which are typically associated with an immunosuppressive microenvironment (128–131). In contrast, the low NTRS group exhibited higher proportions of T cells CD4 memory activated, T cells follicular helper, MCs activated, NK cells resting, and Macrophages M0, which are generally linked to an immune-activated state (132–135). These differences in immune cell distribution may have reflected the ability of the high NTRS group to evade immune surveillance through immune escape mechanisms, thereby promoting tumor progression. Further analysis revealed that MCs resting, Macrophages M2, Monocytes, and B cells naive were significantly positively correlated with the risk score, while T cells CD4 memory activated and T cells follicular helper were negatively correlated with the risk score. These evidences reflected that the TME in the high NTRS group might be more immunosuppressive, providing favorable conditions for tumor cell immune escape. Additionally, the high NTRS group showed significantly higher stromal, immune, ESTIMATE, and TIDE scores, further supporting the hypothesis of an immunosuppressive microenvironment in the high NTRS group. The elevated TIDE score, often associated with enhanced immune escape capabilities of tumor cells (136, 137), indicated that the high NTRS group might evade immune system attacks through multiple mechanisms. During this analysis, we observed a positive correlation between SOX9 expression and patient prognosis, but a negative correlation with immune checkpoint markers. This suggested that the primary mechanism by which SOX9 improved patient outcomes might be independent of or only weakly related to the immune system, instead operating through direct effects on tumor cell-intrinsic biological behaviors. Future studies should further investigate the precise molecular mechanisms through which SOX9 regulates malignant phenotypes in GC cells and the immune microenvironment. These findings indicate that clinical practice should comprehensively consider both the tumor biological characteristics and immune microenvironment features reflected by SOX9 expression levels to develop more precise treatment strategies.
KEGG and GSEA enrichment analyses revealed striking differences in metabolic pathways between both groups. The high NTRS group’s DEGs were primarily mapped to metabolic pathways such as neuroactive ligand-receptor interaction, calcium signaling pathways, and cAMP signaling pathways. The activation of these pathways might have promoted metabolic reprogramming in tumor cells (138–140), enabling them to adapt to nutrient deprivation and hypoxic conditions in the TME, thereby sustaining their proliferation and survival. Notably, these metabolic reprogramming features were highly consistent with the metabolic activity of the previously mentioned C1 NDUFAB1+ subtype. This metabolic reprogramming not only provided energy for tumor cells but also supplied precursor molecules for biosynthesis, supporting the maintenance of their malignant phenotype (141). Finally, we conducted extensive drug sensitivity screening and found that the high NTRS group was more sensitive to drugs such as AZD.2281, MG.132, and Pazopanib. These findings helped identify novel potential therapeutic targets, elucidate resistance mechanisms, and provide preliminary clues for future drug repositioning or combination therapy research.
This study revealed that the unique metabolic characteristics of the C1 NDUFAB1+ subtype not only supported its malignant phenotype but also potentially mediated treatment resistance. The concurrent activation of both oxidative phosphorylation and glycolysis pathways in this subtype enabled it to switch energy production modes in response to microenvironmental stresses, which likely represented a key mechanism underlying the failure of conventional chemotherapy and targeted therapies. Particularly noteworthy was the significant activation of the glutathione metabolic pathway, which might have further promoted drug resistance by scavenging reactive oxygen species and protecting tumor cells from chemotherapeutic damage (142, 143). The hyperactive TCA cycle and mitochondrial function in the C1 NDUFAB1+ subtype suggested its potential adaptive resistance to glycolysis-targeted therapies. Furthermore, this study revealed that high-NTRS patients exhibited significant immunosuppressive microenvironment characteristics, manifested by increased infiltration of immunosuppressive cells such as Macrophages M2 and MCs resting, accompanied by suppressed T cell function. This immunosuppressive state was closely associated with the unique metabolic reprogramming features of the C1 NDUFAB1+ subtype, indicating that targeting metabolic-immune interactions might represent a novel strategy to overcome immune evasion. Based on these findings, we proposed the following therapeutic strategies: First, targeting key metabolic reprogramming pathways (such as oxidative phosphorylation and glycolysis) and developing drugs targeting NDUFAB1 or ELK4 specifically disrupted the metabolic hub of this subtype, effectively inhibiting its energy supply and consequently limiting its proliferation. Second, combined application of mitochondrial inhibitors with glycolysis blockers could help overcome monotherapy resistance caused by “metabolic compensation.” Additionally, combining ELK4 targeting with TME feature-based prognostic models, along with co-administration of immune checkpoint inhibitors and metabolic modulators, could effectively counteract immunosuppression, overcome immune evasion, and enhance anti-tumor immune responses. Finally, patients were classified by detecting characteristic metabolic markers of the C1 NDUFAB1+ subtype (including NDUFAB1 and ELK4 expression levels as well as metabolic enzyme activity profiles), which guided personalized treatment selection. A risk scoring model incorporating metabolic-immune features (integrating parameters such as mitochondrial function indicators and immunosuppressive cell infiltration levels) was explored in combination with other emerging modalities (144–146) to predict patient responses to conventional chemotherapy and immunotherapy, as well as prognosis.
This study provided new targets and strategies for the precision treatment of GC, particularly demonstrating significant advantages in combining targeted metabolic pathways with immunomodulatory therapies. However, there are some limitations. First, the study was constrained by the sample size and geographical origin, which may not fully reflect the broad heterogeneity of GC, especially the differences across various molecular subtypes. Second, the existing ST data and resolution limitations (147) made it difficult to precisely analyze the interactions between individual tumor cells and specific stromal cells. Additionally, although in vitro experiments confirmed the regulatory role of ELK4, further in vivo studies were needed to validate its function in the intact TME. Due to the prolonged experimental timeline and complex sample acquisition procedures, the collection, processing, and staining analysis of samples could not be completed within the short term. Subsequent studies needed to establish prospective cohorts incorporating systematic staining analyses of both NDUFAB1 and ELK4 into the research design. Finally, although we built a risk prediction framework, its clinical application still requires validation through large-scale prospective cohort studies to ensure its clinical utility and reliability. Future studies should further explore the molecular characteristics and clinical translational value of the C1 NDUFAB1+ subtype and its key regulator ELK4 in different GC subtypes, and deepen the understanding of the metabolic-immune interaction mechanisms through multi-omics integrated analysis.
Conclusion
This study successfully deciphered the high heterogeneity of GC through scRNA-seq technology and identified the C1 NDUFAB1+ subtype as a critical cell population. In-depth analysis revealed that this subtype exhibited unique metabolic reprogramming characteristics, including significantly activated oxidative phosphorylation and glycolysis pathways, which potentially drove malignant tumor progression by providing energy and biosynthetic precursors. To further investigate the microenvironmental features of this subtype, we employed ST technology to elucidate its spatial distribution pattern in tumor tissues, discovering its enrichment at the tumor-stroma interface with observed strong spatial interactions. This distinct spatial distribution pattern not only enhanced the understanding of tumor-immune microenvironment interactions but also held potential clinical translational value, such as predicting tumor progression or therapeutic response based on spatial distribution features. In vitro experiments demonstrated that ELK4, as a key molecule of the C1 NDUFAB1+ subtype, significantly promoted the proliferation, migration, and invasion of GC cells. Based on these findings, we constructed a prognostic risk scoring model, which provided an important tool for patient stratification and personalized treatment. In summary, this study systematically elucidated the core mechanisms through which ELK4 mediated tumor progression, metabolic reprogramming, and immune evasion in C1 NDUFAB1+ subtype, offering crucial theoretical foundations for developing precision therapeutic strategies targeting this subtype.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.
Author contributions
YS: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. WN: Conceptualization, Methodology, Software, Writing – original draft, Writing – review & editing. ZX: Data curation, Formal Analysis, Visualization, Writing – original draft, Writing – review & editing. XW: Conceptualization, Data curation, Investigation, Methodology, Software, Writing – review & editing. WL: Conceptualization, Project administration, Validation, Visualization, Writing – review & editing. ZKL: Formal Analysis, Project administration, Validation, Visualization, Writing – review & editing. ZHL: Data curation, Investigation, Methodology, Software, Writing – review & editing. ZDL: Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This study received support from the Shandong Provincial Traditional Chinese Medicine Science and Technology Project (2020M015), the Shandong Provincial Key Research and Development Program (2021CXGC010510), and the Fifth Batch of National Excellent Traditional Chinese Medicine Clinical Talents Training Program (Letter No. 1 of the National Administration of Traditional Chinese Medicine, 2022).
Acknowledgments
Figure 1, as a process summary diagram, was created using Figdraw (Image ID: UYARY27400), and we acknowledged Figdraw for its assistance in this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1591123/full#supplementary-material
Supplementary Figure 1 | Analysis of inferCNV. (A, B) Heatmaps compared the expression distribution of CNVs between ECs and EPCs (upper) and between ECs and four tumor cell subtypes (lower), where red represented amplification and orange represented deletion.
Supplementary Figure 2 | Spatial expression pattern analysis of C1 NDUFAB1+ subtype. (A, B) The ST feature maps displayed the first cell types inferred at selected points on ST 2 slide and the spatial expression pattern of the C1 NDUFAB1+ subtype. (C) The ST feature maps visualized the spatial expression of nCount-Spatial, nFeature-Spatial, G2/M.Score, and S.Score for all cell types in ST 2 slide.
References
1. Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, and Lordick F. Gastric cancer. Lancet. (2020) 396:635–48. doi: 10.1016/S0140-6736(20)31288-5
2. Liu R, Liu J, Cao Q, Chu Y, Chi H, Zhang J, et al. Identification of crucial genes through WGCNA in the progression of gastric cancer. J Cancer. (2024) 15:3284–96. doi: 10.7150/jca.95757
3. Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2024) 74:229–63. doi: 10.3322/caac.21834
4. Qiu MZ, Cai MY, Zhang DS, Wang ZQ, Wang DS, Li YH, et al. Clinicopathological characteristics and prognostic analysis of Lauren classification in gastric adenocarcinoma in China. J Transl Med. (2013) 11:58. doi: 10.1186/1479-5876-11-58
5. Chia NY and Tan P. Molecular classification of gastric cancer. Ann Oncol. (2016) 27:763–69. doi: 10.1093/annonc/mdw040
6. Gong EJ, Lee JY, Bae SE, Park YS, Choi KD, Song HJ, et al. Characteristics of non-cardia gastric cancer with a high serum anti-Helicobacter pylori IgG titer and its association with diffuse-type histology. PloS One. (2018) 13:e0195264. doi: 10.1371/journal.pone.0195264
7. Waldum HL and Fossmark R. Types of gastric carcinomas. Int J Mol Sci. (2018) 19:4109. doi: 10.3390/ijms19124109
8. Li YN, Xie B, Zhang Y, He MH, Xing Y, Mu DM, et al. Advances and key focus areas in gastric cancer immunotherapy: A comprehensive scientometric and clinical trial review (1999-2023). World J Gastroenterol. (2023) 29:5593–617. doi: 10.3748/wjg.v29.i40.5593
9. Xu H, Li T, Shao G, Wang W, He Z, Xu J, et al. Evaluation of neoadjuvant immunotherapy plus chemotherapy in Chinese surgically resectable gastric cancer: a pilot study by meta-analysis. Front Immunol. (2023) 14:1193614. doi: 10.3389/fimmu.2023.1193614
10. Charalampakis N, Economopoulou P, Kotsantis I, Tolia M, Schizas D, Liakakos T, et al. Medical management of gastric cancer: a 2017 update. Cancer Med. (2018) 7:123–33. doi: 10.1002/cam4.1274
11. Wang SF, Chang YL, Fang WL, Li AF, Chen CF, Yeh TS, et al. Growth differentiation factor 15 induces cisplatin resistance through upregulation of xCT expression and glutathione synthesis in gastric cancer. Cancer Sci. (2023) 114:3301–17. doi: 10.1111/cas.15869
12. Wang X, Teng F, Kong L, and Yu J. PD-L1 expression in human cancers and its association with clinical outcomes. Onco Targets Ther. (2016) 9:5023–39. doi: 10.2147/OTT.S105862
13. Wu X, Gu Z, Chen Y, Chen B, Chen W, Weng L, et al. Application of PD-1 blockade in cancer immunotherapy. Comput Struct Biotechnol J. (2019) 17:661–74. doi: 10.1016/j.csbj.2019.03.006
14. Wang BC, Zhang ZJ, Fu C, and Wang C. Efficacy and safety of anti-PD-1/PD-L1 agents vs chemotherapy in patients with gastric or gastroesophageal junction cancer: a systematic review and meta-analysis. Med (Baltimore). (2019) 98:e18054. doi: 10.1097/MD.0000000000018054
15. Cousin S, Guegan JP, Shitara K, Palmieri LJ, Metges JP, Pernot S, et al. Identification of microenvironment features associated with primary resistance to anti-PD-1/PD-L1 + antiangiogenesis in gastric cancer through spatial transcriptomics and plasma proteomics. Mol Cancer. (2024) 23:197. doi: 10.1186/s12943-024-02092-x
16. Li F, Yu J, Pan T, Feng H, Li J, Yu B, et al. BPTF drives gastric cancer resistance to EGFR inhibitor by epigenetically regulating the C-MYC/PLCG1/perk axis. Adv Sci (Weinh). (2023) 10:e2303091. doi: 10.1002/advs.202303091
17. Lei ZN, Teng QX, Tian Q, Chen W, Xie Y, Wu K, et al. Signaling pathways and therapeutic interventions in gastric cancer. Signal Transduct Target Ther. (2022) 7:358. doi: 10.1038/s41392-022-01190-w
18. Yang L, He YT, Dong S, Wei XW, Chen ZH, Zhang B, et al. Single-cell transcriptome analysis revealed a suppressive tumor immune microenvironment in EGFR mutant lung adenocarcinoma. J Immunother Cancer. (2022) 10:e003534. doi: 10.1136/jitc-2021-003534
19. Jiang X, Zhou J, Giobbie-Hurder A, Wargo J, and Hodi FS. The activation of MAPK in melanoma cells resistant to BRAF inhibition promotes PD-L1 expression that is reversible by MEK and PI3K inhibition. Clin Cancer Res. (2013) 19:598–609. doi: 10.1158/1078-0432.CCR-12-2731
20. Qin X, Liu C, Zhou Y, and Wang G. Cisplatin induces programmed death-1-ligand 1(PD-L1) over-expression in hepatoma H22 cells via Erk/MAPK signaling pathway. Cell Mol Biol (Noisy-Le-Grand). (2010) 56 Suppl:OL1366–72.
21. Zhang X, Zhang P, Cong A, Feng Y, Chi H, Xia Z, et al. Unraveling molecular networks in thymic epithelial tumors: deciphering the unique signatures. Front Immunol. (2023) 14:1264325. doi: 10.3389/fimmu.2023.1264325
22. Liu Y, Li C, Lu Y, Liu C, and Yang W. Tumor microenvironment-mediated immune tolerance in development and treatment of gastric cancer. Front Immunol. (2022) 13:1016817. doi: 10.3389/fimmu.2022.1016817
23. He D, Chen M, Chang L, Gu J, Liu F, Gao X, et al. De novo pyrimidine synthesis fuels glycolysis and confers chemoresistance in gastric cancer. Cancer Lett. (2022) 549:215837. doi: 10.1016/j.canlet.2022.215837
24. Ilkhani K, Bastami M, Delgir S, Safi A, Talebian S, and Alivand MR. The engaged role of tumor microenvironment in cancer metabolism: focusing on cancer-associated fibroblast and exosome mediators. Anticancer Agents Med Chem. (2021) 21:254–66. doi: 10.2174/1871520620666200910123428
25. Bin YL, Hu HS, Tian F, Wen ZH, Yang MF, Wu BH, et al. Metabolic reprogramming in gastric cancer: trojan horse effect. Front Oncol. (2021) 11:745209. doi: 10.3389/fonc.2021.745209
26. Li Y, Hu X, Lin R, Zhou G, Zhao L, Zhao D, et al. Single-cell landscape reveals active cell subtypes and their interaction in the tumor microenvironment of gastric cancer. Theranostics. (2022) 12:3818–33. doi: 10.7150/thno.71833
27. Li X, Sun Z, Peng G, Xiao Y, Guo J, Wu B, et al. Single-cell RNA sequencing reveals a pro-invasive cancer-associated fibroblast subgroup associated with poor clinical outcomes in patients with gastric cancer. Theranostics. (2022) 12:620–38. doi: 10.7150/thno.60540
28. Sun K, Xu R, Ma F, Yang N, Li Y, Sun X, et al. scRNA-seq of gastric tumor shows complex intercellular interaction with an alternative T cell exhaustion trajectory. Nat Commun. (2022) 13:4943. doi: 10.1038/s41467-022-32627-z
29. Tang J, Wei W, Xu Y, Chen K, Miao Y, Fan W, et al. CXC chemokine receptor 4 - mediated immune modulation and tumor microenvironment heterogeneity in gastric cancer: Utilizing multi-omics approaches to identify potential therapeutic targets. Biofactors. (2025) 51:e2130. doi: 10.1002/biof.2130
30. Ma F, Wang L, Chi H, Li X, Xu Y, Chen K, et al. Exploring the therapeutic potential of MIR-140-3p in osteoarthritis: targeting CILP and ferroptosis for novel treatment strategies. Cell Prolif. (2025) e70018. doi: 10.1111/cpr.70018
31. Xu Q, Liu C, Wang H, Li S, Yan H, Liu Z, et al. Deciphering the impact of aggregated autophagy-related genes TUBA1B and HSP90AA1 on colorectal cancer evolution: a single-cell sequencing study of the tumor microenvironment. Discov Oncol. (2024) 15:431. doi: 10.1007/s12672-024-01322-4
32. Zhai X, Zhang H, Xia Z, Liu M, Du G, Jiang Z, et al. Oxytocin alleviates liver fibrosis via hepatic macrophages. JHEP Rep. (2024) 6:101032. doi: 10.1016/j.jhepr.2024.101032
33. Day BW, Stringer BW, Spanevello MD, Charmsaz S, Jamieson PR, Ensbey KS, et al. ELK4 neutralization sensitizes glioblastoma to apoptosis through downregulation of the anti-apoptotic protein Mcl-1. Neuro Oncol. (2011) 13:1202–12. doi: 10.1093/neuonc/nor119
34. Zhu Z, Guo Y, Liu Y, Ding R, Huang Z, Yu W, et al. ELK4 promotes colorectal cancer progression by activating the neoangiogenic factor LRG1 in a noncanonical SP1/3-dependent manner. Adv Sci (Weinh). (2023) 10:e2303378. doi: 10.1002/advs.202303378
35. Lin Z, Wang F, Yin R, Li S, Bai Y, Zhang B, et al. Single-cell RNA sequencing and immune microenvironment analysis reveal PLOD2-driven Malignant transformation in cervical cancer. Front Immunol. (2024) 15:1522655. doi: 10.3389/fimmu.2024.1522655
36. Li X, Lin Z, Zhao F, Huang T, Fan W, Cen L, et al. Unveiling the cellular landscape: insights from single-cell RNA sequencing in multiple myeloma. Front Immunol. (2024) 15:1458638. doi: 10.3389/fimmu.2024.1458638
37. Sun L, Shao W, Lin Z, Lin J, Zhao F, and Yu J. Single-cell RNA sequencing explored potential therapeutic targets by revealing the tumor microenvironment of neuroblastoma and its expression in cell death. Discov Oncol. (2024) 15:409. doi: 10.1007/s12672-024-01286-5
38. Zhao F, Hong J, Zhou G, Huang T, Lin Z, Zhang Y, et al. Elucidating the role of tumor-associated ALOX5+ mast cells with transformative function in cervical cancer progression via single-cell RNA sequencing. Front Immunol. (2024) 15:1434450. doi: 10.3389/fimmu.2024.1434450
39. Shao W, Lin Z, Xiahou Z, Zhao F, Xu J, Liu X, et al. Single-cell RNA sequencing reveals that MYBL2 in Malignant epithelial cells is involved in the development and progression of ovarian cancer. Front Immunol. (2024) 15:1438198. doi: 10.3389/fimmu.2024.1438198
40. Hou M, Zhao Z, Li S, Zhang Z, Li X, Zhang Y, et al. Single-cell analysis unveils cell subtypes of acral melanoma cells at the early and late differentiation stages. J Cancer. (2025) 16:898–916. doi: 10.7150/jca.102045
41. Li H, Bian Y, Xiahou Z, Zhao Z, Zhao F, and Zhang Q. The cellular signaling crosstalk between memory B cells and tumor cells in nasopharyngeal carcinoma cannot be overlooked: Their involvement in tumor progression and treatment strategy is significant. J Cancer. (2025) 16:288–314. doi: 10.7150/jca.101420
42. Lin L, Zou J, Pei S, Huang W, Zhang Y, Zhao Z, et al. Germinal center B-cell subgroups in the tumor microenvironment cannot be overlooked: Their involvement in prognosis, immunotherapy response, and treatment resistance in head and neck squamous carcinoma. Heliyon. (2024) 10:e37726. doi: 10.1016/j.heliyon.2024.e37726
43. Zhang Y, Zhao Z, Huang W, Kim BS, Lin L, Li X, et al. Pan-cancer single-cell analysis revealing the heterogeneity of cancer-associated fibroblasts in skin tumors. Curr Gene Ther. (2024). doi: 10.2174/0115665232331353240911080642
44. Jin W, Zhang Y, Zhao Z, and Gao M. Developing targeted therapies for neuroblastoma by dissecting the effects of metabolic reprogramming on tumor microenvironments and progression. Theranostics. (2024) 14:3439–69. doi: 10.7150/thno.93962
45. Nie W, Zhao Z, Liu Y, Wang Y, Zhang J, Hu Y, et al. Integrative single-cell analysis of cardiomyopathy identifies differences in cell stemness and transcriptional regulatory networks among fibroblast subpopulations. Cardiol Res Pract. (2024) 2024:3131633. doi: 10.1155/2024/3131633
46. Huang W, Kim BS, Zhang Y, Lin L, Chai G, and Zhao Z. Regulatory T cells subgroups in the tumor microenvironment cannot be overlooked: Their involvement in prognosis and treatment strategy in melanoma. Environ Toxicol. (2024) 39:4512–30. doi: 10.1002/tox.24247
47. Ge Q, Zhao Z, Li X, Yang F, Zhang M, Hao Z, et al. Deciphering the suppressive immune microenvironment of prostate cancer based on CD4+ regulatory T cells: Implications for prognosis and therapy prediction. Clin Transl Med. (2024) 14:e1552. doi: 10.1002/ctm2.1552
48. Ding Y, Zhao Z, Cai H, Zhou Y, Chen H, Bai Y, et al. Single-cell sequencing analysis related to sphingolipid metabolism guides immunotherapy and prognosis of skin cutaneous melanoma. Front Immunol. (2023) 14:1304466. doi: 10.3389/fimmu.2023.1304466
49. Zhao F, Jiang X, Li Y, Huang T, Xiahou Z, Nie W, et al. Characterizing tumor biology and immune microenvironment in high-grade serous ovarian cancer via single-cell RNA sequencing: insights for targeted and personalized immunotherapy strategies. Front Immunol. (2024) 15:1500153. doi: 10.3389/fimmu.2024.1500153
50. An Y, Zhao F, Jia H, Meng S, Zhang Z, Li S, et al. Inhibition of programmed cell death by melanoma cell subpopulations reveals mechanisms of melanoma metastasis and potential therapeutic targets. Discov Oncol. (2025) 16:62. doi: 10.1007/s12672-025-01789-9
51. Wang J, Zhao F, Zhang Q, Sun Z, Xiahou Z, Wang C, et al. Unveiling the NEFH+ Malignant cell subtype: Insights from single-cell RNA sequencing in prostate cancer progression and tumor microenvironment interactions. Front Immunol. (2024) 15:1517679. doi: 10.3389/fimmu.2024.1517679
52. Ni G, Sun Y, Jia H, Xiahou Z, Li Y, Zhao F, et al. MAZ-mediated tumor progression and immune evasion in hormone receptor-positive breast cancer: Targeting tumor microenvironment and PCLAF+ subtype-specific therapy. Transl Oncol. (2025) 52:102280. doi: 10.1016/j.tranon.2025.102280
53. Luo S, Wang L, Xiao Y, Cao C, Liu Q, and Zhou Y. Single-cell RNA-sequencing integration analysis revealed immune cell heterogeneity in five human autoimmune diseases. Bio Integration. (2023) 4:145–59. doi: 10.15212/bioi-2023-0012
54. Zhou W, Lin Z, and Tan W. Deciphering the molecular landscape: integrating single-cell transcriptomics to unravel myofibroblast dynamics and therapeutic targets in clear cell renal cell carcinomas. Front Immunol. (2024) 15:1374931. doi: 10.3389/fimmu.2024.1374931
55. Liu P, Xing N, Xiahou Z, Yan J, Lin Z, and Zhang J. Unraveling the intricacies of glioblastoma progression and recurrence: insights into the role of NFYB and oxidative phosphorylation at the single-cell level. Front Immunol. (2024) 15:1368685. doi: 10.3389/fimmu.2024.1368685
56. Xing J, Cai H, Lin Z, Zhao L, Xu H, Song Y, et al. Examining the function of macrophage oxidative stress response and immune system in glioblastoma multiforme through analysis of single-cell transcriptomics. Front Immunol. (2023) 14:1288137. doi: 10.3389/fimmu.2023.1288137
57. Zhao Z, Ding Y, Tran LJ, Chai G, and Lin L. Innovative breakthroughs facilitated by single-cell multi-omics: manipulating natural killer cell functionality correlates with a novel subcategory of melanoma cells. Front Immunol. (2023) 14:1196892. doi: 10.3389/fimmu.2023.1196892
58. Talevich E, Shain AH, Botton T, and Bastian BC. CNVkit: genome-wide copy number detection and visualization from targeted DNA sequencing. PloS Comput Biol. (2016) 12:e1004873. doi: 10.1371/journal.pcbi.1004873
59. Zhao ZJ, Wei DP, Zheng RZ, Peng T, Xiao X, and Li FS. The gene coexpression analysis identifies functional modules dynamically changed after traumatic brain injury. Comput Math Methods Med. (2021) 2021:5511598. doi: 10.1155/2021/5511598
60. Lin Z, Fan W, Yu X, Liu J, and Liu P. Research into the mechanism of intervention of SanQi in endometriosis based on network pharmacology and molecular docking technology. Med (Baltimore). (2022) 101:e30021. doi: 10.1097/MD.0000000000030021
61. He Y, Luo Z, Nie X, Du Y, Sun R, Sun J, et al. An injectable multi-functional composite bioactive hydrogel for bone regeneration via immunoregulatory and osteogenesis effects. Adv Compos Hybrid Mater. (2025) 8:128. doi: 10.1007/s42114-025-01213-4
62. Zhao ZJ, Zheng RZ, Wang XJ, Li TQ, Dong XH, Zhao CY, et al. Integrating lipidomics and transcriptomics reveals the crosstalk between oxidative stress and neuroinflammation in central nervous system demyelination. Front Aging Neurosci. (2022) 14:870957. doi: 10.3389/fnagi.2022.870957
63. Zhao Z, Li T, Dong X, Wang X, Zhang Z, Zhao C, et al. Untargeted metabolomic profiling of cuprizone-induced demyelination in mouse corpus callosum by UPLC-orbitrap/MS reveals potential metabolic biomarkers of CNS demyelination disorders. Oxid Med Cell Longev. (2021) 2021:7093844. doi: 10.1155/2021/7093844
64. Wang Y, Zhao ZJ, Kang XR, Bian T, Shen ZM, Jiang Y, et al. lncRNA DLEU2 acts as a miR-181a sponge to regulate SEPP1 and inhibit skeletal muscle differentiation and regeneration. Aging (Albany NY). (2020) 12:24033–56. doi: 10.18632/aging.104095
65. Li XY, Zhao ZJ, Wang JB, Shao YH, Hui-Liu, You JX, et al. m7G methylation-related genes as biomarkers for predicting overall survival outcomes for hepatocellular carcinoma. Front Bioeng Biotechnol. (2022) 10:849756. doi: 10.3389/fbioe.2022.849756
66. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
67. Zhang P, Wang L, Liu H, Lin S, and Guo D. Unveiling the crucial role of glycosylation modification in lung adenocarcinoma metastasis through artificial neural network-based spatial multi-omics single-cell analysis and Mendelian randomization. BMC Cancer. (2025) 25:249. doi: 10.1186/s12885-025-13650-x
68. Zhang P, Wen B, Gong J, Liu Z, Zhang M, Zhou G, et al. Clinical prognostication and immunotherapy response prediction in esophageal squamous cell carcinoma using the DNA damage repair-associated signature. Environ Toxicol. (2024) 39:2803–16. doi: 10.1002/tox.24155
69. Lin Z, Fan W, Sui X, Wang J, and Zhao J. Necroptosis-related lncRNA signatures for prognostic prediction in uterine corpora endometrial cancer. Reprod Sci. (2023) 30:576–89. doi: 10.1007/s43032-022-01023-9
70. Lin Z, Sui X, Jiao W, Wang Y, and Zhao J. Exploring the mechanism and experimental verification of puerarin in the treatment of endometrial carcinoma based on network pharmacology and bioinformatics analysis. BMC Complement Med Ther. (2022) 22:150. doi: 10.1186/s12906-022-03623-z
71. Zou J, Lin Z, Jiao W, Chen J, Lin L, Zhang F, et al. A multi-omics-based investigation of the prognostic and immunological impact of necroptosis-related mRNA in patients with cervical squamous carcinoma and adenocarcinoma. Sci Rep. (2022) 12:16773. doi: 10.1038/s41598-022-20566-0
72. Zhang P, Wang D, Zhou G, Jiang S, Zhang G, Zhang L, et al. Novel post-translational modification learning signature reveals B4GALT2 as an immune exclusion regulator in lung adenocarcinoma. J Immunother Cancer. (2025) 13:e010787. doi: 10.1136/jitc-2024-010787
73. Du H, Li S, Lu J, Tang L, Jiang X, He X, et al. Single-cell RNA-seq and bulk-seq identify RAB17 as a potential regulator of angiogenesis by human dermal microvascular endothelial cells in diabetic foot ulcers. Burns Trauma. (2023) 11:tkad020. doi: 10.1093/burnst/tkad020
74. Fu Y, Zhang J, Liu Q, Yang L, Wu Q, Yang X, et al. Unveiling the role of ABI3 and hub senescence-related genes in macrophage senescence for atherosclerotic plaque progression. Inflammation Res. (2024) 73:65–82. doi: 10.1007/s00011-023-01817-w
75. Lin Z, Zou J, Sui X, Yao S, Lin L, Wang J, et al. Necroptosis-related lncRNA signature predicts prognosis and immune response for cervical squamous cell carcinoma and endocervical adenocarcinomas. Sci Rep. (2022) 12:16285. doi: 10.1038/s41598-022-20858-5
76. Zhao J, Zou J, Jiao W, Lin L, Wang J, and Lin Z. Construction of N-7 methylguanine-related mRNA prognostic model in uterine corpus endometrial carcinoma based on multi-omics data and immune-related analysis. Sci Rep. (2022) 12:18813. doi: 10.1038/s41598-022-22879-6
77. Zhao J, Jiao W, Sui X, Zou J, Wang J, and Lin Z. Construction of a prognostic model of luteolin for endometrial carcinoma. Am J Transl Res. (2023) 15:2122–39.
78. Li J, Liang L, Xiu L, Zeng J, Zhu Y, An J, et al. Establishment of a molecular risk model for the prognosis of cervical cancer based on microRNA expression. Ann Transl Med. (2022) 10:125. doi: 10.21037/atm-21-6451
79. Zhu Y, Liang L, Zhao Y, Li J, Zeng J, Yuan Y, et al. CircNUP50 is a novel therapeutic target that promotes cisplatin resistance in ovarian cancer by modulating p53 ubiquitination. J Nanobiotechnology. (2024) 22:35. doi: 10.1186/s12951-024-02295-w
80. Tu H, Hu Q, Ma Y, Huang J, Luo H, Jiang L, et al. Deciphering the tumour microenvironment of clear cell renal cell carcinoma: Prognostic insights from programmed death genes using machine learning. J Cell Mol Med. (2024) 28:e18524. doi: 10.1111/jcmm.18524
81. Zhao ZJ, Chen D, Zhou LY, Sun ZL, Wang BC, and Feng DF. Prognostic value of different computed tomography scoring systems in patients with severe traumatic brain injury undergoing decompressive craniectomy. J Comput Assist Tomogr. (2022) 46:800–07. doi: 10.1097/RCT.0000000000001343
82. Zheng RZ, Zhao ZJ, Yang XT, Jiang SW, Li YD, Li WJ, et al. Initial CT-based radiomics nomogram for predicting in-hospital mortality in patients with traumatic brain injury: a multicenter development and validation study. Neurol Sci. (2022) 43:4363–72. doi: 10.1007/s10072-022-05954-8
83. Zheng R, Zhuang Z, Zhao C, Zhao Z, Yang X, Zhou Y, et al. Chinese admission warning strategy for predicting the hospital discharge outcome in patients with traumatic brain injury. J Clin Med. (2022) 11:974. doi: 10.3390/jcm11040974
84. Luo Z, He Z, Qin H, Chen Y, Qi B, Lin J, et al. Exercise-induced IL-15 acted as a positive prognostic implication and tumor-suppressed role in pan-cancer. Front Pharmacol. (2022) 13:1053137. doi: 10.3389/fphar.2022.1053137
85. Wang G, Li Y, Pan R, Yin X, Jia C, She Y, et al. XRCC1: a potential prognostic and immunological biomarker in LGG based on systematic pan-cancer analysis. Aging (Albany NY). (2024) 16:872–910. doi: 10.18632/aging.205426
86. Yan Z, Liu Y, and Yuan Y. The plasticity of epithelial cells and its potential in the induced differentiation of gastric cancer. Cell Death Discov. (2024) 10:512. doi: 10.1038/s41420-024-02275-x
87. Fan T, Jiang L, Zhou X, Chi H, and Zeng X. Deciphering the dual roles of PHD finger proteins from oncogenic drivers to tumor suppressors. Front Cell Dev Biol. (2024) 12:1403396. doi: 10.3389/fcell.2024.1403396
88. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, and Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2018) 68:394–424. doi: 10.3322/caac.21492
89. Allemani C, Matsuda T, Di Carlo V, Harewood R, Matz M, Niksic M, et al. Global surveillance of trends in cancer survival 2000-14 (CONCORD-3): analysis of individual records for 37 513–025 patients diagnosed with one of 18 cancers from 322 population-based registries in 71 countries. Lancet. (2018) 391:1023–75. doi: 10.1016/S0140-6736(17)33326-3
90. Li GC, Jia XC, Zhao QC, Zhang HW, Yang P, Xu LL, et al. The expression of epidermal growth factor receptor 1 and human epidermal growth factor receptor 2 based on tumor location affect survival in gastric cancer. Med (Baltimore). (2020) 99:e20460. doi: 10.1097/MD.0000000000020460
91. Feng J, Zhang P, Wang D, Li Y, and Tan J. New strategies for lung cancer diagnosis and treatment: applications and advances in nanotechnology. biomark Res. (2024) 12:136. doi: 10.1186/s40364-024-00686-7
92. Zhang F, Wang H, Yu J, Yao X, Yang S, Li W, et al. LncRNA CRNDE attenuates chemoresistance in gastric cancer via SRSF6-regulated alternative splicing of PICALM. Mol Cancer. (2021) 20:6. doi: 10.1186/s12943-020-01299-y
93. Li J, Sun J, Zeng Z, Liu Z, Ma M, Zheng Z, et al. Tumour-associated macrophages in gastric cancer: From function and mechanism to application. Clin Transl Med. (2023) 13:e1386. doi: 10.1002/ctm2.1386
94. Shi T, Zhang Y, Wang Y, Song X, Wang H, Zhou X, et al. DKK1 promotes tumor immune evasion and impedes anti-PD-1 treatment by inducing immunosuppressive macrophages in gastric cancer. Cancer Immunol Res. (2022) 10:1506–24. doi: 10.1158/2326-6066.CIR-22-0218
95. Sorin M, Karimi E, Rezanejad M, Yu MW, Desharnais L, McDowell S, et al. Single-cell spatial landscape of immunotherapy response reveals mechanisms of CXCL13 enhanced antitumor immunity. J Immunother Cancer. (2023) 11:e005545. doi: 10.1136/jitc-2022-005545
96. Yuan Q, Deng D, Pan C, Ren J, Wei T, Wu Z, et al. Integration of transcriptomics, proteomics, and metabolomics data to reveal HER2-associated metabolic heterogeneity in gastric cancer with response to immunotherapy and neoadjuvant chemotherapy. Front Immunol. (2022) 13:951137. doi: 10.3389/fimmu.2022.951137
97. Xie J, Yang A, Liu Q, Deng X, Lv G, Ou X, et al. Single-cell RNA sequencing elucidated the landscape of breast cancer brain metastases and identified ILF2 as a potential therapeutic target. Cell Prolif. (2024) 57:e13697. doi: 10.1111/cpr.13697
98. Yuan Z, Li B, Liao W, Kang D, Deng X, Tang H, et al. Comprehensive pan-cancer analysis of YBX family reveals YBX2 as a potential biomarker in liver cancer. Front Immunol. (2024) 15:1382520. doi: 10.3389/fimmu.2024.1382520
99. Liu P, Deng X, Zhou H, Xie J, Kong Y, Zou Y, et al. Multi-omics analyses unravel DNA damage repair-related clusters in breast cancer with experimental validation. Front Immunol. (2023) 14:1297180. doi: 10.3389/fimmu.2023.1297180
100. Zhang X, Song W, Gao Y, Zhang Y, Zhao Y, Hao S, et al. The role of tumor metabolic reprogramming in tumor immunity. Int J Mol Sci. (2023) 24:17422. doi: 10.3390/ijms242417422
101. Liu Y, Zhao Y, Song H, Li Y, Liu Z, Ye Z, et al. Metabolic reprogramming in tumor immune microenvironment: Impact on immune cell function and therapeutic implications. Cancer Lett. (2024) 597:217076. doi: 10.1016/j.canlet.2024.217076
102. Goodall GJ and Wickramasinghe VO. RNA in cancer. Nat Rev Cancer. (2021) 21:22–36. doi: 10.1038/s41568-020-00306-0
103. Durgeau A, Virk Y, Corgnac S, and Mami-Chouaib F. Recent advances in targeting CD8 T-cell immunity for more effective cancer immunotherapy. Front Immunol. (2018) 9:14. doi: 10.3389/fimmu.2018.00014
104. Runger TM, Klein CE, Becker JC, and Brocker EB. The role of genetic instability, adhesion, cell motility, and immune escape mechanisms in melanoma progression. Curr Opin Oncol. (1994) 6:188–96. doi: 10.1097/00001622-199403000-00012
105. Oakes SA. Endoplasmic reticulum stress signaling in cancer cells. Am J Pathol. (2020) 190:934–46. doi: 10.1016/j.ajpath.2020.01.010
106. Sancho P, Barneda D, and Heeschen C. Hallmarks of cancer stem cell metabolism. Br J Cancer. (2016) 114:1305–12. doi: 10.1038/bjc.2016.152
107. Zhang S, Zhao H, Chen Y, and Zhang Y. GNL3 regulates SIRT1 transcription and promotes hepatocellular carcinoma stem cell-like features and metastasis. J Oncol. (2022) 2022:1555670. doi: 10.1155/2022/1555670
108. Vlashi E and Pajonk F. Cancer stem cells, cancer cell plasticity and radiation therapy. Semin Cancer Biol. (2015) 31:28–35. doi: 10.1016/j.semcancer.2014.07.001
109. He S and Tang S. WNT/beta-catenin signaling in the development of liver cancers. BioMed Pharmacother. (2020) 132:110851. doi: 10.1016/j.biopha.2020.110851
110. Baluapuri A, Wolf E, and Eilers M. Target gene-independent functions of MYC oncoproteins. Nat Rev Mol Cell Biol. (2020) 21:255–67. doi: 10.1038/s41580-020-0215-2
111. Krenz B, Lee J, Kannan T, and Eilers M. Immune evasion: An imperative and consequence of MYC deregulation. Mol Oncol. (2024) 18:2338–55. doi: 10.1002/1878-0261.13695
112. Stacy AE, Jansson PJ, and Richardson DR. Molecular pharmacology of ABCG2 and its role in chemoresistance. Mol Pharmacol. (2013) 84:655–69. doi: 10.1124/mol.113.088609
113. Salagacka-Kubiak A, Zawada D, Saed L, Kordek R, Jelen A, and Balcerczak E. ABCG2 gene and ABCG2 protein expression in colorectal cancer-in silico and wet analysis. Int J Mol Sci. (2023) 24:10539. doi: 10.3390/ijms241310539
114. Katoh M. Canonical and non-canonical WNT signaling in cancer stem cells and their niches: Cellular heterogeneity, omics reprogramming, targeted therapy and tumor plasticity (Review). Int J Oncol. (2017) 51:1357–69. doi: 10.3892/ijo.2017.4129
115. Qi F, Li J, Qi Z, Zhang J, Zhou B, Yang B, et al. Comprehensive metabolic profiling and genome-wide analysis reveal therapeutic modalities for hepatocellular carcinoma. Res (Wash D C). (2023) 6:36. doi: 10.34133/research.0036
116. Xu X, Zhang M, Xu F, and Jiang S. Wnt signaling in breast cancer: biological mechanisms, challenges and opportunities. Mol Cancer. (2020) 19:165. doi: 10.1186/s12943-020-01276-5
117. Arora P, Cuevas BD, Russo A, Johnson GL, and Trejo J. Persistent transactivation of EGFR and ErbB2/HER2 by protease-activated receptor-1 promotes breast carcinoma cell invasion. Oncogene. (2008) 27:4434–45. doi: 10.1038/onc.2008.84
118. Grisaru-Granovsky S, Salah Z, Maoz M, Pruss D, Beller U, and Bar-Shavit R. Differential expression of protease activated receptor 1 (Par1) and pY397FAK in benign and Malignant human ovarian tissue samples. Int J Cancer. (2005) 113:372–78. doi: 10.1002/ijc.20607
119. Greenberg DL, Mize GJ, and Takayama TK. Protease-activated receptor mediated RhoA signaling and cytoskeletal reorganization in LNCaP cells. Biochemistry. (2003) 42:702–09. doi: 10.1021/bi027100x
120. Sedda S, Marafini I, Caruso R, Pallone F, and Monteleone G. Proteinase activated-receptors-associated signaling in the control of gastric cancer. World J Gastroenterol. (2014) 20:11977–84. doi: 10.3748/wjg.v20.i34.11977
121. Fujimoto D, Hirono Y, Goi T, Katayama K, Hirose K, and Yamaguchi A. Expression of protease activated receptor-2 (PAR-2) in gastric cancer. J Surg Oncol. (2006) 93:139–44. doi: 10.1002/jso.20420
122. Pawar NR, Buzza MS, and Antalis TM. Membrane-anchored serine proteases and protease-activated receptor-2-mediated signaling: co-conspirators in cancer progression. Cancer Res. (2019) 79:301–10. doi: 10.1158/0008-5472.CAN-18-1745
123. Arora P, Ricks TK, and Trejo J. Protease-activated receptor signalling, endocytic sorting and dysregulation in cancer. J Cell Sci. (2007) 120:921–28. doi: 10.1242/jcs.03409
124. Rad SM, Langroudi L, Kouhkan F, Yazdani L, Koupaee AN, Asgharpour S, et al. Transcription factor decoy: a pre-transcriptional approach for gene downregulation purpose in cancer. Tumour Biol. (2015) 36:4871–81. doi: 10.1007/s13277-015-3344-z
125. Lourenco C, Resetca D, Redel C, Lin P, MacDonald AS, Ciaccio R, et al. MYC protein interactors in gene transcription and cancer. Nat Rev Cancer. (2021) 21:579–91. doi: 10.1038/s41568-021-00367-9
126. Mao X, Ji T, Liu A, and Weng Y. ELK4-mediated lncRNA SNHG22 promotes gastric cancer progression through interacting with EZH2 and regulating miR-200c-3p/Notch1 axis. Cell Death Dis. (2021) 12:957. doi: 10.1038/s41419-021-04228-z
127. Zheng L, Xu H, Di Y, Chen L, Liu J, Kang L, et al. ELK4 promotes the development of gastric cancer by inducing M2 polarization of macrophages through regulation of the KDM5A-PJA2-KSR1 axis. J Transl Med. (2021) 19:342. doi: 10.1186/s12967-021-02915-1
128. Qin H, Luo Z, Sun Y, He Z, Qi B, Chen Y, et al. Low-intensity pulsed ultrasound promotes skeletal muscle regeneration via modulating the inflammatory immune microenvironment. Int J Biol Sci. (2023) 19:1123–45. doi: 10.7150/ijbs.79685
129. Deng Y, Shi M, Yi L, Naveed KM, Xia Z, and Li X. Eliminating a barrier: Aiming at VISTA, reversing MDSC-mediated T cell suppression in the tumor microenvironment. Heliyon. (2024) 10:e37060. doi: 10.1016/j.heliyon.2024.e37060
130. Xie H, Xi X, Lei T, Liu H, and Xia Z. CD8(+) T cell exhaustion in the tumor microenvironment of breast cancer. Front Immunol. (2024) 15:1507283. doi: 10.3389/fimmu.2024.1507283
131. Feng X, Luo Z, Zhang W, Wan R, Chen Y, Li F, et al. Zn-DHM nanozymes enhance muscle regeneration through ROS scavenging and macrophage polarization in volumetric muscle loss revealed by single-cell profiling. Adv Funct Mater. (2025), 2506476. doi: 10.1002/adfm.202506476
132. Wang JX, Choi S, Niu X, Kang N, Xue H, Killam J, et al. Lactic acid and an acidic tumor microenvironment suppress anticancer immunity. Int J Mol Sci. (2020) 21:8363. doi: 10.3390/ijms21218363
133. Chen H, Zuo H, Huang J, Liu J, Jiang L, Jiang C, et al. Unravelling infiltrating T-cell heterogeneity in kidney renal clear cell carcinoma: Integrative single-cell and spatial transcriptomic profiling. J Cell Mol Med. (2024) 28:e18403. doi: 10.1111/jcmm.18403
134. He S, Su L, Hu H, Liu H, Xiong J, Gong X, et al. Immunoregulatory functions and therapeutic potential of natural killer cell-derived extracellular vesicles in chronic diseases. Front Immunol. (2023) 14:1328094. doi: 10.3389/fimmu.2023.1328094
135. Xu K, Wu Y, Chi H, Li Y, She Y, Yin X, et al. SLC22A8: An indicator for tumor immune microenvironment and prognosis of ccRCC from a comprehensive analysis of bioinformatics. Med (Baltimore). (2022) 101:e30270. doi: 10.1097/MD.0000000000030270
136. Chen F, Shen L, Wang Y, Chen Y, Pan X, Liang H, et al. Signatures of immune cell infiltration for predicting immune escape and immunotherapy in cervical cancer. Aging (Albany NY). (2023) 15:1685–98. doi: 10.18632/aging.204583
137. Mardis ER. Neoantigens and genome instability: impact on immunogenomic phenotypes and immunotherapy response. Genome Med. (2019) 11:71. doi: 10.1186/s13073-019-0684-0
138. Situ Y, Zhang J, Liao W, Liang Q, Lu L, Xu Q, et al. SHMT as a potential therapeutic target for renal cell carcinoma. Front Biosci (Landmark Ed). (2023) 28:196. doi: 10.31083/j.fbl2809196
139. Gong C, Yang M, Long H, Liu X, Xu Q, Qiao L, et al. IL-6-driven autocrine lactate promotes immune escape of uveal melanoma. Invest Ophthalmol Vis Sci. (2024) 65:37. doi: 10.1167/iovs.65.3.37
140. Prabhakaran HS, Hu D, He W, Luo G, and Liou YC. Mitochondrial dysfunction and mitophagy: crucial players in burn trauma and wound healing. Burns Trauma. (2023) 11:tkad029. doi: 10.1093/burnst/tkad029
141. Zhang H, Zhai X, Liu Y, Xia Z, Xia T, Du G, et al. NOP2-mediated m5C Modification of c-Myc in an EIF3A-Dependent Manner to Reprogram Glucose Metabolism and Promote Hepatocellular Carcinoma Progression. Res (Wash D C). (2023) 6:184. doi: 10.34133/research.0184
142. Niu B, Liao K, Zhou Y, Wen T, Quan G, Pan X, et al. Application of glutathione depletion in cancer therapy: Enhanced ROS-based therapy, ferroptosis, and chemotherapy. Biomaterials. (2021) 277:121110. doi: 10.1016/j.biomaterials.2021.121110
143. Yang Z, Chen Y, Miao Y, Yan H, Chen K, Xu Y, et al. Elucidating stearoyl metabolism and NCOA4-mediated ferroptosis in gastric cancer liver metastasis through multi-omics single-cell integrative mendelian analysis: advancing personalized immunotherapy strategies. Discov Oncol. (2025) 16:46. doi: 10.1007/s12672-025-01769-z
144. Tang S, Gokbag B, Fan K, Shao S, Huo Y, Wu X, et al. Synthetic lethal gene pairs: Experimental approaches and predictive models. Front Genet. (2022) 13:961611. doi: 10.3389/fgene.2022.961611
145. Song CY, Zhu JY, Huang W, and Lu YQ. Development and validation of a predictive model for the assessment of potassium-lowering treatment among hyperkalemia patients. World J Emerg Med. (2023) 14:198–203. doi: 10.5847/wjem.j.1920-8642.2023.048
146. Deng X, Yan Y, Zhan Z, Xie J, Tang H, Zou Y, et al. A glimpse into the future: Integrating artificial intelligence for precision HER2-positive breast cancer management. IMetaOmics. (2024) 1:e19. doi: 10.1002/imo2.19
Keywords: gastric cancer, scRNA-seq, spatial transcriptomics, tumor microenvironment, metabolic reprogramming, immune evasion
Citation: Sun Y, Nie W, Xiahou Z, Wang X, Liu W, Liu Z, Lin Z and Liu Z (2025) Integrative single-cell and spatial transcriptomics uncover ELK4-mediated mechanisms in NDUFAB1+ tumor cells driving gastric cancer progression, metabolic reprogramming, and immune evasion. Front. Immunol. 16:1591123. doi: 10.3389/fimmu.2025.1591123
Received: 10 March 2025; Accepted: 09 June 2025;
Published: 04 July 2025.
Edited by:
Lilong Zhang, Renmin Hospital of Wuhan University, ChinaReviewed by:
Xinpei Deng, Sun Yat-Sen University Cancer Center (SYSUCC), ChinaXiang-Xu Wang, Air Force Medical University, China
Yuquan Chen, Monash University, Australia
Kunyu Wang, Chinese Academy of Medical Sciences and Peking Union Medical College, China
Hui Dong, Shanxi Medical University, China
Kui Wang, Shandong University, China
Copyright © 2025 Sun, Nie, Xiahou, Wang, Liu, Liu, Lin and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhaidong Liu, emRsaXU3NEAxNjMuY29t; Zhiheng Lin, bGluemhpaGVuZ0BzaHV0Y20uZWR1LmNu
†These authors have contributed equally to this work and share first authorship