The TCGA-STAD dataset was downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) (access on 2024) as a transcriptome dataset for gastric cancer (GC). In this dataset, primary tumor samples from GC patients at T1 and T2 stages were extracted, with samples lacking T1/T2 staging information excluded. Based on the lymph node (LN) metastasis status of these samples (NX samples were excluded), the GC patient samples were divided into LN-positive samples (N1, N2, N3) as the NP group and LN-negative samples (N0) as the NO group, consisting of 49 and 40 samples, respectively (Samples with follow-up periods of less than 30 days during clinical evaluation and those lacking survival information were excluded). The resulting dataset was named TCGA-GC-LN. Simultaneously, somatic mutation data, tumor mutation burden (TMB), and clinical data were downloaded. The GSE84433 (GPL6947) dataset was downloaded from the GEO as a validation cohort. A total of 28 tumor samples from GC patients in T1/T2 stages with LN-positive metastasis were selected as the NP group, while 18 tumor samples from GC patients in T1/T2 stages with LN-negative metastasis were selected as the NO group. Additionally, 254 calcium signaling pathway-related genes (CSP-RGs) (Supplementary Table S1) were obtained from KEGG database.

DEGs between the NP and NO groups in the TCGA-GC-LN dataset were identified using the DESeq2 package (v 1.38.0) (Huang A. et al., 2024), using thresholds of P < 0.05 and |log₂ FC| > 0.5. A volcano plot of these DEGs was generated using the ggplot2 package (v 3.4.4) (Gustavsson et al., 2022), with annotations for the top 10 upregulated and downregulated genes. A heatmap of the top 10 upregulated and downregulated genes was generated using the ggplot2 package (v 3.4.4).

DEGs and CSP-RGs were integrated using the ggvenn package (v 0.1.9) (Zhou et al., 2025) to identify the candidate genes. Subsequently, functional enrichment analysis of the candidate genes was carried out with the clusterProfiler package (v 4.7.1.003) (Chen and Boutros, 2011), including Gene Ontology (GO) (P < 0.05) and KEGG enrichments (P < 0.05). The top 5 terms for each category and the top 11 KEGG pathways were displayed based on P-value. Additionally, a protein-protein interaction (PPI) network was established based on candidate genes using the STRING database (confidence score ≥0.4).

Univariate Cox regression analysis was performed on the candidate genes in the TCGA-GC-LN dataset using the survival package (v 3.5–3) (Ni et al., 2023) to identify genes linked to overall survival (OS) (HR ≠ 1, P < 0.2). The results were visualized using a forest plot created with the forestplot package (v 3.1.1). The proportional hazards (PH) assumption test (P > 0.05) was then applied to the univariate Cox analysis results, and genes passing the PH assumption were selected as candidate biomarkers. Additionally, in both the TCGA-GC-LN and GSE84433 datasets, the “survv_cutpoint” function from the “survminer” R package was employed to ascertain the optimal cut-off value for each gene expression level using the maximum selection rank statistical method. Subsequently, based on the determined optimal cut-off value for each candidate biomarker expression, the samples were categorized into high expression and low expression groups. Survival analysis, conducted using the survival package (v 3.5–3), compared survival differences between 2 groups (P < 0.05), and Kaplan-Meier (KM) survival curves were plotted using the survminer package (v 0.4.9) (Shi et al., 2023). Genes that exhibited significant survival differences between the 2 groups and consistent survival trends were selected as biomarkers.

Besides, In the TCGA-GC-LN dataset (excluding samples with missing clinical information), the differences (P < 0.05) in biomarkers expression between different clinical subgroups were compared using the Wilcoxon test.

The reference gene set used was the “c2.cp.kegg.v7.4.symbols.gmt” file from MSigDB. Differential analysis of all TCGA-GC-LN dataset samples was carried out with the DESeq2 package (v 1.38.0) to detect differences between the NP and NO groups. The log₂FC was calculated, and the genes were sorted by log₂FC from largest to smallest. GSEA was performed using the clusterProfiler package (v 4.7.1.003) to identify enriched pathways (|NES| > 1, FDR <0.25, and P < 0.05). The top 5 most significant pathways were visualized based on P-value.

