Identification of prognostic genes associated with sphingosine-1-phosphate in gastric cancer to construct a risk mode

Introduction

Background

According to Global Cancer Statistics 2020, gastric cancer (GC) is one of the five most common malignant tumors and the fourth leading cause of cancer-related deaths worldwide (1). With an aging population, the incidence and mortality of GC are increasing year by year, and the characteristics of high metastasis, high tumor heterogeneity and high chemotherapy resistance are the main reasons for poor prognosis and low survival rate (2,3). Risk factors include Helicobacter pylori infection, age, genetic factors, poor diet and smoking (4,5). In China, more than 80% of GC patients have advanced disease at diagnosis (6). Although surgical resection is an effective treatment for GC and recent advances in chemotherapy have improved progression-free survival and overall survival, the prognosis of patients remains poor, with the 5-year survival rate for patients with advanced GC being approximately 20% (7). The main diagnostic methods for GC are invasive endoscopy and biopsy, which are expensive. Some of the serum biomarkers identified so far have poor specificity. Hence, the search for early specific biomarkers is crucial for the treatment and prognosis of GC.

Rationale and knowledge gap

Sphingolipids are biologically active signaling lipids that play key roles in many cellular processes. Sphingomyelin is converted into ceramide via the action of sphingomyelinase, and ceramide can be further hydrolyzed into sphingosine by ceramidase (8). Sphingosine, in the presence of adenosine triphosphate (ATP) and through the catalysis of sphingosine kinases SPHK1/SPHK2, is phosphorylated to generate sphingosine-1-phosphate (S1P), a process accompanied by the production of adenosine diphosphate (ADP) (9). On one hand, S1P exerts biological effects including “promoting fibroblast and endothelial cell migration, as well as angiogenesis”, thereby inducing an “anti-apoptotic” response. Subsequently, S1P can be dephosphorylated back to sphingosine by sphingosine-1-phosphate phosphatases (SGPP1/2), forming a metabolic cycle. On the other hand, S1P undergoes cleavage by sphingosine-1-phosphate lyase (SGPL1) to generate hexadecenal and phosphoethanolamine, ultimately exerting “pro-apoptotic” effects (Figure 1). In a mouse xenograft tumor model, ceramide metabolism to S1P is metabolized through the function of SPHK1 and SPHK2, which play an important role in chemotherapy and radiotherapy for GC through S1PR-dependent signaling (10). However, to date, there is a lack of systematic studies on S1P-related genes (SRGs) as new prognostic indicators for GC.

Figure 1 Regulatory patterns of S1P metabolism, genes, and phenotypes. ADP, adenosine diphosphate; ATP, adenosine triphosphate; S1P, sphingosine-1-phosphate.

Objective

In this study, we used GC data from The Cancer Genome Atlas (TCGA) database and screened seven genes related to S1P through a series of bioinformatics methods. We constructed a risk model and performed independent prognostic analyses, enrichment analyses, immune microenvironment analyses and drug sensitivity analyses. These provided a theoretical basis for research into the prognosis of GC, associated with SRGs. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-172/rc).

Methods

Data extraction

