Gene expression data for stomach adenocarcinoma (STAD) were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), comprising expression profiles from 448 patients. The GSE26901 Series Matrix File was downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html), with the annotation file being GPL6947, and included expression data from 109 patients. Similarly, the GSE62254 Series Matrix File was downloaded from GEO, annotated with GPL570, and contained expression data from 300 patients. In addition, single-cell RNA sequencing data from GSE183904 were retrieved from GEO, encompassing 36 samples with complete single-cell expression profiles, including 10 controls and 26 tumor samples.

Gastric cancer subtypes were classified through consensus clustering based on the expression profiles of ArMGs. The clustering process was repeated 50 times, with 90% of the samples included in each iteration. The optimal number of clusters was determined by locating the inflection point in the cumulative distribution function curve of the consensus score.

Differential expression analysis among the three sample groups was conducted using the R package “Limma”. Differentially expressed genes (DEGs) were identified with the thresholds of P < 0.05 and |log₂FC| > 0.585. Volcano plots and heatmaps were generated to visualize the DEGs.

Intersecting genes were used to construct a prognostic model using the Lasso regression. A risk score was calculated for each patient by weighting the expression values of the selected genes with their corresponding Lasso regression coefficients. Patients were stratified into low-risk and high-risk groups based on the median risk score. Survival differences between the two groups were assessed using Kaplan-Meier analysis and compared using the log-rank test. The predictive value of the risk score was further evaluated using Lasso regression and stratified analyses. The accuracy of the model was examined using receiver operating characteristic (ROC) curve analysis.

A points-based scoring system was developed using a multivariate regression model. Points were assigned to each level of the influencing factors according to the magnitude of their regression coefficients, reflecting their contribution to the outcome variable. The total score for each patient was obtained by summing these points, which allowed the prediction of the clinical outcome (Cheng et al., 2024).

The chemotherapeutic sensitivity of each tumor sample was predicted using the R package “oncoPredict”, with reference to the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). A regression approach was applied to estimate the IC50 values for each chemotherapy drug, and the GDSC training set was used for 10-fold cross-validation to assess regression and prediction accuracy. All analyses were conducted with default parameters, including the ComBat algorithm for batch effect correction and averaging of expression values for duplicate gene symbols.

In this study, tumor mutational burden (TMB) was calculated as the mutation frequency, defined by the ratio of the number of variants to the exon length for each tumor sample. Microsatellite instability (MSI) scores for TCGA patients were obtained from a previously published study (Bonneville et al., 2017). Patient-specific neoantigens were predicted using NetMHCpan v3.0 (Huang and Fu, 2019).

The CIBERSORT algorithm was used to estimate the relative proportions of 22 immune cell types within the tumor microenvironment based on patient gene expression data. Correlation analyses were then performed to evaluate associations between gene expression levels and immune cell abundancies.

Gene sets were obtained from the Molecular Signatures Database (MSigDB, version 7.0). GSVA was performed to calculate enrichment scores for each gene set in individual samples, enabling the evaluation of potential biological functional differences across sample groups.

Patients were stratified into high-risk and low-risk groups based on the prognostic model. GSEA was then conducted to identify signaling pathways that were differentially enriched between the two groups. Background gene sets for pathway annotation were sourced from the MSigDB database. Significantly enriched gene sets (adjusted P < 0.05) were ranked according to their enrichment scores. GSEA was employed to link tumor classification with underlying biological functions.

The single-cell RNA sequencing data were first preprocessed and normalized using the Seurat package. The t-distributed stochastic neighbor embedding (t-SNE) algorithm was then applied to visualize the spatial distribution of distinct cell clusters. Cell types within the tissues were annotated based on known marker genes retrieved from the CellMarker database and supported by published literature.

Human GC cell lines (AGS and HGC-27) were obtained from the Cell Bank of the Chinese Academy of Sciences or the Shanghai Institute of Biochemistry and Cell Biology. Cells were cultured in RPMI-1640 medium (Invitrogen, Grand Island, USA) supplemented with 10% fetal bovine serum (FBS) and maintained at 37 °C in a humidified atmosphere with 5% CO₂.

Actively growing cells were digested and seeded into 6-well plates. Once cells reached approximately 50%–60% confluence, transfection was performed according to the manufacturer’s instructions using Lipofectamine 2000 reagent (Invitrogen, CA, USA) and Opti-MEM I reduced-serum medium (Gibco, CA, USA). Cells were transfected with either siRNA (100 nM) or negative control siRNA (si-NC, 100 nM). After 6 h of transfection, the medium was replaced with RPMI-1640 supplemented with 10% FBS.

