Weighted Gene Coexpression Network Analysis Identifies Neutrophil-Related Molecular Subtypes and Their Clinical Significance in Gastric Cancer

Introduction

Gastric cancer (GC), which originates from epithelial cells of the gastric mucosa, is a common malignant tumor of the digestive tract.¹ It is the fifth most common cancer worldwide and the third leading cause of cancer death.² During the division and proliferation of GC cells, cell subtypes with different biological functions and associated with different prognoses can be generated due to differences in the tumor microenvironment (TME), demonstrating the high heterogeneity of GC.³ Currently, surgery and systemic chemotherapy are the mainstays of treatment for GC patients, but the 5-year survival rate of patients with advanced GC remains less than 30%.⁴ Despite the increasingly frequent use of immunotherapies and targeted therapies, the prognosis of patients and the effectiveness of various treatment modalities remain different even within the same GC stage due to significant tumor heterogeneity.^4–6 Therefore, more precise methods for tumor typing and prognostic assessment of GC, along with research into the relationships between different molecular subtypes of GC and patient clinical characteristics, are needed to develop personalized treatment approaches for patients.

As a cancer with high heterogeneity, GC can be divided into distinct subtypes based on morphology and histology.^7–9 With the intensive development of high-throughput sequencing technologies, molecular classification systems have been further proposed by entities such as The Cancer Genome Atlas (TCGA) and the Asian Cancer Research Group (ACRG).^10,11 However, these systems focus mainly on genetic mutations and expression patterns within tumor cells, with little consideration of the composition and functional status of immune cells in the TME and their role in cancer immunotherapy; thus, the prognosis of GC patients cannot yet be accurately assessed.

Tumor-associated neutrophils (TAN) are key components of the TME and have dual roles. On the one hand, TAN can promote tumor development by promoting angiogenesis, extracellular matrix remodeling, metastasis and immunosuppression. On the other hand, anti-tumor TAN kills tumor cells through direct cytotoxicity and activation of adaptive immune response.^12,13 Vimentin (VIM) and RNA binding motif single stranded interacting protein 1 (RBMS1) have been shown to be immunosuppressive molecules, with VIM involved in the interaction between immune cells and endothelial cells in the TME,¹⁴ while RBMS1 has been shown to reflect the unique TME and promote tumor progression in ocular melanoma.¹⁵ Regulator of G-protein signaling 2 (RGS2), an endogenous regulator of a G protein-coupled receptor (GPCR), has been proposed as a novel biomarker for gastric cancer and is associated with CD8+ T cell infiltration.¹⁶ However, the relationship between neutrophils and gastric cancer typing, as well as these three genes, remains unclear.

Therefore, this study aimed to identify immune-related molecular subtypes and their clinical significance in GC. First, by single-cell analysis, we established that neutrophil subpopulations accounted for the highest percentage of the total cell population in GC samples. Then, neutrophil-related GC molecular subtypes (C1 and C2) associated with different prognostic outcomes were identified using high-dimensional weighted gene coexpression network analysis (hdWGCNA), univariate Cox regression analysis and a nonnegative matrix factorization (NMF) algorithm. The unique clinical, immunological and molecular features of these subtypes were further comprehensively characterized. In addition, key prognostic genes VIM, RBMS1 and RGS2 in GC were successfully identified based on the module genes, differentially expressed genes (DEGs), and random survival forest analysis. Finally, the roles and functional mechanisms of the key genes in GC were explored by clinical samples and cellular experiments. In this study, we explored the heterogeneous immune microenvironment phenotypes of GC and the clinical significance of these phenotypes through analysis of multiple datasets and identified distinct neutrophil-related molecular subtypes and key prognostic genes with both prognostic and therapeutic significance. Thus, our findings provide new insights and strategies for immunotherapy and individualized treatment of GC.

Materials and Methods

Data Collection

The original processed mRNA expression data for stomach adenocarcinoma (STAD) were acquired from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/), including 36 normal tissue samples and 412 STAD samples. Subsequently, multiple GC-related datasets were obtained from the Gene Expression Omnibus (GEO) public database, available at https://www.ncbi.nlm.nih.gov/geo/info/datasets.html. First, the GSE163558 dataset, containing data for three primary tumour samples (PT), one adjacent non-tumoural sample (NT), was selected for single-cell analysis. Second, as validation sets, GSE84437 (n = 433), GSE66229 (n = 300) and GSE15459 (n = 192) contained expression profiles and clinical information of STAD patients.

Single-Cell Analysis

First, the expression profile was converted using a Seurat package (version 4.1.1),¹⁷ excluding genes with low expression levels (nFeature_RNA > 300 and percent.mt < 20). Subsequently, the data were subjected to preprocessing steps, including standardization, homogenization, principal component analysis (PCA), and batch effect correction, using Harmony. Through ElbowPlot, the optimal number of principal components (PCs) was determined. Next, the t-distributed stochastic neighbor embedding (t-SNE) method was employed to analyze and initially reveal the positional relationship between each pair of clusters. The Celldex package was used to annotate these cell clusters, associating them with various cell types closely related to disease progression. The FindAllMarkers function was utilized to identify the specific marker genes for each cell subcluster, with the min.pct parameter set to 0.25 and the logfc.threshold parameter set to 1. Finally, key marker genes distinguishing each cell subcluster were screened with min.pct. p_val_adj < 0.05 and |avg_log2FC| > 1.