Gene expression data, clinical data and survival data for patients with GC were retrieved from TCGA database [National Cancer Institute Genomic Data Commons (GDC), https://portal.gdc.cancer.gov/]. The training set, TCGA-GC, contained 375 GC samples and 32 normal samples, of which 350 had complete survival data. The GSE62254 dataset was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and included 300 GC samples with complete survival information and gene expression data. This dataset was used as a validation set (11). A total of 15 SRGs were identified from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/). These were SPHK2, S1PR3, S1PR1, AKT1, S1PR5, EZR, SPNS2, S1PR2, GPR6, MFSD2B, S1PR4, PIK3CB, PIK3CG, RAC1 and SPHK1. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential expression analysis and weighted gene co-expression network analysis (WGCNA)

Differential expression analysis was carried out between the GC and normal samples from the TCGA-GC dataset via the ‘DESeq2’ package (version 1.38.0) to obtain differentially expressed genes (DEGs) (12). Screening conditions were |log₂ fold change (FC)| >1 and P_adj<0.05. Volcano plots were generated using the ggplot2 package (version 3.3.5) (13), while heat maps were created using the pheatmap package (version 1.0.12) (14). The heat maps showed the top 10 up-regulated genes and the top 10 down-regulated genes in the log₂FC ranking. In order to evaluate the differences in SRGs between the GC and normal group, a single-sample gene set enrichment analysis (ssGSEA) score was calculated for the SRGs via the ‘GSVA’ package (version 1.42.0), using the TCGA-GC dataset (15). A WGCNA was then carried out via the ‘WGCNA’ package (version 1.71) (16). All samples from the TCGA-GC dataset clustered together and outliers were eliminated. Next, an optimal soft threshold was determined to ensure the co-expression network adhered to the distribution of a scale-free network. A dynamic tree cutting algorithm was employed to generate distinct modules, each containing a minimum of 100 genes, with a mergeCutHeight of 0.25. The modules with the highest correlation to the ssGSEA score for SRGs were considered to be key modules for subsequent analysis. The genes contained in these modules were key module genes.

Enrichment analysis and protein-protein interaction (PPI) network construction

The candidate genes were derived from the intersection between DEGs and key module genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were also conducted via ‘cluster Profiler’ (version 4.2.2). The significance level was P<0.05 (17). In order to explore the network interactions between candidate genes, a PPI network was established, with an interaction score >0.9, via the search tool for the retrieval of interacting genes (STRING) database (https://string-db.org/).

Risk model construction, evaluation and verification

Candidate genes were incorporated into a univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) algorithm to determine the prognostic genes, using the TCGA-GC dataset. The ‘survival’ (version 3.3-1) program (18) was used to conduct the univariate Cox regression analysis and the screening conditions were hazard ratio (HR) ≠1 and P<0.01. The LASSO algorithm was performed using the ‘glmnet’ package (version 4.1-2) (19). The formula for calculation of the risk score was: Riskscore=∑i=1n(coefficient(genei)∗espr(genei)). Based on the median value of the risk score, patients were separated into low- and high-risk groups. The ‘survminer’ program (version 0.4.9) (20) was used to plot Kaplan-Meier (K-M) curves. The receiver operating characteristic (ROC) curves (1-, 3-, 5-year) were plotted using ‘survivalROC’ (21), to assess the risk model. The GSE62254 dataset was regarded as an external verification set.

Independent prognostic analysis

The clinical information from the TCGA-GC dataset, which included risk score, cancer stage, histological information, tumor (T)-stage, age, metastasis (M)-stage, gender and node (N)-stage, was incorporated into the univariate and multivariate Cox regression analysis (P<0.05) and proportional hazards (PH) hypothesis test (P>0.05). A nomogram was constructed on the basis of independent prognostic factors using the ‘rms’ program (version 6.3-0) (22). A calibration curve was plotted using ‘rms’. A decision curve analysis (DCA) was performed to predict the effect of the nomogram.

Functional enrichment analysis

A gene set enrichment analysis (GSEA) was performed to investigate the high-risk and low-risk groups in the TCGA-GC dataset using the ‘clusterProfiler’ package (version 4.2.2) (17). The screening conditions were |normalized enrichment score (NES)| >1 and P<0.05.

Comprehensive immune analysis

The ESTIMATE algorithm was used to calculate the stromal score, immune score and estimate score in the TCGA-GC dataset. The Wilcoxon test was used to compare the difference in the three scores between the high- and low-risk groups (P<0.05). The content of 22 immune cell types was calculated using the CIBERSORT algorithm. The difference in level of immune cell infiltration between the GC and normal group was calculated by the Wilcoxon test, using the TCGA-GC dataset for the 22 immune cell types.

Construction of a competitive endogenous RNA (ceRNA) network

The TargetScan, miRanda and miRDB databases were used to predict microRNAs (miRNAs) connected with the prognostic genes, via the R package ‘multiMiR’ (version 2.3) (23). The final miRNAs were derived from the intersection between the three databases. The Starbase database was used to predict long non-coding RNAs (lncRNAs) associated with miRNAs. The ‘Cytoscape’ software was used to establish the ‘lncRNA-miRNA-messenger RNA (mRNA)’ network (24).

Drug sensitivity analysis

A total of 138 chemotherapy/targeted therapy drugs were identified from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). In order to explore the predictive ability of these drugs, the half-maximal inhibitory concentration (IC50) of the drugs was calculated using the pRRophetic algorithm. The IC50 was compared between the GC and normal group for each drug.

Cell transfection

Cells were quantified during the logarithmic growth phase to assess population dynamics. Seed 6×10⁵ cells into each well of a six-well plate and incubate overnight in a cell culture incubator to allow for proper cell adhesion. Removed the original culture medium and replaced it with 1.75 mL of fresh DMEM medium that was free of serum and antibiotics. Prepared two 1.5 mL centrifuge tubes and added 120 µL of OPTI-MEM medium to each. To tube A, added 5 µL of Lipofectamine 3000, and to tube B, added 5 µL of siRNA. Incubated both tubes at room temperature for 5 minutes. Subsequently, gently mixed the contents of tubes A and B and incubated the mixture at room temperature for an additional 15 minutes. Transfered the resulting transfection complex to the six-well plate. Placed the plate back into the cell culture incubator and allowed transfection to proceed for 6 hours. After this period, replaced the medium with fresh DMEM supplemented with 10% fetal bovine serum (FBS) and antibiotics. The transfected cells were suitable for routine phenotypic analysis after 48 hours.

Reverse transcription quantitative polymerase chain reaction (RT-qPCR) analysis

Complementary DNA (cDNA) synthesis was performed in PCR tubes using a 20 µL reaction mixture containing 2.0 µg total RNA, 1.0 µL Oligo(dT)₁₈ primer, 4.0 µL 5× reaction buffer, 2.0 µL 10 mM dNTP Mix, 1.0 µL RiboLock™ RNase Inhibitor, and 1.0 µL RevertAid™ Reverse Transcriptase, with nuclease-free water added to a final volume of 20 µL. Reactions were carried out in a thermal cycler under the following conditions: 42 °C for 60 min followed by 70 °C for 5 min.