The oncogenic pathways were obtained from the literature (Li et al., 2024) (Supplementary Table S2), and the ssGSEA algorithm was used to calculate the oncogenic pathway scores for both the NP group and the NO group across all samples of the TCGA-GC-LN dataset. The differences in scores between the 2 groups were compared using the Wilcoxon test (P < 0.05), and the results with significant differences were presented.

The ssGSEA algorithm was used to evaluate the proportions of 28 immune cell types in the TCGA-GC-LN dataset. The Wilcoxon test (P < 0.05) was applied to identify differential scores between the NP and NO groups. Spearman correlation analysis was then performed on all TCGA-GC-LN samples to explore the relationships between biomarkers and differential immune cell types, as well as among the differential immune cell types themselves, using the stats package (v 4.2.2) (Aslam et al., 2023) (|cor| > 0.30, P < 0.05) (v 2.2.9). The TIDE platform was used to calculate the TIDE score for patients in both the NP and NO groups, with differences compared using the Wilcoxon test (P < 0.05). The ESTIMATE package (v 1.0.13) (Bi et al., 2020) was then employed to calculate immune, stromal, and ESTIMATE scores for both groups, and differences between them were assessed.

Immune checkpoint genes were obtained from the literature (Supplementary Table S3) (Liu et al., 2023), and the differences in their expression between the NP group and the NO group were compared using the Wilcoxon test (P < 0.05).

In the TCGA-GC-LN dataset (Samples lacking somatic mutation information were excluded), somatic mutation data for each patient were obtained, and the Maftools package (v 2.14.0) (Xu et al., 2021) was used to generate a waterfall plot displaying the top 20 most frequently mutated genes. Next, mutation data from the 2 groups were used for TMB analysis, and TMB score was calculated using the maftools package (v 2.14.0). The TMB scores between the 2 groups were compared using the Wilcoxon test (P < 0.05).

A total of 138 chemotherapeutic drugs were sourced from the GDSC database, and the pRRophetic package (v 0.5) (Yan et al., 2022) was utilized to calculate the IC50 values for patients in the TCGA-GC-LN dataset. Subsequently, the Wilcoxon test (P < 0.05) was employed to explore differences in IC50 values between the 2 groups. The top 3 drugs were displayed, sorted by P-value.

The expression of biomarkers in different cell types of GC was analyzed using the TISCH database. Additionally, the HR associated with biomarkers mutations in different cancer types were presented.

Transcription factors (TFs) associated with biomarkers were predicted through the mirnet database. The TF-mRNA network was visualized. MicroRNAs (miRNAs) related to biomarkers were predicted using the microcosm database. The StarBase database was then used to predict long non-coding RNAs (lncRNAs) targeting the identified miRNAs. The lncRNA-miRNA-mRNA network was visualized using ggalluvial package (v 3.9.1) (Shi et al., 2023).

To assess RET mRNA expression levels, gastric cancer samples at T1 and T2 stages, both with (N = 21) and without (N = 21) lymph node metastasis, were collected from patients who underwent curative surgical resection with lymph node dissection in 2024 at Tianjin Medical University Cancer Hospital, Tianjin, China. None of the patients received neoadjuvant therapy prior to gastrectomy. The Institutional Research Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (approval no. bc2022220) sanctioned all experimental procedures involving these samples. Informed consent was obtained for the collection of all samples in accordance with the Declaration of Helsinki.

Total RNA was extracted using RNAiso Plus (TaKaRa, Japan), and complementary DNA (cDNA) was synthesized with the TaKaRa Reverse Transcription Kit (TaKaRa, Japan). Quantitative real-time PCR was conducted on the QuantStudio 5 system (Applied Biosystems, Foster City, CA, USA) to quantify mRNA levels, employing TB Green Premix Ex TaqTM II (TaKaRa, Japan). β-actin served as the reference gene for data normalization. The primer sequences were as follows: RET:5′- GTC​CTC​TTG​CTC​CAC​TTC​AAC​G/CCT​GGC​AGT​TTT​CCA​CAC​AGA​C -3’; β-actin: 5′-ATA​GCA​CAG​CCT​GGA​TAG​CAA​CGT​AC/CAC​CTT​CTA​CAA​TGA​GCT​GCG​TGT​G-3′