Ligand–Receptor Interaction Analysis

The CellPhoneDB database (version: 4.0)¹⁸ was employed to assess the significance of the ligand‒receptor interactions identified in the single-cell data. Within this analysis, all of the cells’ cluster labels were randomly assigned with 1000 permutations, and the average expression level of each receptor and the average expression level of each ligand in the interaction cluster in which the receptor interacts with the ligand were calculated. Lastly, several particularly interesting ligand–receptor pairings were chosen for graphic presentation.

hdWGCNA

hdWGCNA package (version 0.2.2)¹⁹ provides built-in functions for network inference, gene module identification, functional gene enrichment analysis, statistical testing of network reproducibility, and data visualization. hdWGCNA is not only suitable for analysis of traditional single-cell RNA sequencing (scRNA-seq) data but can also handle long-read single-cell data for homotypic level network analysis. SetupForWGCNA function was used to construct the co-expression network of cell expression genes, and soft threshold was set at four. The hierarchical clustering tree was generated by the PlotDendrogram function. The GetMEs function was used to extract the expression level of the modular eigengene (ME) of each module. Additionally, module hub genes were identified and acquired by setting the n_hubs parameter to 500.

Gene Functional Enrichment Analysis

The Metascape database (www.metascape.org) was used for comprehensive functional annotation analysis of key gene sets to fully understand their functional relevance. Gene Ontology (GO) term enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for specific gene sets. A minimum overlap ≥ 3 and P ≤ 0.01 were considered to indicate statistical significance.

Classification of Subtypes

Univariate Cox regression analysis was performed with the survival package to assess the correlation between each gene within the module genes and overall survival (OS). Then, the candidate gene set was clustered using unsupervised learning through the NMF program The accuracy and stability of the subtype assignments were verified by using the same set of candidate genes in three external validation sets. mRNA expression data for the abovementioned immune genes were analyzed to identify significant differences in expression between subtypes.

Immune Cell Infiltration Analysis

The microenvironment cell populations (MCP) algorithm was utilized to assess the abundance of eight immune and two nonimmune stromal cell populations. As an alternative method for assessing immune infiltration, single-sample gene set enrichment analysis (ssGSEA) was employed to measure the enrichment of specific gene sets inside a single sample. In addition, six other immune cell populations were estimated using the GSVA-R software package. Furthermore, the ESTIMATE algorithm was applied to calculate immune scores and matrix scores, inferring the degrees of matrix and immune cell infiltration in cancer tissues.

Subtype GSEA

