Gene expression profiles for 407 samples (32 normal stomach tissue samples and 375 stomach adenocarcinoma samples) and corresponding clinical and survival information were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). After excluding samples with missing survival time and survival time less than 30 days, the entire set, including 337 STAD cases, was randomly divided into a training set (n = 169, Supplementary Table S1) and a testing set (n = 168, Supplementary Table S2). The training set was utilized to build the risk model, and the testing set and entire set were used to verify the risk model. A total of 279 cellular senescence-related genes were obtained from the CellAge database (https://genomics.senescence.info/cells/, Supplementary Table S3).

Differentially expressed lncRNAs (DElncRNAs) and cellular senescence-related genes (DECSGs) between 32 normal stomach tissue samples and 375 stomach adenocarcinoma samples were obtained with adjusted p < 0.05 and | log2-fold change (FC)| > 1 using the R package limma (Ritchie et al., 2015). Next, Pearson’s correlation analysis was performed to screen cellular senescence-related lncRNAs (CSlncRNAs) based on DECSGs and lncRNAs with |R|>0.4 and p < 0.001. Finally, DECSlncRNAs were obtained by overlapping DElncRNAs and CSlncRNAs.

Univariate Cox regression analysis was performed to obtain potential prognostic DECSlncRNAs with the threshold of p < 0.05. Then, LASSO Cox regression analysis was executed to reduce overfitting lncRNAs with 10-fold cross-validation and 1,000 repeated times. The risk score of each STAD patient was calculated based on the expression levels and regression coefficients of cellular senescence-related lncRNAs. The formula was as follows: Risk score = β _lncRNA1 × exp _lncRNA1 + ß _lncRNA2 × exp _lncRNA2 + … + ß _lncRNAn × exp _lncRNAn. Patients were divided into high- and low-risk groups based on the median risk score.

To compare overall survival (OS) between the high- and low-risk groups in training, testing , and the entire set, Kaplan–Meier survival analysis was performed using the R packages survminer and survival (Zhao et al., 2021). Subgroup analysis for the OS of STAD patients was also performed based on clinicopathological characteristics. ROC curves were constructed, and the area under the curve (AUC) values were calculated using the R package survival ROC (Heagerty et al., 2000).

We performed univariate and multivariate Cox regression analyses to explore whether the risk score could be an independent prognostic factor for STAD patients. Age, gender, grade, clinical stage, tumor size (T), distance metastasis (M), lymph node metastasis (N), and risk score were included for analysis. To predict the survival of STAD patients at 1, 3, and 5 years, a nomogram integrating the risk score and clinicopathological factors was created using the R package rms (Xu et al., 2021). Calibration curves were plotted to detect the predictive performance of the nomograms for OS.

PCA is a commonly used unsupervised learning method that can reduce the dimension of multidimensional data and extract the main feature components of the data (Ringnér, 2008). To explore the distribution of high- and low-risk patients, we performed PCA based on the whole-genome, CSlncRNAs, and the CSLPS, including 15 cellular senescence-related lncRNAs.

GSEA software (version 4.1.0) was utilized to explore potential biological functions in high- and low-risk groups. The c2. cp.kegg.v7.4. symbols.gmt were used for annotated gene sets. A total of one thousand permutations were performed, and the normalized enrichment score (NES) was calculated based on the Affymetrix chip platform. Normal p-value < 0.05 and a false discovery rate (FDR q-value) < 0.25 were regarded as significantly enriched (Subramanian et al., 2005).

The ESTIMATE algorithm was utilized to explore the difference in the TME between high- and low-risk groups. XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT algorithms were utilized to analyze the correlation between immune components and risk scores based on CSLPS. Moreover, single-sample GSEA (ssGSEA) was used to analyze 16 infiltrating immune cells and 13 immune functions between high- and low-risk groups using the R package gsva (Hänzelmann et al., 2013). We also performed immune checkpoint-related gene differential expression analysis between the two subgroups.

The mutation data of STAD patients were also downloaded from the TCGA database and analyzed using the R package Maftools (Mayakonda et al., 2018). The TMB difference between the high- and low-risk groups was compared. Then, patients were divided into high- and low-TMB groups according to the best cut-off TMB values. Survival analysis was performed based on tumor mutation burden status and risk score. Furthermore, we downloaded the microsatellite status data from TCIA (http://tcia.at/) and compared the differences between high- and low-risk groups. The tumor immune dysfunction and exclusion (TIDE) algorithm was utilized to explore the immunotherapy response in STAD patients using the website (http://tide.dfci.harvard.edu/).

R software (version 4.1.2) was used for statistical analyses. The Wilcoxon test was used to compare clinicopathological characteristics, TME, TMB, MSI, and TIDE scores between high- and low-risk groups. The Kaplan–Meier curve was used to compare survival between different groups. Univariate and multivariate Cox regression analyses were used to analyze independent prognostic factors. ROC curves were used to assess the predictive power of the CSLPS. p < 0.05 was considered statistically significant. *, p < 0.05; **p < 0.01; ***p < 0.001.

The research process is shown in Figure 1. From the TCGA database, we downloaded 32 normal stomach tissue samples and 375 stomach adenocarcinoma samples. Then, the 279 genes obtained from the CellAge database were compared in 32 normal stomach tissue samples and 375 stomach adenocarcinoma samples to explore the expression of cellular senescence-related genes in STAD patients. Among them, 47 genes were upregulated, whereas 23 were downregulated (Supplementary Table S4). A total of 595 cellular senescence-related lncRNAs were obtained from coexpression analysis based on 70 cellular senescence-related genes via the criteria |R|>0.4 and p < 0.001 (Figure 2A; Supplementary Table S5). Then, 393 DECSlncRNAs were identified by overlapping with 3625 DElncRNAs in STAD (Figures 2B,C; Supplementary Table S6).

After excluding samples with a survival time of less than 30 days, 337 out of 375 STAD samples were set as the entire set. No statistical differences in clinicopathological factors were observed between the training and testing sets (Table 1). Univariate Cox regression analysis showed that 29 lncRNAs were significantly associated with OS (p < 0.05, Figure 3A). To reduce overfitting of lncRNAs, LASSO Cox regression analysis was performed, and 15 of the 29 lncRNAs were chosen to construct CSLPS based on 1,000 times 10-fold cross-validation and the optimal value (Figures 3B,C; Table 2). The Sankey diagram displayed that six were protective and nine were risk lncRNAs (Figure 3D). Among them, AC007405.2, AL033527.2, AL136115.1, AL355574.1, REPIN1-AS1, and TYMSOS are potential protective factors, but AL391152.1, AC005165.1, AC083902.1, AC104695.3, AC129507.1, ADAMTS9-AS1, AL121748.1, LINC00460, and SCAT1 are underlying hazardous indicators. Risk score was calculated according to the formula: risk score = (−0.0056 × REPIN1-AS1) + (−0.0779 × AL355574.1) + (0.1341 × AC104695.3) + (−0.5803 × AL033527.2) + (0.7814 × AC083902.1) + (−0.1343 × TYMSOS) + (0.0129 × LINC00460) + (0.0104 × AC005165.1) + (−0.6483 × AL136115.1) + (−0.2790 × AC007405.2) + (0.8521 × AL391152.1) + (0.0659 × SCAT1) + (0.7297 × AC129507.1) + (1.8695 × AL121748.1) + (−0.9227 × ADAMTTS9-AS1).

According to the median value of the risk score, STAD patients were divided into high- and low-risk groups. As shown in Figures 4A,B, the risk score was positively associated with the number of deaths. The Kaplan–Meier survival analysis indicated that STAD patients in the high-risk group had significantly shorter OS time than those in the low-risk group (Figure 4C). The 1-, 3-, and 5-year AUC values were 0.741, 0.819, and 0.865, respectively (Figure 4D). Moreover, the AUC value of the risk score at 1 year was higher than that of age, gender, grade, stage, T, M, and N (Figure 4E). At the same time, we performed the same analysis in two validation sets. Similar results were observed in the testing and entire sets (Figures 4F–O). Taken together, our established CSLPS shows good performance in predicting survival outcomes of STAD patients.

To further explore whether the CSLPS was associated with the clinicopathological features of STAD patients, we performed a subgroup survival analysis. The subgroups were divided by age (>65 years or ≤ 65 years), sex (female or male), grade (G1-2 or G3), M stage (M0 or M1), N stage (N0 or N1-3), TNM stage (stages I–II or stages III–IV), and T stage (T1-2 or T3-4). We found that the OS time of high-risk group STAD patients was significantly shorter than that of low-risk group STAD patients in all the subgroups (Figure 5).

We conducted univariate and multivariate Cox regression analyses to explore whether CSLPS could be an independent prognostic factor for STAD patients. Univariate Cox regression analysis showed that age (HR = 1.019, 1.001–1.037, p = 0.036), stage (HR = 1.496, 1.199–1.867, p < 0.001), N stage (HR = 1.315, 1.120–1.545, p < 0.001), and risk score (HR = 1.468, 1.288–1.673, p < 0.001) predicted worse OS (Figure 6A). Furthermore, multivariate Cox regression analysis verified that the risk score (HR = 1.542, 1.330–1.787, p < 0.001) based on CSLPS was an independent prognostic factor in STAD patients (Figure 6B). To further improve the predictive value of the CSLPS in STAD patients, we constructed a nomogram taking into account age, gender, stage, grade, T, N, M, and risk score to predict OS at 1, 3, and 5 years (Figure 6C). The 1-, 3-, and 5-year calibration curves demonstrated good agreement between predicted and observed OS (Figure 6D).

PCA visualization analysis based on the whole genome and CSlncRNAs showed that the distribution of the high-risk group and the low-risk group was scattered (Figures 7A,B), while visualization analysis based on the 15 lncRNAs in CSLPS showed that the high- and low-risk groups had significantly different distributions (Figure 7C). PCA further verified the grouping ability of CSLPS, including 15 CSlncRNAs. Next, GSEA was utilized to explore the potential biological functions of patients in high- and low-risk groups based on the CSLPS. The results suggested that the high-risk group was associated with the Toll-like receptor signaling pathway, JAK/STAT signaling pathway, NOD-like receptor signaling pathway, chemokine signaling pathway, and cytokine–cytokine receptor interaction (Figures 8A–E). In contrast, the low-risk group was related to glycosylphosphatidylinositol GPI anchor biosynthesis (Figure 8F).

To explore the relevance of our established CSLPS to the immune landscape, we first explored differences in the TME between high- and low-risk groups. ESTIMATE analysis showed that the high-risk group had higher stromal, immune, and ESTIMATE scores (Figure 9A). Then, we explored the correlation between risk scores and immune cell infiltration. A bubble chart based on seven different algorithms showed that the risk score was positively correlated with myeloid dendritic cells, cancer-associated fibroblasts, M2 macrophages, B cells, hematopoietic stem cells, T cell CD8⁺, and mast cells while negatively correlated with NK cells, M1 macrophage, T cell CD4⁺ Th1, and T cell CD4⁺ Th2 (all p < 0.05, Figure 9B, Supplementary Table S7). In addition, ssGESA was applied to explore the difference between the two subgroups of 16 immune cells and 13 immune-related pathways. We found that B cells, DCs, iDCs, macrophages, mast cells, neutrophils, NK cells, pDCs, T helper cells, TIL, Treg, CCR, parainflammation, and type Ⅱ IFN response were more enriched in the high-risk group, while the MHC class Ⅰ is less enriched in the high-risk group (Figures 9C,D). Finally, we analyzed the expression levels of the immune checkpoint-related genes between the two subgroups. The results showed that TNFSF14, CD28, CD276, TNFSF18, CD80, CD40LG, BTLA, LAIR1, NRP1, CD86, TNFRSF8, CD200, CD48, PDCD1LG2, and CD200R1 genes were more highly expressed in the high-risk group, while TNFSF9 and TNFRSF14 were lower expressed in the high-risk group (Figure 9E). The aforementioned findings indicate that high-risk group patients present an immunosuppressive microenvironment.

TMB and MSI were considered predictive biomarkers of tumor immunotherapy response (Cristescu et al., 2018; Liu et al., 2019). As shown in Figure 10A, the most common type of mutation in high- and low-risk group patients was missense mutation, and the top three mutated genes were TTN, TP53, and MUC16. Intriguingly, TTN and MUC16 were more likely to be mutated in the low-risk group than in the high-risk group. Moreover, TMB was negatively associated with risk scores, and STAD patients in high-risk groups had a lower TMB than those in low-risk groups (Figures 10B,C). Survival analysis showed that STAD patients with a high TMB had better outcomes than those with a low TMB (Figure 10D), and the risk score reduced the prognostic value in the high-TMB group according to the survival analysis combined TMB and risk score (Figure 10E). In addition, the low-risk group had a lower proportion of patients with microsatellite stability (MSS) and a higher proportion of patients with high microsatellite instability (MSI-H) (Figure 10F; Supplementary Table S8). In addition, the TIDE score was a novel valuable predictive biomarker for tumor immunotherapy response, and patients with lower TIDE scores could benefit from immunotherapy and have a longer survival time (Jiang et al., 2018). Interestingly, we found that STAD patients in the low-risk group had lower TIDE scores and T-cell dysfunction scores than those in the high-risk group, but no statistical difference between the two subgroups in T cell exclusion was found (Figure 10G). The aforementioned results suggest that the low-risk group of STAD patients may be more effective for immunotherapy.