The human gastric cancer cell line AGS was procured from the American Type Culture Collection (Virginia, USA) and cultured in F12K medium supplemented with 10% fetal bovine serum (FBS). The HGC-27 cell line was obtained from the National Infrastructure of Cell Line Resource (Beijing, China) and maintained in RPMI-1640 medium, also supplemented with 10% FBS. Both cell lines were incubated at 37 °C in a humidified atmosphere containing 5% CO₂. All cell lines utilized in this study underwent short tandem repeat (STR) validation to confirm their authenticity.

Initially, the cells were seeded into six-well plates. Subsequently, cell transfection was conducted in accordance with the protocol provided by Lipofectamine™ 3000 (Invitrogen, New York, United States). Following an incubation period of 6 h, a complete medium was added to the cells. The empty pSIH1-puro vector served as the corresponding negative control. The shRNA target sequences were as follows: AGG​GTC​GGA​TTC​CAG​TTA​AAT (shRET#1) and ACC​GCT​GGT​GGA​CTG​TAA​TAA (shRET#2).

The arrangement of colonies was employed to assess the proliferative capacity of the cells. A total of 2,000 cells were seeded into a 6-well plate and incubated at 37 °C with 5% CO₂. After a period of 2 weeks, the cells were fixed and stained using 4% paraformaldehyde and crystal violet, respectively.

For the migration assay, a serum-free cell suspension was introduced into the upper chambers of 24-well Transwell plates (8-μm pore size, Corning, NY, USA) at a density of 4 × 10^4 cells per well. The lower chambers were supplemented with culture medium containing 20% FBS. The invasion assay was conducted using chambers pre-coated with Matrigel (BD Biosciences, USA). Following a 24-h incubation period, the inserts were fixed with 4% paraformaldehyde for 30 min and subsequently stained with 0.4% crystal violet for another 30 min. Cell numbers were quantified by averaging counts from three randomly selected fields.

R (v 4.2.2) was utilized to conduct the bioinformatics analysis. Differences between the 2 groups were assessed using the Wilcoxon test (P < 0.05). Each experiment was independently performed three times. All experimental data were presented as the mean values ±SD. Data analysis was performed using Prism 8.0 software (GraphPad Software), in which *, **and *** indicated P < 0.05, P < 0.01, P < 0.001, respectively.

A total of 2,154 DEGs were identified in the TCGA-GC-LN dataset, including 1,334 upregulated genes and 820 downregulated genes in the NP group (Figures 1A,B). By intersecting DEGs and CSP-RGs, 51 candidate genes were selected (Figure 1C). Enrichment analysis revealed that the candidate genes were strongly linked to 1,068 GO terms (P < 0.05) and 91 KEGG pathways (P < 0.05). For instance, GO terms included “calcium ion transport” and “calmodulin binding” (Figure 1D; Supplementary Table S4). The KEGG pathways primarily included “MAPK signaling pathway” and “cAMP signaling pathway” (Figure 1E; Supplementary Table S5). In the PPI network, a total of 51 proteins were found to interact with each other (such as CALML5 and ADCY1) (Figure 1F).

A total of 14 candidate biomarkers were identified through univariate Cox regression analysis (Figure 2A) and the PH assumption test. In both the TCGA-GC-LN and GSE84433 datasets, the high expression level of RET was associated with poor survival rates and significant differences. Therefore, RET was defined as biomarkers (Figures 2B,C). In addition, RET showed significant differences between the T1 and T2 stages (Figure 2D).

GSEA analysis revealed 36 significantly different pathways between the NP and NO groups (Figure 3A; Supplementary Table S6). Additionally, 2 oncogenic pathways, the “EMT” pathway and the “CSCs activity” pathway, showed significant differences between the NP and NO groups (P < 0.05) (Figure 3B). These results suggested that signaling pathways played a crucial role in tumor progression, particularly oncogenic pathways.

The Heatmap demonstrated the percentages of 28 immune cell types in NP and NO groups (Figure 4A). The scores of 14 immune cell types showed significant differences, including mast cells (P < 0.05) (Figure 4B). The majority of the differential immune cell types exhibited positive correlations, with the strongest correlation observed between monocytes and MDSCs (cor = 0.90, P < 0.001) (Figure 4C; Supplementary Table S7). In addition, Spearman correlation analysis revealed that RET exhibited the strongest correlation with mast cells (cor = 0.55, P < 0.05) (Figure 4D).

The TIDE score did not show a significant difference between the NP and NO groups (Figure 4E). However, the stromal score (P = 0.0051), immune score (P = 0.009), and ESTIMATE Score (P = 0.0024) were significantly elevated in the NP group (Figure 4F). Additionally, the expression of 13 immune checkpoint genes (BTLA, CD209, CD27, CD28, CD40, CD86, HAVCR2, HLA-DOA, HLA-DPB1, HLA-DQA1, THFRSF4, TNFSF18, and TNFSF4) was significantly higher in the NP group (P < 0.05) (Figure 4G). These findings indicated that the tumor microenvironment significantly influenced tumor cells, emphasizing the need to give greater consideration to the tumor microenvironment when developing cancer treatment strategies.

TP53 was the gene with the highest mutation rate in both the NP and NO groups, with mutations including missense mutations (Figure 5A). While no significant difference in TMB score was observed, the TMB score in the NP group was slightly elevated, indicating a higher frequency of gene mutations in this group (Figure 5B).

In total, 15 drugs demonstrated significant differences between the NP and NO groups. AKT inhibitor VIII and BIBW2992 had higher IC50 values in the NP group, whereas NU.7441 showed a higher IC50 value in the NO group (Figure 5C). These results suggested that these drugs might have played an important role in patients with GC LN metastasis.

The expression of RET in GC single cells was analyzed using the TISCH database. The results indicated that RET was expressed in multiple cell types, including DC, mast cells, fibroblasts, myofibroblasts, and gland mucous (Figure 6A). In addition, RET mutations increase the risk of cancers such as KIRC, while decreasing the risk of cancers such as UVM (Figure 6B). Fibroblasts were one of the major sources of the extracellular matrix (ECM). During the early stages of gastric cancer development, fibroblasts secreted a large amount of ECM components, such as collagen. Abnormal accumulation and remodeling of the ECM disrupted the normal growth, differentiation, and signaling of gastric epithelial cells, thus creating conditions conducive to the development of early GC.

A total of 10 TFs were found to be associated with RET, including TP53 (Figure 7A). Additionally, 5 miRNAs and 6 lncRNAs were predicted. The identified regulatory networks included MALAT1-hsa-miR-185-5p-RET and AC007952.4--hsa-miR-185-5p-RET (Figure 7B). These potential regulatory mechanisms had implications for the development of GC LN metastasis.

In both the TCGA-GC-LN and GSE84433 datasets, RET was significantly upregulated in the NP group (Figures 8A,B). The RT-qPCR results of clinical samples also indicated a significant upregulation of RET in the NP group (Figure 8C), consistent with the bioinformatics analysis, further highlighting the importance of RET.

To examine the role of RET in GC cells, we employed shRNAs to deplete endogenous RET in HGC-27 and AGS cell lines. The efficiency of RET knockdown was validated through RT-PCR analysis (Figure 9A). RET-depleted GC cells exhibited a significantly reduced growth rate compared to control cells (Figure 9B). Furthermore, RET knockdown in GC cells markedly impaired cell migration and invasion (Figure 9C). Taken together, these findings indicate that RET plays a crucial role in promoting the progression of GC cells.