The MSigDB database (http://www.gsea-msigdb.org/gsea/downloads.jsp) provides a large set of predefined genes for subsequent GSEA, encompassing GO term and KEGG pathway gene sets. The log2 (fold change) (log2FC) values for each gene across GC subtypes were computed with the limma package and analyzed with the GSEA function in the clusterProfiler package. Finally, among the results for each subtype, the 10 pathways with the highest enrichment according to the normalized enrichment score (NES) values were selected for visualization.

Performance Verification of Subtypes

The limma package was applied to identify DEGs between GC subtypes, and the criterion for differential gene expression was P < 0.05 after correction (this criterion was used only for the analysis described in Results section 3.5). The 30 genes with the highest log2FC values in each subtype were selected to generate a classifier based on 60 genes for predicting the subtype of new samples. Then, the nearest template prediction (NTP) algorithm was utilized to predict the subtypes associated with the 60 genes in the dataset, and the results were compared with the results of classification based on the NMF algorithm to evaluate the accuracy and consistency of the two methods.

Immunotherapy and Drug Sensitivity Analyses

Subnetwork Mappings in Alignment of Pathways (SubMAP) is an algorithm based on the similarity of gene expression profiles that was used in this study to indirectly predict the potential response of GC subtypes to immunotherapy by quantifying the similarity between our subtypes and those of previously published samples from patients who had undergone immunotherapy. In addition, We used the R software package “pRRophetic” for tumor drug sensitivity prediction. The half-maximal inhibitory concentration (IC50) values for specific chemotherapeutic agents were predicted by regression analysis, and the accuracy of the regression model and the predictions was tested.

WGCNA

Using the WGCNA-R package, coexpression networks of the top 5000 genes based on variance in the bulk RNA-seq dataset were constructed, with a soft threshold of 3. Based on the correlations in gene expression, a weighted adjacency matrix was constructed and further transformed into a topological overlap matrix (TOM) to estimate network connectivity. The TOM was used as input for hierarchical clustering to construct a dendrogram. In the dendrogram, Different branches represent different gene modules, and genomes of the same module have similar expression patterns. Through clustering, genes were divided according to their expression patterns into several modules represented by different colors.

Random Survival Forest Screening of Key Genes

The randomForestSRC function was used for feature selection. The Monte Carlo simulation was run with 500 iterations (nrep = 500), which helped to stabilize the model and assess the importance of the variables. Next, the importance score of each gene was extracted and ranked, and the genes with a relative importance score greater than 0.4 were selected as the final key genes.

Analysis of Regulatory Networks of Important Genes

The R package “RcisTarget” was used to predict transcription factors by identifying transcription factor binding motifs that were enriched in specific gene sets. A NES was calculated based on each motif in the database. To estimate the enrichment of each motif in a specific gene set, the area under the curve (AUC) for each pair of motif–motif sets was first calculated by the recovery curve of the sequence based on the gene set. NESs were calculated to identify enriched motifs based on the AUC distribution of all motifs in the gene set. We used rcistarge.hg19.motifdb.cisbpont.500bp for the Gene-motif rankings database.

Clinical Specimen Collection

Thirty pairs of GC tissues and adjacent normal tissues were obtained from the First Affiliated Hospital of Ningbo University between 2022 and 2024. None of the patients received preoperative chemotherapy or radiotherapy, and the diagnosis of all samples had a clear pathologic basis. The tissues were immersed in RNA-fixer Reagent (Bioteke, Beijing, China) and then stored in a refrigerator at −80°C until RNA was extracted. All patients did not receive radiotherapy, chemotherapy, or other treatments prior to surgery. The Ethics Committee of the First Affiliated Hospital of Ningbo University approved this study (IRB No. KY20220101). All participants signed an informed consent form and were informed of the details of the study.

Cell Culture and Transfection

GES-1, HGC-27 and NUGC-3 were acquired from the Shanghai Institute of Biochemistry and Cell Biology of the Chinese Academy of Sciences. All cells were cultured in RPMI-1640 medium (Invitrogen, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (FBS) at 37°C in 5% CO₂. Cells were inoculated into 6-well plates and transfected when the growth density reached 60%. The siRNA was transfected into cells according to the instructions of Lipofectamine 2000 (Invitrogen), and the cells were collected after 24 h. Quantitative real-time polymerase chain reaction (qRT‒PCR) was performed to verify the transfection effect. The sequences of specific siRNAs were constructed by GenePharma Co., Ltd. (Shanghai, China) and listed in Table S1.

Plate Colony Formation Assay

The transfected cells were cultured for 24 h, then inoculated into 6-well plates at a density of 500 cells/well and cultured for 14 d. After washing with PBS for 2 times, the cells were fixed with 4% paraformaldehyde (Solarbio, Beijing, China) for 15 min, stained with 0.1% crystal violet and photographed for colony counting.

Transwell Migration Assay

Cells were collected 24h after transfection, resuspended with serum-free RPMI 1640 medium and then 200 uL of suspension (1~1.5×10⁵ cells/mL) was added to the upper chamber, and then the bottom chamber was filled with 600 uL of medium containing 30% FBS. The cells were incubated in a cell culture incubator for 24 h, fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and the images were collected under an inverted microscope.

RNA Extraction and qRT‒PCR Analysis

Total cellular RNA was extracted using TRIzo1 (Ambion, Carlsbad, CA), and cDNA synthesis was performed using the GoScript RT system (Promega, Madison, WI, USA) according to the manufacturer’s instructions. qRT‒PCR was performed using the GoTaq-qPCR Master Mix (Promega) kit with cDNA as the template with 40 cycles of amplification. Four replicate wells were set up for each sample, and relative gene expression was calculated by the 2^−ΔΔCt method with Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the internal reference gene. The sequences of the primers used for qRT‒PCR in this study are listed in Table S1.

Statistical Analysis

The R language (version 4.0) or Statistical Product and Service Solutions (SPSS) 19.0 software was used for all analyses. P < 0.05 was statistically significant.

Results

Preliminary Processing of Single-Cell Expression Profile Data

This analysis focused on the expression profiles of GC-related tissue samples. In this analysis, only cells meeting the criteria of nFeature_RNA > 300 and percent.mt < 20 within the expression profile were retained. The feature expression levels of 14,976 cells were included in the subsequent analysis (Figure S1A and B). Gene expression profiles within the sample were characterized, highlighting the 5 genes with the highest normalized variance (Figure S1C).

Cluster Analysis of Subtypes in Single-Cell Sample Subtypes

Through dimensionality reduction of the data for the 20 selected genes by PCA, it was found that each gene exhibited unique score values across various dimensions (Figure S1D). However, dimensionality reduction between samples revealed batch effects between samples (Figure S1E). Harmony analysis was subsequently applied to further eliminate batch effects and minimize dimensionality (Figure S1F). Through ElbowPlot, the optimal number of PCs was determined to be 10 (Figure S1G). Then, 18 subtypes were acquired through t-SNE (Figure 1A).

Figure 1 Cell type annotation based on single-cell data. (A) t-SNE analysis resulted in the identification of 18 cell clusters. (B) The 18 clusters were clustered into 8 cell types with the R package SingleR. (C) Proportions of neutrophil subtypes in tumor samples and control samples.

Annotation of Clustering Subtypes

We annotated each subtype using the R package SingleR, representing the eighteen subtypes to neutrophils, CD8+ T cells, epithelial cells, B cells, monocytes, fibroblasts, endothelial cells, and NK cells (Figure 1B). Among these subtypes, the neutrophil subtype accounted for the largest proportion of cancer samples (Figure 1C). Further subsets of the neutrophil subtype were analyzed sequentially by PCA, Harmony, ElbowPlot and FindClusters (Figure S2A–D). The percentage of cancer samples in subtypes C1, C2, C3, and C4 was greater than 99.64% (the median proportion), leading to the classification of these subtypes as cancer cell clusters (STAD_Clu) (Figure S2E). Furthermore, trajectory analysis demonstrated that these newly defined cancer cell clusters exhibited distinct branches (Figure S2F–I).

Identification of receptor‒ligand Interaction Pairs

The software package CellPhoneDB was used to analyze the ligand‒receptor interactions of features within the single-cell expression profile. Subsequently, specific pairs of ligands and receptors were chosen for presentation (Figure 2A and B). The Endothelial_cells|Epithelial_cells and Monocytes|Fibroblasts clusters exhibited high interaction scores with CD74_APP and CD74_COPA.

Figure 2 Analysis of receptor‒ligand pairs. (A) Bubble plot showing scores for ligand‒receptor interactions between cells. (B) Cell network displaying the strength of interactions between specific cell types.

hdWGCNA

hdWGCNA was performed to construct the coexpression network of genes in the neutrophil cluster. The parameter group.by was set to cluster subtype and custom subtype to facilitate the construction of the coexpression network and the exploration of biomarkers for disease development. The soft threshold was determined with the function “TestSoftPowers” and was set to 4 (Figure 3A). In our research, two gene modules were identified: the turquoise and blue modules (Figure 3B). Further analysis between modules and ME levels (Figure 3C and D) revealed that the ME levels in the turquoise module, were higher in tumor cell subtypes (Figure 3E). In the turquoise module, genes with a percent expression value > 75 and an average expression value > 1.5 were identified. Analysis of this module using the GetHubGenes function revealed 500 module genes.

Figure 3 hdWGCNA reveals key gene modules in the neutrophil subtypes. (A) Calculation of pairwise correlations of input features with soft-threshold weighted correlations (β). (B) Cluster dendrogram showing two key gene modules: turquoise and blue. (C and D) Visualization of the differences in ME expression in each module across subtypes. (E) Expression characteristics of the turquoise and blue module genes in the neutrophil subtypes.

Functional Analysis of Candidate Gene Sets

Pathway analysis of the 500 module genes was conducted using the Metascape database, revealing that these module genes predominantly exhibited enrichment in the ficolin-1-rich granule, secretory granule membrane, focal adhesion and cell activation pathway, among other pathways (Figure S3A). Concurrently, protein interaction network analysis was performed on the genes in this gene cluster via Cytoscape software (Figure S3B).

Identification of Prognostic Genes and Molecular Subtypes

In order to identify prognostic genes among the module genes, clinical information from GC patients was collected, and GC feature genes were screened via the univariate Cox regression feature selection algorithm. Univariate Cox regression analysis identified 20 prognostic genes (P < 0.01) (Table S2). The NMF consensus clustering method was applied to the TCGA dataset containing GC samples; the samples according to the expression profiles of the 20 candidate genes, and the optimal k value was determined. Following comprehensive evaluation, k = 2 was chosen as the most accurate classification (Figure 4A). Subsequent independent validation of the dataset of GC samples with three external validation sets using the previously mentioned K2 classification similarly revealed two distinct molecular subtypes. A significant difference in prognosis was observed in the TCGA dataset, with C1 associated with a better survival probability than C2 (Figure 4B). In addition, a similar difference was observed in the subtype analysis of the three validation sets (GSE84437, GSE66229, and GSE15459), with the survival time in the C1 subgroup significantly longer than that in the C2 subgroup (Figure 4C–E).

Figure 4 Identification of GC prognosis-related genes and molecular subtypes. (A) NMF consensus clustering analysis showed that k = 2 was the optimal number of clusters. (B) Survival analysis based on the TCGA dataset. (C–E) Analysis of three external validation sets confirmed the differences in prognosis between the subtypes.

Analysis of The immunological Features of the Subtypes

Further investigations were conducted to determine whether the subtypes exhibited distinct immune characteristics. Initially, the immune response was quantified through ssGSEA (Figure 5A), and differential expression analysis was then performed to identify subtype-specific immune signatures. The results indicated clear differences in immune characteristics between the two subgroups, with significantly greater numbers of B cells, T cells, and T helper cells in the C1 subgroup than in the C2 subgroup (Figure 5B). Furthermore, the NRF2 and PI3K pathway activation scores of immune cells were greater in the C1 subtype than in the C2 subtype (Figure 5C). Subsequently, we employed the ESTIMATE algorithm to calculate immune and matrix scores. Differences were observed in these scores between the two subgroups, with the scores of C1 significantly elevated compared to those of C2 (Figure 5D and E).

Figure 5 Immunological characterization of the GC subtypes. (A) Pathway-specific heatmap for metaclusters in TCGA. (B) Differences in the immune characteristics of the subtypes. (C) Box plots showing the differences in key pathway activation scores between the subtypes. (D) Analysis of differences in immune scores between the subtypes. (E) Analysis of differences in stromal scores between the subtypes. (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

Evaluation of Immune Infiltration by the MCP Algorithm and ssGSEA

The MCP-counter tool and ssGSEA were employed to investigate concentrated immune infiltration based on immune-related gene expression. The abundances of 16 immune-related cell types were visualized in a heatmap (Figure 6A). The levels of immune infiltration differed significantly between the two subgroups, suggesting that the immune characteristics of the classification were reproducible. In addition, the infiltration of Th2 cells, Th17 cells, and Tregs was significantly greater in the C1 subtype than in the C2 subtype, while the results for B cells and Tregs were consistent with the previous enrichment analysis results for immune-related features (Figure 6B). Further study of the associations between the subtypes and the expression of 12 potentially targetable immune checkpoint genes selected based on inhibitor drugs that are licensed for use in particular cancer types or are presently undergoing clinical trials. It was indicated that C1 exhibited higher expression of LAG3, CD274, PDCD1, CD247 and other immunotherapy targets (Figure 6C).

Figure 6 Assessment of subtype immunity levels and molecular characterization. (A) Heatmap showing the differences in the abundance of 16 immune-related cell types assessed using the MCP-counter algorithm and ssGSEA between two GC subtypes, C1 and C2. (B) Box plots showing the differences in the abundance of immune cell and stromal cell populations between the subtypes. (C) Expression levels of immune checkpoint genes in the two clusters. p values are indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001. (D) GO analysis of pathways with higher subtype scores. (E) KEGG analysis of pathways with higher subtype scores.

Evaluation of the Molecular Characteristics of the Subtypes

We further quantitatively analyzed GO term and KEGG pathway enrichment via ssGSEA and identified several related pathways with obvious differences in enrichment between the two subgroups. Among these pathways, the main signaling pathways with high C1 subtype scores in the GO term analysis were GOCC IMMUNOGLOBULIN COMPLEX and GOMF ANTIGEN BINDING (Figure 6D). According to the KEGG pathway analysis, the pathways with relatively high C1 subtype scores were AUTOIMMUNE THYROID DISEASE and GRAFT VERSUS HOST DISEASE, among others (Figure 6E).

Performance Verification of the Molecular Subtypes

We ranked the genes in each subtype according to log2FC values and selected the top 30 genes, which had the most significant changes in expression in that subtype, to generate a 60-gene classifier with high accuracy and potential for clinical application (Figure 7A and B). By utilizing drug sensitivity data from the GDSC database and the R package “pRRophetic”, the chemosensitivity of each cancer sample was predicted to further investigate the sensitivity of the subtypes to common chemotherapy drugs. Compared with C2 subtype, C1 subtype showed significant differences in sensitivity to imatinib, pazopanib, midostaurin, dasatinib and other drugs (Figure 7C). Further exploration of mutation profiles across the subtypes revealed a significantly greater proportion of TTN mutations in patients with the C1 subtype than in those with the C2 subtype (Figure 7D). Moreover, multigroup comparisons revealed significant differences in neoantigen levels and mutation burdens between the subtypes (Figure 7E and F). Moreover, the sensitivity of the two subtypes anticancer immunotherapy was further predicted based on an immunotherapy dataset, revealing that the C1 subtype was associated with greater sensitivity to immunotherapy (Figure 7G). Analysis of tumor immune dysfunction and rejection showed differences between subtypes, among which the Responder and Exclusion statuses showed significant differences between the subtypes (Figure S4A–D).

Figure 7 Identification of classifiers and prediction of immunotherapy and chemotherapeutic drug sensitivity. (A) Heatmap of the expression levels of the 60 genes in the classifier. (B) Consistency between the NMF and NTP algorithms. (C) Associations between GC subtypes and targeted drug sensitivity. (D) Heatmap of mutations in GC subtype feature genes. (E and F) Differences in neoantigen levels and mutation burdens between the GC subtypes. (G) Prediction of sensitivity to antitumor immunotherapy.

Coexpression Network Determined by RNA-Seq Analysis in GC

WGCNA was performed to construct the coexpression network of genes in the GC cohort. The molecular subtypes were further utilized as traits to construct WGCNA networks and explore biomarkers of GC development (Figure 8A and B). The soft threshold β was determined as 3 by the function “sft$powerEstimate”. Subsequently, gene modules were identified according to the TOM. A total of 9 gene modules were identified in this analysis: the black (n = 243), blue (n = 1061), green (n = 308), magenta (n = 179), pink (n = 198), purple (n = 114), red (n = 261), turquoise (n = 2211), and yellow (n = 425) modules. The module characteristic genes of the two subtypes were analyzed, and the expression of the 9 module feature genes was found to differ significantly between the subtypes (Figure 8C and D). Further analysis of the relationships between modules and traits revealed that the turquoise module had the highest correlation with sample traits (cor = 0.43, P = 4e-19) (Figure 8E and F). Functional analysis of the 2,211 genes in the turquoise module revealed the enrichment of these genes in pathways such as extracellular matrix, tube morphogenesis, tissue morphogenesis, and epithelial cell differentiation (Figure 8G).

Figure 8 WGCNA identifies key coexpression modules in the GC RNA-seq data. (A) Sample clustering to identify outliers. (B) Cluster dendrogram of coexpression network modules. The different colors represent different modules. (C and D) The expression of characterized genes in all 9 modules differed significantly between the subtypes. *P < 0.05, ****P < 0.0001. (E) Correlation analysis between gene modules and sample traits. (F) Correlation analysis between gene significance for group and module membership in the turquoise module. (G) Pathway enrichment analysis of genes in the turquoise module.

Differential Expression Analysis Between the Subtypes

The expression levels of genes were compared between the two subgroups via the R package limma (|logFC| > 0.585 and adj. p.AL < 0.05). Differential expression analysis between the C2 and C1 subgroups revealed 220 downregulated genes and 477 upregulated genes, for a total of 697 DEGs (Figure 9A). Six of these genes were identified as prognostic genes by subtype analysis (Figure 9B). To further identify key genes affecting GC, random survival forest analysis was performed on these 6 genes in the TCGA-STAD cohort, and genes with a relative importance score > 0.4 were identified as the final marker genes (Figure 9C). Kaplan–Meier analysis showed significant correlations for three of these genes: VIM, RBMS1 and RGS2 (Figure 9D–F). The expression of these three key genes at the single-cell level was also determined (Figure 9G–I).

Figure 9 Differential expression analysis between subtypes and identification of key prognostic genes in GC. (A) Volcano plot of differentially expressed genes. The blue dots represent downregulated genes; the gray dots represent genes with no difference in expression; and the red dots represent upregulated genes. (B) Venn diagram showing the identified prognostic genes. (C) Random survival forest analysis was used to filter key genes and rank them by importance with the criterion of relative importance > 0.4. (D–F) Kaplan‒Meier survival analysis showed that the expression of RBMS1 (D), RGS2 (E), and VIM (F) was significantly associated with the survival rate of patients with GC. (G) The level of RBMS1, RGS2 and VIM in different identity. (H) The expression distribution of marker genes. (I) Bubble plot displaying the normalized mean level of markers.

Analysis of the Regulatory Network of the Key Genes

We selected 3 key genes for analysis and found that they are regulated by common mechanisms such as multiple transcription factors. To gain a deeper understanding of how these transcription factors regulate the target genes, we performed enrichment analysis using cumulative recovery curves (Figure S5A). The results showed that the motif with the highest Normalized Enrichment Score (NES) was cisbp__M1503, with an NES of 5.94. 2 genes were significantly enriched in this motif, RBMS1 and RGS2. Finally, all the enriched motifs and the corresponding transcription factors were presented (Figure S5B).

Validation of the Differential Expression and Biological Function of the Key Genes in GC

To verify the differential expression, we examined the expression levels of VIM, RBMS1 and RGS2 in GC cell lines and tissues using qRT-PCR. The results showed that compared with the normal gastric epithelial cell line GES-1, the expression of VIM and RBMS1 was significantly up-regulated and the expression of RGS2 was significantly down-regulated in GC cells (Figure 10A). Consistent with the cellular results, the expression levels of VIM and RBMS1 were significantly upregulated and the expression level of RGS2 was significantly downregulated in GC tissues (Figure 10B). Combined with the clinicopathological characteristics of GC patients, VIM levels in GC tissues were significantly correlated with tumour diameter (P = 0.016) (Table S3), RBMS1 was correlated with tumour diameter (P = 0.038) and venous invasion (P = 0.038) (Table S4), and RGS2 was correlated with tumour diameter (P = 0.003) (Table S5).

Figure 10 Validation of differential expression and biological functions of key genes in GC. (A) Expression of VIM, RBMS1 and RGS2 in GC cells. (B) Expression levels of VIM, RBMS1 and RGS2 in GC tissues and adjacent normal tissues (n = 30). (C) Measurement of transfection efficiency of si-RBMS1 and si-RGS2. (D) Effect of silencing RBMS1 and RGS2 expression on GC cell proliferation. (E) Effect of silencing RBMS1 and RGS2 expression on GC cell migration. (F) Effect of silencing RBMS1 and RGS2 expression on GC cell colony formation. Data are expressed as mean ± SD, n = 3. *P < 0.05, **P < 0.01, ***P < 0.001.

To further validate the biological functions of the key genes, we selected RBMS1 and RGS2 for functional validation. The siRNAs of RBMS1 and RGS2 were transfected in two GC cells, AGS and HGC-27. qRT-PCR showed that the gene expression levels of both the si-RBMS1 and si-RGS2 groups were lower than those of the si-NC group (Figure 10C). Then, cellular experiments showed that silencing the expression of RBMS1 significantly inhibited the proliferation and migration of GC cells. On the contrary, silencing the expression of RGS2 significantly promoted the proliferation and migration of GC cells (Figure 10D–F).

Discussion

Recently, studies based on molecular typing have become increasingly critical in the diagnosis, treatment and prognostic assessment of cancer. For example, Dai et al identified subtypes of head and neck squamous cell carcinoma (HNSCC) with distinct prognoses and treatment vulnerabilities.²⁰ Chen et al developed a classification system based on metabolic gene expression profiles to characterize the molecular features of hepatocellular carcinoma (HCC).²¹ Wu et al categorized colorectal cancer (CRC) samples into four subtypes with distinct TME characteristics through comprehensive multiomics analysis and artificial intelligence-powered spatial analysis.²² Overall, these and other such studies are dedicated to providing clinicians with more precise diagnostic tools and therapeutic strategies through in-depth exploration of the molecular features and heterogeneity of tumors, as well as to providing new ideas and directions for drug development and clinical trial design. Although numerous genome- and transcriptome-based GC classifications have been proposed in recent years,^10,11 no consensus has been reached on the molecular taxonomy associated with tumor immunity in GC. In this study, starting with immune cell subpopulations, we combined scRNA-seq and bulk RNA-seq datasets, identified immune molecular subtypes of GC by hdWGCNA and machine learning algorithms, investigated the clinical significance of the different subtypes, and then identified and validated key prognostic genes associated with these GC molecular subtypes.

The landscape of cancer treatment is shifting from traditional methods to precision medicine, with detailed immunotyping allowing researchers to identify the functions and interactions of different immune cells in the TME, thus providing crucial insight for developing new treatment strategies.^23,24 In this study, eight cell categories, dominated by neutrophils, were identified by analysis of the scRNA-seq dataset (Figure 1A). Neutrophils, a crucial component of the TME, play dual roles in tumor promotion and inhibition in GC.^25–27 By regulating key molecular markers, neutrophils are expected to contribute significantly to antitumor immunity, increasing the efficacy of cancer immunotherapy.^12,25 Therefore, we performed an in-depth analysis of neutrophil subtypes and identified four STAD clusters with significant tumorigenic properties (Figure S2E). Additionally, there were frequent and close interactions and communication networks connecting the cell subgroups, which may facilitate the discovery of new therapeutic targets (Figure 2A and B). hdWGCNA can be used to construct coexpression networks based on scRNA-seq data, aiming to identify coexpressed gene modules, explore the associations between the gene coexpression network and phenotypes, and reveal the core genes in the network interactions.²⁸ Our study revealed that the turquoise module was significantly associated with tumor cell subtype and is a potential biomarker related to neutrophils (Figure 3E). Notably, the key genes of this module have a wide range of biological functions involving several key pathways, including the ficolin-1-rich granule, secretory granule membrane, focal adhesion, and cell activation pathways (Figure S3A), which provides insight into the biological regulatory mechanisms of GC. Module genes can be used to classify tumors into distinct molecular subtypes based on gene expression patterns, helping to reveal the biological differences in tumors and predict disease progression and treatment responses.²⁹ Based on the prognostic genes among the module genes, our two GC molecular subtypes (C1 and C2) were associated with significant differences in prognosis, with the C1 subtype associated with a better survival probability (Figure 4B). This finding provides the foundation for the development of precision medicine and personalized treatment strategies.

We conducted an in-depth comprehensive analysis of the molecular subtypes of GC, aiming to explore the potential mechanisms of the subtypes in impacting the prognosis of GC patients. First, subtypes C1 and C2 significantly different immune characteristics. Compared with those of the C2 subtype, the immune characteristics of B cells, T cells, T helper cells, etc, were significantly more prominent in the C1 subtype (Figure 5B). Furthermore, C1 exhibited higher NRF2 and PI3K pathway activation scores and a higher immune score than C2 (Figure 5C), a finding that may have significant implications for the development of new treatment strategies targeting the immune system in patients with GC.^21,30 Second, the better prognosis of patients with the C1 subtype may be attributed to the combined effect of the high infiltration rates of Th2 cells, Th17 cells, and Tregs and high expression levels of immune checkpoint genes (Figure 6B and C). Immunotherapy targeting immune checkpoints has recently become a new strategy for GC treatment. For example, by inhibiting specific immune checkpoints such as PD-1, the T-cell response to tumor cells can be enhanced, resulting in anticancer effects.³¹ Third, multiple pathways related to immune processes differed significantly between the two subtypes (Figure 6D and E). The C1 subtype was associated with high activity in pathways such as GOCC IMMUNOGLOBULIN COMPLEX and GRAFT VERSUS HOST DISEASE, revealing the active immune response and autoimmune characteristics of this subtype. Understanding the immune pathways that are specifically upregulated in distinct subtypes can provide valuable insight for the development of treatments targeting specific subtypes.

Accurate subtype classification can facilitate the development of personalized therapeutic regimens targeting specific subtypes and improve the treatment outcome and quality of life of GC patients. In this study, we developed a predictive classification model based on key DEGs between different subtypes of GC that demonstrated excellent advantages in terms of accuracy and clinical application potential (Figure 7A and B). First, compared with C2 subtype, C1 subtype showed significant differences in sensitivity to imatinib, pazopanib, midostaurin, dasatinib and other drugs. Although tyrosine kinase inhibitors such as imatinib have limited anti-gastric cancer effects when used alone, studies have shown that they can be used as effective chemosensitizing agents to target the PDGF/PDGFR signaling pathway in tumor cells and stromal cells in disease progression and angiogenesis, thereby enhancing the efficacy of mainstream anticancer drugs such as 5FU.³² In addition, the proportion of mutations in genes such as TTN was significantly greater in patients with the C1 subtype than in patients with the C2 subtype, and the neoantigen levels and mutation burdens were also significantly different between subtypes. The higher neoantigen levels and mutation burdens may imply that the tumors of patients with the C1 subtype were more immunogenic and thus more likely to stimulate an immune response.³³ Tumor rejection usually involves the recognition and clearance of tumor cells by immune cells, such as T cells.³⁴ Our study demonstrated that the C1 and C2 subtypes also differed in terms of tumor rejection and response, with the C1 subtype more sensitive to immunotherapy. Finally, to identify the coexpression network of genes in the GC cohort, we conducted WGCNA and identified the turquoise module as the module most strongly related to GC classification (Figure 8E and F); this module is thus a potential biomarker of GC progression. Therefore, our model is highly important for gaining an in-depth understanding of biological mechanisms, discovering new therapeutic targets, and developing new drugs.

The development of GC involves alterations in the expression of multiple genes that may predict cancer progression or changes in tumor cell function, resulting in heterogeneous tumor subtypes.³⁵ Here, we successfully identified three key prognostic genes, VIM, RBMS1, and RGS2, in neutrophil-related molecular subtypes of GC (Figure 9D–F). VIM encodes a type III intermediate filament protein that is highly expressed in a variety of cancers and serves as a biomarker for epithelial–mesenchymal transition (EMT), which is involved in cell migration, invasion, adhesion, and signaling.^36,37 Multiple lines of evidence suggest that VIM is also a potential tumor marker, exhibiting potent therapeutic targeting ability.^36,38,39 RBMS1 belongs to the single-stranded DNA/RNA binding protein family, is involved in translation regulation, tumor immunity, cancer metastasis, and other biological processes, and has potential as a new therapeutic target.^40–42 RGS2 is aberrantly expressed in various cancers and is significantly associated with clinical prognosis.^43,44 RGS2 has various functions, including promoting angiogenesis, regulating apoptosis, mediating cancer cell dormancy, and promoting tumor recurrence.^45–47 In summary, VIM, RBMS1, and RGS2 can individually participate in tumor development through various mechanisms. Transcription factors can control the transcription process of chromatin and participate in regulating gene expression, thus playing an essential role in the pathogenesis of GC.⁴⁸ In our study, VIM and RBMS1 were significantly overexpressed in GC cells and tissues, whereas RGS2 was significantly downregulated in GC cells and tissues (Figure 10A and B). Furthermore, we verified the biological functions of RBMS1 and RGS2 in gastric carcinogenesis, demonstrating good performance as potential therapeutic targets for GC. Therefore, VIM, RBMS1, and RGS2 are key genes related to neutrophils in GC and could play important roles in the prognosis or treatment of this disease. In our study, it was found that multiple transcription factors such as NR2F1 and FOXP3 co-regulate the expression of VIM, RBMS1 and RGS2, showing good performance as potential therapeutic targets for GC. The abnormalities of these transcription factors are closely related to the occurrence of GC, and may participate in the evolution of GC by regulating the immune microenvironment, directly regulating tumor cells, and participating in the regulation of signaling pathways.^49,50

We used scRNA-seq and hdWGCNA to identify for the first time two molecular subtypes of neutrophil-associated GC with different prognostic outcomes and clinical significance. It is worth mentioning that compared with other reported molecular subtypes of GC, our proposed neutrophil-associated GC subtype achieved significant breakthroughs and integration at multiple levels. Specifically, our study not only in-depth explored the complex relationship between GC and immune infiltration, but also systematically analyzed the differences between these two subtypes in immune checkpoint expression, signaling pathway regulation, and tumor mutation burden. More importantly, we successfully correlated these neutrophil-associated GC subtypes with sensitivity to immunotherapy and chemotherapeutic agents. On this basis, three new key genes (VIM, RBMS1 and RGS2) related to the prognosis of GC patients were further screened from the differentially expressed genes of the two GC subtypes. These results are of great significance for predicting the success of GC immunotherapy and making treatment decisions.

Conclusion

In conclusion, in this study, two GC neutrophil-related molecular subtypes with different prognostic outcomes and clinical significance were identified. On this basis, three new key genes (VIM, RBMS1 and RGS2) related to the prognosis of GC patients were further screened from the differentially expressed genes of the two GC subtypes. Overall, these results may provide guidance for the advancement of innovative immunotherapy approaches for GC and improve patient prognosis.