RNA-sequencing (RNA seq) datasets, metadata and clinical information for 375 STAD samples and 32 control samples were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The RNA-seq dataset and metadata files were processed using the Perl script (version 1.7.8) (https://bioperl.org/) to obtain the gene expression matrix, followed by ID transformation to obtain the corresponding gene names. STAD samples with indeterminate clinical stages and survival time <30 days were excluded.

Human senescence genes were obtained from the Human Ageing Genomic Resources database (https://genomics.senescence.info/) (Table SI). To create the senescent gene expression matrix file, the limma R package (version 3.56.0) (https://bioconductor.org/packages/release/bioc/html/limma.html) was used to determine the intersection of senescent genes with the gene expression matrix. Pearson correlation analysis was applied to analyze the co-expression of lncRNAs and senescence-related genes. The expression of lncRNAs and senescence-related genes was defined as function variables, and the correlation between them was tested to obtain correlation coefficients and P-values. The SenRLs were extracted with |Correlation Coefficient|>0.4 as well as P<0.001 (26). Cytoscape (version 3.6.1) (https://cytoscape.org/) was used to visualize the coexpression network.

The sample name, survival time and survival status were the only information kept in the clinical data of STAD cases that were obtained from the TCGA database. The survival times were sorted in ascending order to remove samples with unclear survival times. Next, survival status was sorted in ascending order, and samples with unclear survival status were deleted. The processed survival information was merged with the SenRLs expression data to produce SenRLs survival data for STAD for subsequent analysis.

Using R software packages (version 4.3.1; RStudio, Inc.), patients within the SenRLs survival data file were split into training and test sets at random. Univariate Cox regression was performed to screen STAD prognostic SenRLs based on P<0.05. Kaplan-Meier analyses were also conducted requiring KM <0.05 and plotted forest plots using the survival (version 3.4–0; http://cran.r-project.org/package=survival) and survminer (version 0.4.9; http://cran.r-project.org/package=survminer) R packages. Next, these SenRLs were further included in a multivariate Cox regression to get the key lncRNAs applied to create the prognostic risk model and to get the regression coefficients. Risk scores were calculated utilizing the regression coefficients and expression values of SenRLs. The median risk score was used to divide patients into different risk groups. The Sankey diagrams were made using the R tool ggalluvial (version 0.12.5) (https://cran.r-project.org/package=ggalluvial). Survival analysis was then carried out using the survminer (version 0.4.9; http://cran.r-project.org/package=survminer) R package. To examine the survival disparities between the two groups, Kaplan-Meier survival curves were charted. Receiver operating characteristic (ROC) curves were produced using timeROC (version 0.4; http://cran.r-project.org/package=timeROC), and area under the curve (AUC) values were calculated to compare the role of different factors in predicting outcomes. Scatterplot3D (version 0.3–44; http://cran.r-project.org/package=scatterplot3d) R package was applied to create principal component analysis (PCA) plots. Univariate and multivariate Cox regression analyses were carried out on risk value as well as clinicopathological information to clarify whether risk values were independent prognostic indicators.

With the rms (version 6.7–1; http://cran.r-project.org/package=rms) R package, a clinically useful nomogram was developed and validated. Calibration curves were applied to demonstrate the predictive power of the nomogram. To further investigate the biological functions and pathways associated with the prognostic SenRLs signature, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the GSEA software (version 4.1.0; http://www.gsea-msigdb.org/gsea/index.jsp). The analysis was conducted with the following parameters: the phenotypic file was set to h_vs._1, the permutation type was adjusted to phenotype, the number of simulations was fixed to 1,000 and the plot graphs for the top sets of each phenotype were set to 5. Significantly enriched pathways were identified based on a nominal P-value <0.05 and a false discovery rate (FDR) <0.25.

STAD tissue specimens and paired paraneoplastic tissue specimens were acquired from 15 patients treated at Xingtai People's Hospital (Xingtai, China) and identified as STAD by postoperative pathology between March 2025 and September 2025. All included STAD cases were exclusively identified as the intestinal type according to Lauren's classification (1965) (27), with pathological diagnoses established in accordance with the WHO Classification of Tumours of the Digestive System (28). To ensure diagnostic accuracy, the histological subtyping and final diagnosis were independently reviewed and confirmed by two senior pathologists from Xingtai People's Hospital. The patients with STAD included 8 males and 7 females, aged 41 to 68 years. All patients were strictly screened, with no history of other malignant tumors and no prior radiotherapy, chemotherapy, targeted therapy or immunotherapy. The Xingtai People's Hospital Ethics Committee approved the study (approval no. 2025-048). Written informed consent was voluntarily obtained from each patient or their authorized family member. The present study was conducted in accordance with The Declaration of Helsinki. Tissue samples were collected immediately after the operation, quickly frozen with liquid nitrogen and then carefully transported and kept in a refrigerator at −80°C.

Total RNA was extracted using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. First-strand cDNA was synthesized using the Takara PrimeScript RT Kit (Takara Bio, Inc.) following the manufacturer's instructions. RT-qPCR detection employed the SYBR Premix Ex Taq^™ Kit (Takara Bio, Inc.), strictly adhering to the manufacturer's operating procedures. The primers used in this procedure were synthesized by Sangon Biotechnology Co., Ltd., and Table I lists the primer sequences. The following are the thermal cycling conditions employed for qPCR: An initial denaturation step at 95°C for 5 min; followed by 40 cycles, each consisting of 5 sec denaturation at 95°C, 30 sec annealing at 60°C and 30 sec extension at 72°C; concluding with a final extension step at 72°C for 2 min. The relative amounts of genes were analyzed by 2^−∆∆Cq (29) with GAPDH and ACTB as internal reference. A paired t-test was employed to compare the expression levels of these SenRLs between STAD tissues and paired paracancerous tissues.

The immune infiltration landscape of each STAD sample was quantified utilizing CIBERSORT (version 1.03) (https://cibersort.stanford.edu/). For immune cell analysis and visualization using bubble plots the following were utilized: XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS and CIBERSORT on TIMER2.0 (http://timer.cistrome.org/). Following this, the limma, ggplot2 (version 3.4.2) (https://ggplot2.tidyverse.org), ggpubr (version 0.6.0; http://rpkgs.datanovia.com/ggpubr/) and ggextra (version 0.10.1; http://github.com/daattali/ggExtra) R packages were utilized to determine the association between the immune infiltration and the model. According to the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/), the half-maximal inhibitory concentration (IC50) differences across groups were compared, and the results are presented as box plots.

Using the ConsensusClusterPlus (version 1.66.0; http://bioconductor.org/packages/ConsensusClusterPlus/) in R software, consensus clustering was carried out utilizing the expression of key prognostic SenRLs. The pheatmap (version 1.0.12; http://cran.r-project.org/package=pheatmap) R package was used to build heatmaps for clustering. PCA and t-distributed stochastic neighbor embedding (t-SNE) analysis were carried out using the Rtsne (version 0.16; http://cran.r-project.org/package=Rtsne) R package. Utilizing the survival (version 3.5–7; http://cran.r-project.org/package=survival) and survminer R packages, survival analysis was performed for clusters, and survival curves were plotted. Utilizing single sample GSEA with the GSVA (version 1.48.3; http://bioconductor.org/packages/GSVA/) R package, enrichment scores of immunological pathways were calculated. The TME of the two clusters was scored and compared using the estimate (version 1.8.2; http://bioconductor.org/packages/estimate/), limma and BiocManager (version 1.30.22; http://cran.r-project.org/package=BiocManager) R packages. Immunocheckpoint gene analysis was performed for both clusters using ggplot2 and ggpubr R packages.

To further validate the reliability and interpretability of the established prognostic SenRLs signature, multi-model machine learning and SHAP analysis were performed. A total of four distinct classification models were constructed: Random Forest (RF), Gradient Boosting Machine (GBM), Support Vector Machine (SVM) and K-Nearest Neighbors (KNN). Subsequently, SHAP analysis was conducted using the SHAP Python package (version 0.44.0; http://github.com/slundberg/shap).

Survival analyses utilizing the Kaplan-Meier method were performed for different groups or clusters. A weighted test (Renyi test) was used to compare the survival differences between the high and low lncRNAs expression groups. Clinical information was analyzed through GraphPad Prism (version 8.0; Dotmatics) using the χ² test. Statistical analysis was performed utilizing R software (version 4.3.1; RStudio, Inc.). P<0.05 was considered to indicate a statistically significant difference.

Co-expression analysis among cellular senescence-associated genes and lncRNAs revealed 1,619 SenRLs, visualized as a network (Fig. 1A). The expression differences of SenRLs in STAD and normal samples were examined [|log2 fold change|>1 and false discovery rate (FDR)<0.05]. The findings demonstrated that there were variations in the expression levels of 789 SenRLs, with 192 lncRNAs exhibiting low expression and 597 lncRNAs exhibiting high expression. The results were visualized with a volcano plot (Fig. 1B). Heatmaps were constructed to demonstrate the expression density of 100 of these SenRLs (Fig. 1C). The survival of patients with STAD was notably correlated with 11 SenRLs, according to univariate Cox regressions and Kaplan-Meier analysis (P<0.05; KM<0.05; Table SII). A total of six of these were protective lncRNAs [hazard ratio (HR)<1], and the other five were risk-related lncRNAs (HR>1) (Fig. 1D). Heatmaps were used to display the expression of 11 SenRLs in patients with STAD and normal samples (Fig. 1E).

A total of five key prognostic SenRLs (AL139147.1, LINC02057, AC093801.1, AL353804.2 and AC005363.2) were screened by multivariate Cox regression analysis (P<0.05; Table SIII). AL353804.2 and AC005363.2 were protective lncRNAs, and the remaining three were risk-related lncRNAs (Fig. 2A). Kaplan-Meier survival curves were plotted, utilizing the expression of each of these five lncRNAs. The Renyi test was employed to compare survival differences between high- and low-expression groups, given its robustness in handling survival curve crossovers. The high expression group of AL139147.1, LINC02057 and AC093801.1 showed poorer overall survival (OS) (P<0.05), as demonstrated in Fig. 2B-D. The AL353804.2 and AC005363.2 high-expression groups had a prolonged survival time (P<0.05; Fig. 2E and F).

Based on these five lncRNAs, the prognostic SenRLs signature was established. Risk score=(coef 1.81069714813647 × expr AL139147.1) + (coef 0.950768101588694 × expr LINC02057) + (coef 1.17094750602435 × expr AC093801.1) + (coef-1.06513974991373 × expr AL353804.2) + (coef-1.25744132430942 × expr AC005363.2). coef denotes the regression coefficient derived from multivariate Cox regression analysis, and expr indicates the expression level of each SenRL. The allocation of risk score distribution, survival status and expression of the five lncRNAs were determined in the two groups of patients from the TCGA dataset (Fig. 3A-I). The results indicate that the risk score determines the prognosis of STAD, and patients with high-risk values had worse survival status. In the high-risk group, the expression of LINC02057, AL139147.1 and AC093801.1 was notably greater, while the expression of AC005363.2 and AL353804.2 was decreased. Kaplan-Meier survival analysis was performed in the present study for the training set, test set and whole set, respectively (Fig. 3J-L). The results all displayed that the OS of patients in the high-risk group was significantly lower (P<0.05). For the training set, test set and whole set, the 5-year AUC values were 0.828, 0.703 and 0.772, respectively (Fig. 3M-O). The 5-year ROC of the SenRLs signature had notable predictive efficacy.

As demonstrated in the univariate Cox regression analysis (Fig. 4A), the risk score was significantly associated with patient prognosis (P<0.001). Furthermore, the multivariate Cox regression analysis (Fig. 4B), which adjusted for other clinical covariates, confirmed that the risk score remained a significant prognostic indicator (P<0.001). These results suggest that the SenRLs signature is an independent factor affecting patient prognosis. The AUC of the risk score was 0.752, which was greater than that of the other medical factors (Fig. 4C), indicating the notable accuracy of the model. Next, PCA was performed based on all the SenRLs identified via univariate Cox regression and the 5 key prognostic SenRLs, respectively. As shown in Fig. 4D, all the SenRLs could not effectively discriminate between the high and low-risk groups; however, based on the key prognostic SenRLs signature, there was a significant distribution distinction between the two groups (Fig. 4E); therefore, the PCA findings further support the accuracy of the model.

The clinical indicators (including age, gender, tumor grade, and TNM stage) were subdivided, and Kaplan-Meier survival curves were plotted for the subgroups stratified by these clinical indicators. The model was found to apply to clinical factors such as age, sex, grade 2, grade 3, T3 and T4 staging (Fig. 5). However, there was no statistical difference in OS in some subtypes, such as grade 1, T1 and T2. This may be related to the overly refined subgroup typology and the small number of patient cases in the subgroup. Next, the relationship between risk values and clinicopathological parameters was analyzed to learn more about the function of the SenRLs signature in the development of STAD. Patients who developed distant metastases had significantly higher risk values than those without distant metastases (P<0.001; Fig. 6A). In addition, the survival status was lower for patients with STAD with high-risk values (P<0.001; Fig. 6B). Risk values in STAD are notably linked to distant metastases and death, and the SenRLs signature may be involved in STAD disease progression.

To achieve individualized prognosis prediction of patients with STAD, a nomogram was constructed combining risk score and clinical variables (Fig. 7A). There was a high agreement between the expected and actual observed values, as demonstrated by the calibration curves, which revealed that the predicted curves for survival were all around the standard curves (Fig. 7B-D). Therefore, the nomogram combining medical factors and risk score is effective in forecasting patient OS.

The present study quantified the expression of five key prognostic lncRNAs using RT-qPCR in 15 pairs of STAD tissue specimens and paracancerous tissue specimens in order to further test the viability of the SenRLs signature. The results are shown in Fig. 8, and the expression of LINC02057, AL139147.1 and AC093801.1 was significantly greater in STAD than in paracancerous tissues. The expression of AC005363.2 and AL353804.2 was significantly decreased in STAD tissues. The tissue expression profiles of these lncRNAs are consistent with the bioinformatics findings, thus confirming the reliability and accuracy of the above analysis results.

Next the specific signaling pathways involved in prognostic SenRLs was investigated by enrichment analysis. Through GO analysis, it could be seen that the signaling pathways, such as ‘collagen-containing extracellular matrix’, ‘peptide receptor activity’ and ‘glycosaminoglycan binding’, were enriched in the high-risk group. The low-risk group was enriched for signaling pathways such as ‘mRNA splice site selection’, ‘mRNA 5 splice site recognition’ and ‘spliceosomal complex assembly’ (P<0.05; FDR<0.25; Fig. 9A). Further analysis by KEGG showed that signaling pathways such as the ‘intestinal immune network of IgA production’, ‘cytokine-cytokine receptor interaction’ and ‘cell adhesion molecules cams’ were enriched in the high-risk group (P<0.05; FDR<0.25; Fig. 9B). No significant differences in KEGG gene concentration were observed in the low-risk group. The immune response-related signaling pathways were markedly enriched in the high-risk group, suggesting that prognostic SenRLs may influence STAD development through immune pathways. Next, immunological analysis was performed in the model.

The relative proportions of 22 immunological infiltrating cells were analyzed using CIBERSORT, and the results were displayed as a bar chart (Fig. S1). A correlation analysis of immune cells was performed using seven working platforms, including XCELL, TIMER and QUANTISEQ. As shown in the bubble plot in Fig. 10A, the high-risk group had more immune cell infiltration, such as Monocyte, Myeloid dendritic cell, B cell naive, T cell CD8+, Macrophage, NK cell and B cell (P<0.05, Table SIV). The study further analyzed the correlation between risk score and immune cell infiltration. It showed that the risk score was positively correlated with most immune cell infiltration, such as dendritic cells activated, mast cells resting, monocytes and neutrophils (P<0.05; Fig. 10B). Therefore, this demonstrated that the SenRLs signature-based risk score is effective in distinguishing different features of immune cells in STAD.

Immune checkpoint inhibitor (ICI) selectively enhances host immune response to malignancies by modulating the TME (30). It was observed how immunological checkpoints were expressed in patients with STAD. The high-risk group had high expression levels of the majority of immune checkpoint genes, including CD40LG, CD200 and ICOSLG (P<0.05; Fig. 10C). This implies that patients with STAD can be grouped by the SenRLs risk patterns to select appropriate checkpoint inhibitors for them. The IC50 values for the chemotherapeutic medications methotrexate and mitomycin C differed significantly between the two groups, with the high-risk group having higher IC50 values (P<0.01; Fig. 10D and E). Therefore, the SenRLs signature could be used as a potential predictor of drug sensitivity for oncology treatment.

Clustering analysis of STAD samples was performed depending on the SenRLs signature to improve the characterization of the subtypes of cold and hot tumors in the samples. Combining the matrix clustering heatmap, cumulative distribution function (CDF) and area change beneath the CDF plot, the sub-groups at K=2 (the highest intra-group correlation) were selected as the clustering result in this study (Fig. 11A and B and Fig. S2). The samples were divided into 2 groups, named Cluster 1 and Cluster 2, respectively. Both t-SNE and PCA showed that the above division could better classify the samples into two clusters (Fig. 11C and D). The OS of the two clusters was compared, and the OS of Cluster 2 was superior to Cluster 1, but no significant difference was found (P=0.188; Fig. S3). Cluster subtyping cannot be used as a prognostic indicator, and the risk score of the SenRLs signature remains the main OS predictor. In addition, charts were developed to understand the correlation of clustering with risk groups, and according to Fig. 11E, Cluster 1 is linked to high risk, whereas Cluster 2 is linked to low risk.

The examination of the various platforms showed discrimination in the immune cell infiltration of the clusters (P<0.01; Table SV). The enrichment scores for the immunological pathways and immune cells were calculated. These findings revealed that a variety of immune cells, such as regulatory T cells (Treg), T follicular helper cells (Tfh), T helper (Th)1 cells, dendritic cells (DCs), immature DCs (iDCs), plasmacytoid DCs (pDCs) and CD8+ T cells, had elevated enrichment scores in Cluster 2 (P<0.05; Fig. 11F). T cell co-stimulation and other immune function enrichment fractions were notably greater in Cluster 2 (P<0.01; Fig. 11G), suggesting that Cluster 2 is more important in immune cell function regulation. Compared with Cluster 1, Cluster 2 had significantly greater immune, stromal and estimate scores (P<0.01; Fig. 11H-J). The two clusters have different TME, and immunotherapy is dependent on TME (31). Most of the immune checkpoint genes, including PDCD1LG2 (also called programmed death ligand-2, PD-L2), HAVCR2 (also called T-cell Immunoglobulin and Mucin Domain-containing Protein 3, TIM3) and CD274 (also called PD-L1), were upregulated in Cluster 2 (P<0.05; Fig. 11K). A tumor is defined as ‘hot’ if it contains abundant tumor-reactive immune cells, displays an active immune microenvironment and exhibits high expression of immune checkpoints (8,9). In the present study, Cluster 2 showed marked immune cell infiltration, a more active TME (with greater enrichment of immune function-related pathways and higher immune/stromal/estimate scores), and elevated immune checkpoint expression. Based on these canonical immunological features, Cluster 2 is likely to be more responsive to immunotherapy and can therefore be classified as a hot tumor. By contrast, Cluster 1 exhibited immunotherapy resistance and can be categorized as a cold tumor.

To verify the robustness of the model, cross-validation was performed using four different algorithms, and all AUC values were >0.7. The RF and GBM models achieved the highest AUC of 0.849, indicating excellent and stable predictive performance (Fig. S4A). The mean absolute SHAP value of each SenRL served as a metric for its relative importance; AL139147.1 exhibited the highest value (0.0412) (Fig. S4B). This ranking of importance was consistent with the results of the multivariate Cox regression, in which AL139147.1 showed the highest coefficient (1.811), validating the logical consistency of the signature. Additionally, the SHAP waterfall plot visualized the additive prediction process (Fig. S4C). Starting from a base value of E[f(x)]=0.913, the final prediction reached f(x)=0.978. AL139147.1, AC093801.1 and LINC02057 emerged as positive contributors, whereas AC005363.2 and AL353804.2 contributed negatively, a pattern consistent with the risk model. Collectively, these data confirm that the model relies on distinct, biologically plausible feature contributions rather than randomness, substantiating its validity and interpretability.