A total of 30 paired gastric cancer tissues and adjacent normal tissues were collected from surgical patients at the First Affiliated Hospital of Ningbo University from 2022 to 2024. None of the patients had undergone radiotherapy or chemotherapy prior to surgery, and all cancer tissues were confirmed by pathological examination. This study was performed in accordance with Chinese clinical guidelines and relevant regulations, and received approval from the ethics committee of the First Affiliated Hospital of Ningbo University (IRB No. KY20220101). Written informed consent was obtained from all patients.

Total RNA was extracted using TRIzol Reagent (Thermo Fisher Scientific, USA). The extracted RNA was reverse transcribed into complementary DNA (cDNA) using a commercially available reverse transcription kit. Then, quantitative PCR was performed using the GoTaq qPCR master mix (Promega). Thermal cycling conditions included initial pre-denaturation at 95 °C for 5 min, followed by 40 cycles of denaturation at 94 °C for 15 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal control. Relative gene expression was quantified using the ΔCt method, where higher ΔCt values indicate lower expression levels. The primer sequences for the housekeeping gene GAPDH and the key genes ODC1 and ALDH18A1 are listed below: GAPDH(F): 5ʹ-ACCCACTCCTCCACCTTTGAC-3ʹ, (R): 5ʹ-TGTTGCTGTAGCCAAATTCGTT-3ʹ. ODC1(F): 5ʹ-TTTACTGCCAAGGACATTCTGG-3ʹ. (R): 5ʹ-GGAGAGCTTTTAACCACCTCAG-3ʹ. ALDH18A1(F): 5ʹ-GCCCTTCAACCAACATCTTCT-3ʹ. (R): 5ʹ-AGGGGTACAGTGATAAACGGG-3ʹ.

After 24 h of cell transfection, cells were seeded into 96-well plates at a density of 1 × 10⁴ cells per well. At 0, 24, 48, 72, and 96 h of incubation, 10 µL of CCK-8 reagent (Beyotime, Shanghai, China) was added to each well and incubated at 37 °C for 1 h. Optional density (OD) at 450 nm was measured using a multifunctional microplate reader (Thermo Fisher Scientific, MA, USA). Each sample was tested in six replicates.

After 24 h of cell transfection, cells were trypsinized into single-cell suspensions and seeded into 6-well plates at a density of 600 cells per well. After 14 days of incubation, each well was washed twice with phosphate-buffered saline (PBS), fixed with 4% paraformaldehyde (Beyotime, Jiangsu, China) for 20 min, and stained with 0.1% crystal violet solution (Beyotime, Jiangsu, China) for 20 min. The wells were then rinsed with double-distilled water until excess crystal violet was removed. Finally, the wells were photographed, and the number of stained colonies was counted. Each assay was performed in triplicate.

Transwell invasion assays were performed using transwell chambers (Corning, NY, USA) equipped with 8.0 μm-pore polycarbonate membranes coated with Matrigel (100 µg/mL; BD, NJ, USA). After 24 h of cell transfection, cells were resuspended in serum-free RPMI-1640 medium at a concentration of 1 × 10⁵ cells/mL. A total of 500 µL of RPMI-1640 medium containing 40% FBS was added to the lower chamber, while 200 µL of the cell suspension was added to the upper chamber. After incubation for 24–48 h at 37 °C, the membranes were washed with PBS, and the cells on the lower surface were fixed with paraformaldehyde for 20 min and stained with 0.1% crystal violet solution for 20 min. Each assay was performed in triplicate.