For quantitative PCR, reactions were performed using PowerUp™ SYBR™ Green Master Mix (Thermo Fisher Scientific, Massachusetts, USA; Cat# A25742) in 20 µL volumes comprising 10 µL master mix, 1 µL each of forward and reverse primers (10 µM), 1 µL cDNA template, and 7 µL nuclease-free water.

Wound healing assay

Horizontal lines were drawn across the back of a 6-well plate using a marker pen, aligned with a ruler. Cells in the logarithmic growth phase were trypsinized, counted, and resuspended to form a cell suspension. According to the experimental groups, cells were seeded into the 6-well plate at a density of 5×10⁵ cells per well. After the cells reached confluence overnight, a sterile 200 µL pipette tip was used to create scratch wounds perpendicular to the horizontal lines on the back of the plate. The dislodged cells were washed away with 1× phosphate-buffered saline (PBS). Following the group designations, serum-free medium containing the respective drug treatments was added to each well, and the cells were returned to the incubator for continued culture. Images of the wound areas were captured under a microscope at 0 and 24 h post-scratching, and relevant data were recorded.

Statistical analysis

Statistical tests were conducted using R software (version 4.2.2). The Wilcoxon test was used to analyze the difference between two groups. The RT-qPCR data were analyzed for normality using the t-test. Quantitative assays were performed in three independent experiments with three technical repeats. Statistical significance was defined as P<0.05. All analyses were conducted using GraphPad Prism 9.0.

Results

Discovery of 260 candidate genes

A total of 4,484 DEGs were identified between the normal and GC groups from the TCGA-GC dataset. Among these DEGs, 2,349 genes were up-regulated, while 2,135 genes were down-regulated in the GC group (Figure 2A,2B). A WGCNA was performed for the TCGA-GC dataset. Sample clustering demonstrated that there were no outlier samples (Figure S1). An optimum soft-threshold power value of β=5 was used to ensure a scale-free network, when scale-free R² reached 0.85 and average connectivity was close to 0 (Figure 2C). Eighteen modules were generated (Figure 2D). We then identified MEgreen (cor =0.56) and MEtan (cor =0.54) modules that were correlated with the ssGSEA score for the SRGs (Figure 2E). A total of 1,398 key module genes were identified from the MEgreen and MEtan modules. Two hundred and sixty candidate genes were derived from the intersection between the DEGs and the key module genes (Figure S2).

Figure 2 Discovery of 260 candidate genes. (A) Volcano map of differentially expressed genes. Grey dots represent genes with no significant expression, red dots represent significant up-regulation, and green dots represent genes with significant down-regulation. (B) Heat map of differentially expressed genes. The top section is a heat map of the density of up- and down-regulated gene expression in the samples, showing lines for the five quartiles and the mean. The darker the colour of the heat map, the higher the gene expression in the sample (red for high expression, blue for low expression). (C) Screening of soft thresholds. The horizontal axes all represent the weight parameter power value, the vertical axis of the left graph scale-free fit index, i.e., signed R², and the vertical axis of the right graph represents the mean value of the neighbour-joining function of all genes in the corresponding gene module. (D) The top half is a hierarchical clustering dendrogram of genes, and the bottom half is gene modules. (E) Heatmap of correlations between modules and traits. In the heatmap in the middle part, the darker the colour the higher the correlation is higher. Red indicates positive correlation and blue indicates negative correlation. FC, fold change; GC, gastric cancer.

Enrichment analysis and establishment of a PPI network

The GO results indicated that the candidate genes were mainly involved in ‘mononuclear cell differentiation’, ‘external side of plasma membrane’ and ‘extracellular matrix (ECM) structural constituent’ (Figure 3A). The KEGG results demonstrated that the candidate genes were mainly involved in the ‘PI3K Akt signaling pathway’ and ‘cytokine cytokine receptor interaction’ (Figure 3B). Moreover, it was found that the MATN3, CXCL9, ASPN and FCGR1B genes were more important than another set of genes in the PPI network (Figure 3C).

Figure 3 Enrichment analysis and establishment of a PPI network. (A) GO enrichment analysis of candidate genes. Horizontal coordinate is the number of genes involved in GO function, vertical coordinate is the name of GO function. The colour change represents the significance p.value change. (B) KEGG enrichment analysis of candidate genes. The horizontal coordinate is the number of genes involved in KEGG function and the vertical coordinate is the name of KEGG function. Colour changes represent significance P value changes. (C) PPI network of candidate genes. Each point represents a gene. ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.

Risk model prediction of the survival status of GC patients

Nineteen genes were identified by the univariate Cox regression analysis. This analysis showed that BATF2 is a protective gene (HR <1) for GC, whilst MATN3, CDH11, CTHRC1, LAMP5, LPPR4, P4HA3, ADAM12, SPARC, SERPINE1, CST2, FNDC1, GPR176, VCAN, OLFML2B, LOX, FCN1, ASPN and NOX4 are risk genes (HR >1) for GC (Figure 4A). The LASSO algorithm was used to screen seven prognostic genes, MATN3, LAMP5, SERPINE1, BATF2, LOX, FCN1 and ASPN (Figure 4B). The risk model was then established in the TCGA-GC dataset (Figure 4C). The risk score was calculated as:

Riskscore=0.164541036×MATN3+0.019860819×LAMP5​+0.143943288×SERPINE1−0.110590768×BATF2​+0.003098167×LOX+0.079905396×FCN1​+0.011991750×ASPN[1]

Figure 4 Risk model prediction of the survival status of GC patients. (A) Univariate Cox forest plot of prognosis-related genes. (B) Above: plot of correlation coefficient versus lambda in LASSO regression. Bottom: trend of cross-validation error in LASSO regression as a function of lambda. (C,D) Risk curve and survival status in training set. Risk curves with red points for high-risk samples and blue points for low-risk samples. (E) K-M survival curve in training set. (F) Total ROC curve in training set. (G) K-M survival curve in validation set. (H) Total ROC curve in validation set. (I,J) Risk curve and survival status in validation set. Risk curves with red points for high-risk samples and blue points for low-risk samples. AUC, area under the curve; CI, confidence interval; GC, gastric cancer; LASSO, least absolute shrinkage and selection operator; K-M, Kaplan-Meier; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

On the basis of the median value of the risk score, patients were separated into two risk groups (Figure 4D). The K-M curves revealed a worse prognosis for high-risk patients in the TCGA-GC dataset (Figure 4E). The area under the curve (AUC) values in the model that forecasts survival at 1-, 3- and 5-years were 0.64, 0.69 and 0.76, respectively, in the TCGA-GC dataset. This indicated that the survival status of GC patients can be forecasted through the risk model (Figure 4F). Survival status was different between the two groups in the GSE62254 dataset (Figure 4G). The AUC values (1-, 3-, and 5-year) exceeded 0.60 in this dataset (Figure 4H), and as the risk score increased, the survival time of patients decreased and the number of deaths increased (Figure 4I,4J).

Prediction efficiency of the nomogram

The T-stage, risk score, age, N-stage, gender, cancer stage and M-stage were found to be associated with prognosis in GC patients, using the univariate Cox analysis (Figure 5A). These clinical features passed the PH hypothesis test. Cancer stage, age and risk score were also screened as independent prognostic factors, using the multivariate Cox analysis (Figure 5B). A graphical model known as a nomogram was developed to incorporate key variables that influence prognosis. The accuracy of this nomogram in forecasting outcomes for GC was validated through calibration plots, which confirmed its high effectiveness in predicting the disease’s progression (as shown in Figure 5C,5D). Moreover, the results of the DCA reflected the prediction effect of the nomogram (Figure 5E). In the TCGA-GC and GSE62254 datasets, MATN3, LAMP5, SERPINE1, LOX, FCN1 and ASPN showed high levels of expression in the high-risk group, while the expression of BATF2 was low (Figure S3).

Figure 5 Prediction efficiency of the nomogram. (A) Forest plot of univariate Cox regression clinical indicators. The lines on the right side indicate their confidence intervals. t represents the size of the tumour and the degree of infiltration, n represents regional lymph node metastasis or not, and m represents whether the tumour metastasized distally or not. (B) Forest plot of multivariate Cox regression clinical indicators. Variables, P value and HR are shown on the left, and the lines on the right indicate their confidence intervals. n represents regional lymph node metastasis or not. (C) Nomogram of predicted 1-, 3-, and 5-year survival probabilities for patients with GC. (D) Calibration curves for nomogram. The horizontal coordinate is the predicted probability, the vertical coordinate is the actual probability, and the closer the slope of the calibration curve is to 1, the more accurate the model prediction is. (E) DCA decision curves for nomogram. CI, confidence interval; DCA, decision curve analysis; GC, gastric cancer; HR, hazard ratio; N, node; T, tumor.

Functional enrichment analysis of the two groups

The results demonstrated that ‘DNA replication’ and ‘ribosome’ were involved in the low-risk group, whilst ‘neuroactive ligand receptor interaction’, ‘ECM receptor interaction’ and ‘focal adhesion’ were involved in the high-risk group (Figure 6). Through the application of GO and KEGG pathway analyses, we have meticulously examined the functional enrichment of two distinct groups, differentiated by their risk profiles. The low-risk group is significantly associated with core cellular functions like ‘DNA replication’, indicative of their role in maintaining basic cellular health. On the other hand, the high-risk group demonstrates enrichment in pathways that involve intricate, and possibly harmful, cellular interactions, such as the ‘neuroactive ligand receptor interaction’ pathway, which often mediates neuron-to-cell communication (Figure 6).

Figure 6 Functional enrichment analysis of high- and low-risk groups. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Immune microenvironment analysis

Stromal score and estimate score were significantly higher in the high-risk group (P<0.05), when compared with the low-risk group (Figure 7A). Moreover, the GC and normal groups showed significant differences in immune infiltration of nine types of immune cells. When compared with the low-risk group, the high-risk group showed higher levels of infiltration of T cell CD4⁺ memory resting, monocyte, B cell plasma and activated mast cells, and lower levels of immune infiltration of M1 macrophage, T follicular helper cell, M0 macrophage and regulatory T cells (Tregs) (Figure 7B).

Figure 7 Immune microenvironment analysis. (A) Violin plots of intergroup differences in stromal score, immunity score and estimate score. The larger the area of the graph, the greater the probability of distribution around a given value. (B) Above: stacking diagram of immune cell infiltration. Different colours represent different immune cells. Bottom: box-and-line diagram of immune cell expression. Horizontal coordinates are immune cells and vertical coordinates are cellular composition. Red is the control group and blue is the disease group. ns, not significant; *, significant; **, very significant; ***, highly significant; ****, extremely significant. GC, gastric cancer.

ceRNA network establishment

An ‘lncRNA-miRNA-mRNA’ network, which contained 89 lncRNAs, 34 miRNAs, MATN3, LAMP5, LOX and ASPN, was constructed using Cytoscape. The miR497HG miRNA regulated LOX by hsa-miR-24-3p (Figure 8A). The hsa-miR-24-3p is a key microRNA that regulates gene expression by inhibiting protein synthesis through binding to mRNA’s 3'UTR, affecting cellular processes like growth and differentiation. Ninety-nine drugs showed significant differences between the two groups. The IC50 values of 39 drugs were significantly higher in the high-risk group, while the IC50 values of 60 drugs were significantly lower in the high-risk group (Figure 8B), when compared with the low-risk group. Six drugs related to GC were selected for analysis. The results demonstrated that the IC50 values of AZD8055 and docetaxel were significantly higher in the low-risk group than in the high-risk group, while the IC50 values of camptothecin, gemcitabine, paclitaxel and sorafenib were significantly higher in the high-risk group, when compared with the low-risk group (Figure 8C).

Figure 8 The ceRNA network establishment. (A) ceRNA regulatory network diagram. Green represents target miRNAs, purple represents prognostic genes, and orange represents lncRNAs. (B) Differential expression plot of 99 drugs between high and low risk groups. Red represents high-risk group, and blue is low-risk group. (C) Differential expression plots of six GC-related drugs between high and low risk groups. Red represents high risk group, blue is low risk group. *, significant; **, very significant; ***, highly significant; ****, extremely significant. ceRNA, competitive endogenous RNA; IC50, half-maximal inhibitory concentration.

Biological function verification of key genes

To further verify the biological functions of the key genes (ASPN, FCN1 and LAMP5) in the above model, we used siRNA transfection to knock down the expression levels of the related genes in GC cells, and then combined the wound healing assay to detect the effect of the expression levels of the key genes on the migration ability of GC cells. The RT-qPCR results showed that the knockdown effect of the corresponding gene was significant after siRNA transfection (Figure 9A-9C). In the wound healing assay, we observed that the migration activity of GC cells significantly decreased as the levels of key genes decreased (Figure 9D-9G).

Figure 9 Validation of tumor-suppressive effects upon knockdown of key genes. (A-C) RT-qPCR analysis confirming the knockdown efficiency of ASPN, FCN1, and LAMP5 in gastric cancer cell line HGC-27; (D-G) Wound healing assays demonstrated that the migration capacity of HGC-27 cells was significantly inhibited after knockdown of ASPN, FCN1, or LAMP5 (10×). *, significant; **, very significant; ***, highly significant; ****, extremely significant. NC, negative control; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

Discussion

GC is the malignant gastrointestinal tumor with the highest morbidity and mortality rate in the world, with a poor 5-year survival rate (1). The incidence rate in East Asian countries is significantly higher, particularly in China, than in European and American countries. The poor prognosis of GC patients affects their quality of life (1). Studies have demonstrated that the S1P signaling pathway promotes chemoresistance in tumor cells by activating the PI3K-Akt or NF-κB signaling pathways, upregulating drug resistance proteins such as ABCG2 (25,26). Concurrently, S1P facilitates the polarization of tumor-associated macrophages toward the immunosuppressive M2 phenotype, suppresses T-cell infiltration, and contributes to immune evasion. However, existing research has primarily focused on individual cancer types (e.g., breast cancer, lung cancer), with limited exploration of the role of S1P in regulating ECM receptor interactions and its prognostic implications in GC-specific models (27). Therefore, investigations into the molecular biology of GC are vital to improve this prognosis. Analysis of the expression profiles of genes that affect the prognosis of GC may allow patient outcomes to be evaluated and predicted.

This study represents the first investigation to incorporate S1P-associated genes into transcriptomic profiling in the context of GC, revealing a significant association between S1P signaling and the ECM receptor interaction pathway. We identified seven prognostic genes and constructed a risk prediction model, demonstrating that core S1P pathway genes (e.g., ASPN) serve as independent prognostic biomarkers. Stratifying GC patients into high- and low-risk groups using this model revealed that the high-risk group exhibited significantly higher risk scores and stromal scores compared to the low-risk group. Pathway enrichment analysis showed that “DNA replication” and “ribosome” pathways dominated in the low-risk group, whereas “neuroactive ligand-receptor interaction”, “ECM receptor interaction”, and “focal adhesion” pathways were enriched in the high-risk group. Notably, we innovatively discovered that S1P may synergistically regulate the tumor immune microenvironment in the high-risk group via integrin αvβ3-mediated ECM remodeling, a mechanism that has not been systematically reported in previous studies. These findings not only bridge critical gaps in understanding the multifaceted regulatory roles of the S1P pathway in GC but also provide a theoretical foundation for targeting the S1P pathway in combination with chemotherapy or immunotherapy.

It has been reported that low BATF2 expression may lead to poor prognosis in hepatocellular carcinoma (HCC) (2). The BATF2 gene is located on chromosome 11q12-11q13 and is a newly discovered tumor suppressor gene (3). In lung adenocarcinoma, BATF2 deficiency leads to invasion and metastasis (4). Patients with BATF2-negative colorectal cancer have worse pathological differentiation, higher invasion and greater susceptibility to metastasis than patients with BATF2-positive colorectal cancer (5). In this study, we found that BATF2 expression had a protective effect on GC patients and patients with a high level of expression of the BATF2 gene had a better prognosis. Therefore, BATF2 is considered to be a prognostic indicator and potential target for gene therapy in various cancers.

The LAMP5 gene also plays a crucial role in a variety of tumors and is a member of the LAMP family (6). A previous study has shown that LAMP5 is specifically expressed in blastic plasmacytoid dendritic cell tumors (7) and in mixed lineage leukemia-rearranged (MLL-r) leukemias. It is associated with poor prognosis. Down-regulation of LAMP5 inhibits NF-κB (28). In GC, knock down of LAMP5 significantly inhibits the proliferation, invasion and metastasis of GC cells, and increases apoptosis and cell cycle arrest (29).

The serine protease inhibitor clade E member 1 (SERPINE1) protein is a serine protease inhibitor that consists of 379 amino acids and is mainly synthesized and secreted by vascular endothelial cells, adipocytes and platelets (30). In recent years, abnormal levels of SERPINE1 have been found in various types of cancer. Overexpression of SERPINE1 has been observed in breast cancer (31), non-small cell lung cancer (32) and bladder cancer (33). Therefore, SERPINE1 is expected to become a promising new target for cancer diagnosis and treatment.

Our results demonstrate that ECM receptor interaction plays a critical role in the high-risk group. A study has shown that in the high-risk group of idiopathic pulmonary fibrosis, hyperactive pulmonary fibroblasts produce excessive ECM components (e.g., collagen) through ECM receptor interaction, leading to scar tissue formation and progressive pulmonary fibrosis with impaired lung function (34). S1P exerts regulatory control over ECM remodeling by activating its receptors (e.g., S1P1 and S1P2), thereby modulating integrin signaling pathways (35). The S1P-S1P1 axis enhances the activation of integrin αvβ3, promoting adhesion of tumor cells to fibronectin (FN) and accelerating the secretion of ECM-degrading enzymes [e.g., matrix metalloproteinase (MMPs)], which disrupts ECM structural integrity (36). In the high-risk group of tumors, bidirectional interactions between tumor cells and the surrounding ECM via receptor engagement facilitate tumor cell proliferation, survival, migration, and invasion. Additionally, ECM remodeling dynamically regulates the tumor microenvironment, influencing angiogenesis and immune cell recruitment (37). S1P exerts systemic regulatory effects on tumor cell behaviors (proliferation, apoptosis, migration, angiogenesis) by binding to G protein-coupled receptors (S1PRs) (38). Collectively, these mechanisms elucidate how activation of the ECM receptor interaction pathway in the high-risk group may be mechanistically linked to S1P-driven remodeling of the tumor microenvironment.

Analysis of two independent mRNA microarrays from the GSE27342 database and the TCGA cohort allowed the identification of LOX-1 up-regulation in approximately 100 GC tissues. In addition, an association was also found between high LOX-1 expression and depth of invasion, tumor, lymph node metastasis (TNM) and overall survival (15). The LOX-1 protein promotes cell migration and invasion, in vitro, and enhances epithelial-mesenchymal transition (EMT) in GC cells by activating the PI3K/AKT/GSK3β pathway. After LOX-1 combines with oxidized low-density lipoprotein (ox-LDL), it can increase several pro-angiogenic factors, such as VEGF, which leads to tumor growth, invasion and metastasis (39). The VEGF-C protein is mainly found in and secreted by cancer tissues and serves as a strong stimulator that promotes lymphatic angiogenesis and lymph node metastasis. Ma et al. found that plasma levels of ox-LDL were elevated in patients with GC, as measured by an enzyme-linked immunosorbent assay, and were associated with high levels of VEGF-C and lymphangiogenesis. This mechanism works through LOX-1 receptor activation of the NF-κB signaling pathway (40). Therefore, LOX-1 has been classified as an independent predictor of poor patient prognosis. Inhibition of LOX-1 mediates ox-LDL activation and it is a potential therapeutic target to prevent and intervene in early lymph node metastasis in patients with GC.

Ficolin (FCN) family proteins are present in various tissues and FCN1 is mainly a cellular molecule associated with monocytes and neutrophils. A previous study has analyzed the correlation between FCN genes and liver cancer, ovarian cancer and lung cancer, to indicate that there is a link between these genes and tumor development (41). In principle, complement-activating proteins contribute to the development of cancer and may also influence progression of the underlying disease. Sun et al. conducted a comprehensive analysis of FCNs and investigated expression differences, gene mutations, diagnostic value, prognostic value and immune cell infiltration. The results revealed that the expression level of FCN in HCC is significantly lower than in normal liver tissue and is associated with immune cell infiltration. They found that the expression of FCN1 is positively correlated with the risk of immune evasion (42).

Bioactive fibroblasts (ASPNs) are crucial for the growth and progression of tumors. The abnormal expression of ASPN can be used as a potential biomarker for the diagnosis of GC. Ding et al. found that ASPN expression promotes CD44 and MMP-2 expression and induces GC metastasis by stimulating the epidermal growth factor receptor and ERK-CD44/MMP-2 pathway (43). Apoptosis of GC cells is blocked by ASPN, and proliferation and migration are promoted by ASPN binding to LEF-1 (44). We used bioinformatics to show that high expression levels of ASPN in GC leads to poor prognosis of patients. Therefore, ASPN can be used as an important prognostic and diagnostic marker for GC.

For patients with GC, especially those with advanced GC, effective treatments are very limited. Tumor progression is largely influenced by a variety of immune cells and other factors. The type and infiltration level of immune cells around the GC microenvironment, as well as immune-related marker genes, are significantly related to tumor progression. Our results showed that in the high-risk GC group, T cell CD4⁺ memory resting, monocytes, B cell plasma and activated mast cell immune infiltration levels were higher than in the low-risk group, while levels of M1 macrophages, T follicular helper cells, M0 macrophages and Tregs were low in the high-risk group. Macrophages are an important component of tumors and are generally divided into two categories. M1 macrophages can influence anti-tumor immunity and kill tumor cells, while M2 macrophages induce immune suppression. Numerous studies have confirmed that immunotherapy combined with small molecule inhibitors is more effective than immunotherapy alone. High-risk genes MATN3 and ASPN, as ECM components, may promote immune suppression by remodeling the tumor microenvironment to facilitate the infiltration of immunosuppressive cells (e.g., plasma cells and activated mast cells), while simultaneously inhibiting the antitumor activity of M1-polarized macrophages (45,46). The protective gene BATF2 may reverse immune escape by enhancing cytotoxic T-cell function through regulation of the interferon-gamma (IFN-γ) signaling pathway (47). Furthermore, high expression of LAMP5 and LOX could induce pro-inflammatory cytokine secretion [e.g., tumor necrosis factor-alpha (TNF-α), interleukin-6 (IL-6)] via activation of the NF-κB signaling pathway, thereby recruiting immunosuppressive cell subsets (e.g., Tregs, myeloid-derived suppressor cells). These mechanisms collectively highlight how ECM remodeling driven by high-risk genes and inflammatory signaling orchestrates immunosuppressive networks in the tumor microenvironment (48). Our drug response analysis showed that the IC50 values of AZD8055 and docetaxel were significantly higher in the low-risk group than in the high-risk group, and the IC50 values of camptothecin, gemcitabine, paclitaxel and sorafenib were significantly higher in the high-risk group than in the low-risk group. These results suggested that LAMP5, SERPINE1, LOX, FCN1 and ASPN may be promising gene targets and may be useful for the clinical treatment of GC.

In this study, we analyzed the expression levels and prognostic value of GC-related genes using bioinformatics to identify potential predictive biomarkers. However, this study has limitations. Prognostic genes were identified solely through bioinformatics without fully validating in vitro/in vivo mechanistic links between S1P signaling and immune microenvironment/ECM remodeling; future work will verify functions of core genes (MATN3, ASPN) in GC using animal models. Additionally, clinical correlations relied on public datasets lacking chemotherapy/immunotherapy response records and multidimensional features (e.g., histopathology, tumor size), potentially compromising model generalizability. To address these, we will: (I) establish independent cohorts with treatment-response data via clinical collaborations to optimize therapeutic outcome prediction; (II) employ animal models to dissect S1P-integrin regulation of ECM-immune crosstalk; and (III) develop small-molecule inhibitors targeting SRGs for combinatorial immunotherapy, advancing precision medicine in GC.

Conclusions

In summary, this study showed that LAMP5, SERPINE1, LOX, FCN1 and ASPN may represent potential prognostic biomarkers for GC. These genes may play oncogenic roles in gastric tumors, while BATF2 may act as a tumor suppressor. Further studies are needed to explore the functional roles of these genes, particularly in the development of metastasis and cancer progression, to guide clinical direction.

Acknowledgments

We would like to express our sincere gratitude to all individuals and organizations who supported and assisted us throughout this research. Without their support, this research would not have been possible.

Footnote

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

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-172/prf

Funding: This study was supported by the medical and health research project of Zhejiang province (No. 2021KY111) and Zhejiang province scientific research fund of traditional Chinese medicine (No. 2021ZB032).

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

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

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

References

1. Zhang S, Rong P, Chen Q, et al. Suppressor of activator protein-1 regulated by interferon expression in prostate cancer tissues and cells. Life Sci 2019;232:116626. [Crossref] [PubMed]
2. Ma H, Liang X, Chen Y, et al. Decreased expression of BATF2 is associated with a poor prognosis in hepatocellular carcinoma. Int J Cancer 2011;128:771-7. [Crossref] [PubMed]
3. Wang C, Su Y, Zhang L, et al. The function of SARI in modulating epithelial-mesenchymal transition and lung adenocarcinoma metastasis. PLoS One 2012;7:e38046. [Crossref] [PubMed]
4. Liu ZB, Yang Y, Ye XG, et al. Expression and significance of SARI and CCN1 in human colorectal carcinomas. Zhonghua Yi Xue Za Zhi 2011;91:2397-401.
5. Defays A, David A, de Gassart A, et al. BAD-LAMP is a novel biomarker of nonactivated human plasmacytoid dendritic cells. Blood 2011;118:609-17. [Crossref] [PubMed]
6. Beird HC, Khan M, Wang F, et al. Features of non-activation dendritic state and immune deficiency in blastic plasmacytoid dendritic cell neoplasm (BPDCN). Blood Cancer J 2019;9:99. [Crossref] [PubMed]
7. Gracia-Maldonado G, Clark J, Burwinkel M, et al. LAMP-5 is an essential inflammatory-signaling regulator and novel immunotherapy target for mixed lineage leukemia-rearranged acute leukemia. Haematologica 2022;107:803-15. [Crossref] [PubMed]
8. Alkafaas SS, Abdallah AM, Hassan MH, et al. Molecular docking as a tool for the discovery of novel insight about the role of acid sphingomyelinase inhibitors in SARS- CoV-2 infectivity. BMC Public Health 2024;24:395. [Crossref] [PubMed]
9. Alkafaas SS, Elsalahaty MI, Ismail DF, et al. The emerging roles of sphingosine 1-phosphate and SphK1 in cancer resistance: a promising therapeutic target. Cancer Cell Int 2024;24:89. [Crossref] [PubMed]
10. Ogretmen B. Sphingolipid metabolism in cancer signalling and therapy. Nat Rev Cancer 2018;18:33-50. [Crossref] [PubMed]
11. Cristescu R, Lee J, Nebozhyn M, et al. Molecular analysis of gastric cancer identifies subtypes associated with distinct clinical outcomes. Nat Med 2015;21:449-56. [Crossref] [PubMed]
12. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
13. Ito K, Murphy D. Application of ggplot2 to Pharmacometric Graphics. CPT Pharmacometrics Syst Pharmacol 2013;2:e79. [Crossref] [PubMed]
14. Zheng X, Ma Y, Bai Y, et al. Identification and validation of immunotherapy for four novel clusters of colorectal cancer based on the tumor microenvironment. Front Immunol 2022;13:984480. [Crossref] [PubMed]
15. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
16. Nomiri S, Karami H, Baradaran B, et al. Exploiting systems biology to investigate the gene modules and drugs in ovarian cancer: A hypothesis based on the weighted gene co-expression network analysis. Biomed Pharmacother 2022;146:112537. [Crossref] [PubMed]
17. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
18. Lee JH, Jung S, Park WS, et al. Prognostic nomogram of hypoxia-related genes predicting overall survival of colorectal cancer-Analysis of TCGA database. Sci Rep 2019;9:1803. [Crossref] [PubMed]
19. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22.
20. Liu TT, Li R, Huo C, et al. Identification of CDK2-Related Immune Forecast Model and ceRNA in Lung Adenocarcinoma, a Pan-Cancer Analysis. Front Cell Dev Biol 2021;9:682002. [Crossref] [PubMed]
21. Weng W, Chen X, Gong S, et al. Preoperative neutrophil-lymphocyte ratio correlated with glioma grading and glioblastoma survival. Neurol Res 2018;40:917-22. [Crossref] [PubMed]
22. Wu J, Zhang H, Li L, et al. A nomogram for predicting overall survival in patients with low-grade endometrial stromal sarcoma: A population-based analysis. Cancer Commun (Lond) 2020;40:301-12. [Crossref] [PubMed]
23. Ru Y, Kechris KJ, Tabakoff B, et al. The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res 2014;42:e133. [Crossref] [PubMed]
24. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
25. Campos LS, Rodriguez YI, Leopoldino AM, et al. Filamin A Expression Negatively Regulates Sphingosine-1-Phosphate-Induced NF-κB Activation in Melanoma Cells by Inhibition of Akt Signaling. Mol Cell Biol 2016;36:320-9. [Crossref] [PubMed]
26. Li M, Liang RF, Wang X, et al. BKM120 sensitizes C6 glioma cells to temozolomide via suppression of the PI3K/Akt/NF-κB/MGMT signaling pathway. Oncol Lett 2017;14:6597-603. [Crossref] [PubMed]
27. Zhu Q, Zhang K, Cao Y, et al. Adipose stem cell exosomes, stimulated by pro-inflammatory factors, enhance immune evasion in triple-negative breast cancer by modulating the HDAC6/STAT3/PD-L1 pathway through the transporter UCHL1. Cancer Cell Int 2024;24:385. [Crossref] [PubMed]
28. Umeda S, Kanda M, Shimizu D, et al. Lysosomal-associated membrane protein family member 5 promotes the metastatic potential of gastric cancer cells. Gastric Cancer 2022;25:558-72. [Crossref] [PubMed]
29. Siegel RL, Miller KD, Wagle NS, et al. Cancer statistics, 2023. CA Cancer J Clin 2023;73:17-48. [Crossref] [PubMed]
30. Jevrić M, Matić IZ, Krivokuća A, et al. Association of uPA and PAI-1 tumor levels and 4G/5G variants of PAI-1 gene with disease outcome in luminal HER2-negative node-negative breast cancer patients treated with adjuvant endocrine therapy. BMC Cancer 2019;19:71. [Crossref] [PubMed]
31. Sotiropoulos GP, Kotopouli M, Karampela I, et al. Circulating plasminogen activator inhibitor-1 activity: a biomarker for resectable non-small cell lung cancer? J BUON 2019;24:943-54.
32. Becker M, Szarvas T, Wittschier M, et al. Prognostic impact of plasminogen activator inhibitor type 1 expression in bladder cancer. Cancer 2010;116:4502-12. [Crossref] [PubMed]
33. Li C, Zhang J, Wu H, et al. Lectin-like oxidized low-density lipoprotein receptor-1 facilitates metastasis of gastric cancer through driving epithelial-mesenchymal transition and PI3K/Akt/GSK3β activation. Sci Rep 2017;7:45275. [Crossref] [PubMed]
34. Felisbino MB, Rubino M, Travers JG, et al. Substrate stiffness modulates cardiac fibroblast activation, senescence, and proinflammatory secretory phenotype. Am J Physiol Heart Circ Physiol 2024;326:H61-73. [Crossref] [PubMed]
35. Huang K, Liu W, Lan T, et al. Berberine reduces fibronectin expression by suppressing the S1P-S1P2 receptor pathway in experimental diabetic nephropathy models. PLoS One 2012;7:e43874. [Crossref] [PubMed]
36. Kalhori V, Törnquist K. MMP2 and MMP9 participate in S1P-induced invasion of follicular ML-1 thyroid cancer cells. Mol Cell Endocrinol 2015;404:113-22. [Crossref] [PubMed]
37. Gupta R. Epigenetic regulation and targeting of ECM for cancer therapy. Am J Physiol Cell Physiol 2022;322:C762-8. [Crossref] [PubMed]
38. Patmanathan SN, Wang W, Yap LF, et al. Mechanisms of sphingosine 1-phosphate receptor signalling in cancer. Cell Signal 2017;34:66-75. [Crossref] [PubMed]
39. Kapoor P, Deshmukh R. VEGF: A critical driver for angiogenesis and subsequent tumor growth: An IHC study. J Oral Maxillofac Pathol 2012;16:330-7. [Crossref] [PubMed]
40. Ma C, Xie J, Luo C, et al. OxLDL promotes lymphangiogenesis and lymphatic metastasis in gastric cancer by upregulating VEGF C expression and secretion. Int J Oncol 2019;54:572-84. [Crossref] [PubMed]
41. Jang H, Jun Y, Kim S, et al. FCN3 functions as a tumor suppressor of lung adenocarcinoma through induction of endoplasmic reticulum stress. Cell Death Dis 2021;12:407. [Crossref] [PubMed]
42. Sun L, Yu S, Dong C, et al. Comprehensive Analysis of Prognostic Value and Immune Infiltration of Ficolin Family Members in Hepatocellular Carcinoma. Front Genet 2022;13:913398. [Crossref] [PubMed]
43. Ding Q, Zhang M, Liu C. Asporin participates in gastric cancer cell growth and migration by influencing EGF receptor signaling. Oncol Rep 2015;33:1783-90. [Crossref] [PubMed]
44. Zhang Z, Min L, Li H, et al. Asporin represses gastric cancer apoptosis via activating LEF1-mediated gene transcription independent of β-catenin. Oncogene 2021;40:4552-66. [Crossref] [PubMed]
45. Huang Y, Xu X, Lu Y, et al. The phase separation of extracellular matrix protein matrilin-3 from cancer-associated fibroblasts contributes to gastric cancer invasion. FASEB J 2024;38:e23406. [Crossref] [PubMed]
46. Tanaka M. Crosstalk of tumor stromal cells orchestrates invasion and spreading of gastric cancer. Pathol Int 2022;72:219-33. [Crossref] [PubMed]
47. Zong Y, Chang Y, Huang K, et al. The role of BATF2 deficiency in immune microenvironment rearrangement in cervical cancer - New biomarker benefiting from combination of radiotherapy and immunotherapy. Int Immunopharmacol 2024;126:111199. [Crossref] [PubMed]
48. Duprat F, Robles C, Castillo MP, et al. LOX-1 Activation by oxLDL Induces AR and AR-V7 Expression via NF-κB and STAT3 Signaling Pathways Reducing Enzalutamide Cytotoxic Effects. Int J Mol Sci 2023;24:5082. [Crossref] [PubMed]