All statistical analyses were performed using R software (version 4.4.0) (https://www.r-project.org/) or Statistical Package for the Social Sciences (SPSS) 19.0. A P-value < 0.05 was considered statistically significant.

Ten ArMGs were initially retrieved from the GSEA database (https://www.gsea-msigdb.org/gsea/msigdb). A co-expression network was then constructed using gene expression profiles from GC tissues and visualized as a correlation circle plot (Figure 1A). Differential expression analysis between normal and GC tissues revealed significant dysregulation in nine genes: AZIN2, HNF4A, ARG1, ARG2, ASL, ODC1, OTC, AZIN1, and ALDH18A1 (Figure 1B). Finally, a protein-protein interaction (PPI) network for these nine differentially expressed genes was generated using the STRING database (https://string-db.org) and visualized in Cytoscape (Figure 1C).

Consensus clustering based on the expression levels of differentially expressed ArMGs was performed for the molecular classification of GC samples (Figures 2A–C). The results showed that when K = 3, the boundaries among the three subtypes were relatively distinct. Accordingly, the GC dataset was divided into three clusters: Cluster 1 (C1), Cluster 2 (C2), and Cluster 3 (C3). Subsequent survival analysis revealed statistically significant differences in overall survival among the subtypes (P < 0.05; Figure 2D). Notably, the C1 subtype exhibited a more favorable survival probability than the other two subtypes. The distinct arginine metabolism patterns observed among the C1, C2, and C3 subtypes likely arise from the convergence of upstream oncogenic signaling pathways and tumor microenvironmental pressures.

Differential expression analysis was performed based on the subtype groupings using the limma package to identify DEGs between the subtypes. DEGs were selected according to the criteria of P-value < 0.05 and |logFC| > 0.585. Between the C1 and C2 groups, a total of 1,621 DEGs were identified, including 816 upregulated and 805 downregulated genes (Figure 3A). Between the C1 and C3 groups, 399 DEGs were detected, with 315 upregulated and 84 downregulated genes (Figure 3B). Between the C2 and C3 groups, 1,454 DEGs were identified, comprising 859 upregulated and 595 downregulated genes (Figure 3C). The intersection of these three DEG sets yielded 88 overlapping genes (Figure 3D).

Clinical and survival data for patients with GC were obtained from the TCGA database. Univariate Cox regression analysis identified nine genes that were significantly associated with prognosis (P < 0.05; Figure 4A). Seven key prognostic genes were selected using the Lasso regression feature selection algorithm. The TCGA dataset was randomly divided into training and test sets at a ratio of 4:1. Lasso regression analysis was then applied to derive the model coefficients for each sample (Figures 4B–D). The optimal riskscore value was used for subsequent analysis. The risk score formula was defined as: RiskScore = TMEM171× (−0.152323115295076) + SLC5A1 × (−0.121599685780654) + DEGS2 × (−0.0569621957618535) + MGP × 0.0623411051770642 + C7 × 0.064629072356396 + HMGCS2 × 0.137353168245766 + CREB3L3 × 0.154874633934315.

Samples were stratified into high-risk and low-risk groups based on the median risk score. A nomogram was constructed to visually present the results of regression analysis. Logistic regression analysis revealed that the risk score significantly contributed to the predictive scoring process of the nomogram (Figure 5A). In the nomogram, the points scale assigns a score for each variable. For continuous variables (Age, Risk Score), higher original values correspond to higher points. The cut-off for Age (60 years) was determined based on the median age in our cohort and is commonly used in clinical stratification. Categorical clinical variables (Gender, T, N, M stage, Grade, Stage) are displayed with their actual clinical categories on the axis for clarity. Predictive analysis was performed for 3-year and 5-year overall survival (OS) in patients with GC (Figure 5B), and a decision curve analysis (DCA) was generated to evaluate clinical benefit (Figure 5C). Kaplan-Meier survival analysis demonstrated significantly shorter OS in the high-risk groups across both the training and test cohorts (Figures 5D,E). Receiver operating characteristic (ROC) curves further confirmed the predictive performance of the model (Figures 5F,G).

Processed transcriptomic datasets with complete survival information (GSE26901, GSE62254) were obtained from the GEO database for external validation. The established prognostic model was applied to stratify patients with GC into high-risk and low-risk groups. Kaplan-Meier survival analysis was used to evaluate survival differences between the two groups and assess the stability of the prediction model. The results indicated that OS in the high-risk group was significantly lower than that in the low-risk group within the GEO validation datasets (Figures 6A,B). To further verify the predictive accuracy of the model, ROC curve analysis was conducted using the external datasets. The results demonstrated that the model exhibited strong prognostic performance (Figures 6C,D).

The therapeutic efficacy of surgery combined with chemotherapy for early-stage GC is well established. In this study, drug sensitivity data from the GDSC database were combined with the R package “oncoPredict” to estimate the chemotherapy sensitivity of each tumor sample. The relationship between risk score and sensitivity to common chemotherapeutic agents was further examined. The results indicated that the risk score was significantly associated with sensitivity to several drugs, including Camptothecin_1003, Vinblastine_1004, Cisplatin_1005, Cytarabine_1006, Docetaxel_1007, and Gefitinib_1010 (Figure 7A). In addition, immunotherapy response prediction based on the risk score was performed using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. The TIDE score integrates expression signatures of T-cell dysfunction and exclusion to predict clinical response to immune checkpoint blockade therapies, primarily targeting the PD-1/PD-L1 and CTLA-4 axes. The results revealed that patients in the high-risk group exhibited higher TIDE scores, indicating a poorer predicted response to anti-PD-1/CTLA-4 immunotherapy (Figure 7B).

SNP-related data for processed GC samples were downloaded and analyzed. The 30 most frequently mutated genes were selected to compare mutation patterns between high-risk and low-risk patient subgroups. Mutation frequency differences were visualized using the R package “ComplexHeatmap” (Figure 8A). Notably, the mutation frequencies of TTN and several other genes were significantly lower in the high-risk group than in the low-risk group. Next, associations between the prognostic risk score and established immunotherapy biomarkers were assessed. Significant differences were observed between high-risk and low-risk groups in epithelial–mesenchymal transition (EMT) scores, microsatellite instability (MSI), tumor mutation burden (TMB), and tumor neoantigen (Figures 8B–D).

The tumor microenvironment, consisting of cancer-associated fibroblasts, immune cells, extracellular matrix components, signaling molecules, and malignant cells, plays a crucial role in tumor diagnosis, prognosis, and therapeutic response. To investigate the molecular mechanisms linking risk scores to GC progression, correlations between risk scores and immune infiltration profiles were assessed. Comparative analysis between the high-risk and low-risk subgroups revealed significant differences in immune cell composition (Figure 9A). In the high-risk group, decreased infiltration of dendritic cells activated, macrophages M0, mast cells activated, neutrophils, and natural killer (NK) cells resting was observed. Conversely, higher levels of B cells naive, dendritic cells resting, mast cells resting, and monocytes were detected (Figure 9B).

To further explore the signaling pathways associated with high and low-risk groups, the molecular mechanisms by which risk scores may influence tumor progression were investigated. GSVA revealed that differential pathways between the two groups were primarily enriched in HEDGEHOG_SIGNALING, IL6_JAK_STAT3_SIGNALING, and IL2_STAT5_SIGNALING pathways (Figure 10A). GSEA identified additional pathways, including the cAMP signaling pathway, Rap1 signaling pathway, and Wnt signaling pathway (Figure 10B). The molecular interaction network among these pathways is shown in Figure 10C.

Single-cell transcriptomic analysis was performed using the t-SNE algorithm, resulting in the identification of 24 distinct cellular subtypes (Figure 11A). These subtypes were annotated into 14 major cell lineages: venular endothelial cells, pericytes, fibroblasts, parietal cells, chief cells, dendritic cells, macrophages, monocytes, plasma cells, B cells, mast cells, NKT cells, T cells, and epithelial cells (Figure 11B). The expression patterns of the seven prognostic model genes were visualized across all annotated lineages (Figures 11C,D). In addition, the AUCell function was used to perform quantitative analysis of immune-related and metabolic pathways at the single-cell level, demonstrating correlations between model gene expression and specific immune and metabolic signaling pathways (Figure 11E).

A comprehensive review of the literature on the nine ArMGs was conducted in the context of cancer research. Based on the co-expression network analysis between the model genes and ArMGs, ODC1 and ALDH18A1 were selected as target genes for further investigation.

Validation of ODC1 and ALDH18A1 expression and clinical relevance in GC was performed using cell lines and clinical tissue specimens. qPCR analysis revealed that both genes were significantly upregulated in GC tissues compared with adjacent normal tissues (Figure 12A). Two GC cell lines (AGS and HGC-27), which showed significant differential expression of the target genes relative to the normal gastric epithelial cell line GES-1, were selected for subsequent functional assays. To explore the biological roles of ODC1 and ALDH18A1 in GC progression, siRNAs targeting each gene were transfected into AGS and HGC-27 cells. The knockdown efficiency was confirmed by PCR (Figure 12B). The functional consequences of gene silencing were evaluated using cell proliferation, colony formation, and transwell invasion assays. Experimental results demonstrated that knockdown of ODC1 and ALDH18A1 significantly inhibited GC cell proliferation (Figure 12C), colony formation (Figure 12D), and invasive capacity (Figure 12E).