The raw data was downloaded from TCGA and GEO databases. Duplicate samples and those lacking essential clinical characteristics or survival information were removed. The training cohort consisted of 410 STAD samples and 10 gastric normal samples from TCGA. The validation cohort was selected to be the GSE66229 dataset (12), which was verified using the GPL570 platform. To eliminate discrepancies caused by batch effects and ensure research integrity and reliability, COMBAT was used when merging GEO data. In the survival analyses, patients were included based on the availability of survival status and survival time, with a minimum survival time of 30 days. The list of genes associated with autophagy and senescence was obtained from GeneCards datasets (13).

To identify common genes, Venny 2.1.0 was employed. DEGs were detected using the “limma” package in R. A cutoff of |logFC| ≥ 1 and a FDR < 0.05 were applied to determine significant DEGs. A volcano plot was drawn by the “pheatmap” package in R to visually represent the DEGs.

The “clusterProfiler” package was employed to conduct GO and KEGG analyses. These analyses identified the biological functions and pathways associated with the DEASRGs, providing insights into their biological significance. Furthermore, the functional profiles of the different risk groups were assessed using Gene Set Enrichment Analysis (GSEA) to detect any relevant changes.

DEASRGs were utilized to conduct consensus cluster analysis using the “ConsensusClusterPlus” package in R. Employing optimal k-means clustering, STAD patients were categorized into three distinct groups. Principal component analysis (PCA) was implemented to differentiate these clusters. We used the “estimate” package to determine particular scores in tumor tissue for assessing the extent of infiltration by stroma and immune cells.

Univariate Cox analysis was employed to screen core prognostic DEASRGs and we further assessed their copy number variation (CNV) alterations. Subsequently, LASSO regression was utilized to select genes for constructing the prognostic signature. Through the calculation of the following algorithm, we obtained corresponding risk score for every single patient.

The variables n, Coefi, and Expi represent the signature gene number, the risk weighting coefficient index, and the expression level of the signature gene, respectively.

The median risk score in the TCGA cohort was used to distinguish the high-risk group from the low-risk group. KM survival curves were generated to compare prognosis between the groups. We constructed a ROC curve to compare the Area Under the Curve (AUC) value of the risk score and several clinical markers. Risk curves and survival status analyses were performed to evaluate the efficiency of model in high- and low-risk groups. Additionally, we conducted PCA analysis to visualize the distribution of patients.

To enhance predictive power the signature, a nomogram was established incorporating risk scores and clinical features. Variables within the nomogram, including age, gender, M stage, T stage, N stage, clinical stage, and risk score, were assigned points based on their relative prognostic contributions. Individual patient scores were calculated by summing these points. The calibration curve allowed us to evaluate the predictive capability of the nomogram across various survival periods. Decision curve analysis (DCA) was used to evaluate the clinical benefits brought by the model.

To investigate the association between the novel signature and various clinical factors, subgroup analyses were conducted within the TCGA cohort. Available clinicopathological features were extracted for further analysis. Moreover, survival curves were also plotted across distinct clinical subgroups.

Drug sensitivity assessments was performed using data obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) public database (14). The “oncoPredict” package was employed to calculate the half-maximal inhibitory concentration (IC₅₀) values for therapeutic drug.

The CIBERSORT algorithm enabled us to accurately assess the composition of infiltrating immune cell types within patient’s tumor sample (15). To assess potential treatment response based on proportions of immune cell, we examined the expression of immune checkpoint genes within the two subgroups. Furthermore, we employed the tumor immune dysfunction and exclusion (TIDE) algorithm to obtain TIDE scores, dysfunction scores, and exclusion scores for each tumor sample. Tumor purity was evaluated using the ESTIMATE algorithm. The stromal score represents the proportion of stromal cells, while the immune score reflects the balance of immune cell populations. Finally, we visualized the distribution differences between risk groups and microsatellite status (MSI) states using boxplots.

We downloaded somatic mutation data of training group via the UCSC Xena browser. This data was then visualized using waterfall charts from the “maftools” R package. Next, we computed the tumor mutational burden (TMB) score for every single sample and investigated its correlation with risk score.

Associations between these model genes were constructed by STRING database. A PPI network was built using interaction scores higher than 0.15 and P < 0.05 as the significant threshold. Genes with interaction scores greater than 0.15 were selected to construct PPI network. Cytohubba plugin was utilized to estimate the MCC score, Stress score, Degree score and Closeness score. In this analysis, genes with the same score were considered to be sequenced equally.

The raw expression profiling of GSE112302 was retrieved from the GEO dataset. The data pertaining to normal tissue were omitted, while the data corresponding to tumor tissue were utilized for subsequent analysis. The employment of the “Seurat” program was necessary for performing data quality control, PCA, and t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization are all reliant on the utilization of the “Seurat” package. The “SingleR” package was vital for annotating the cell types in each cluster.

The gastric cancer cell lines (SGC7901, AGS, HGC27) and the human normal gastric mucosal cell line GES-1 were procured from Pricella Life Science & Technology Co., Ltd. These cell lines were cultured in RPMI-1640 medium (Pricella Life Science & Technology Co., Ltd) supplemented with 10% Fetal Bovine Serum (FBS, from Shanghai Life-iLab Biotech Co., Ltd.), and containing penicillin and streptomycin.

The growth of AGS and HGC27 cells was assessed using the MTT assay. Cells were seeded at a stable density in a 24-well plate and incubated. Subsequently, 0.5 mg/mL MTT solution (M158055, Aladdin) was added to each well at 24, 72, and 120 hours, respectively. After 3-4 hours of incubation at 37°C, the supernatant was removed, and dimethyl sulfoxide was added to dissolve the formazan precipitate. The absorbance of the resulting solution was measured at 490 nm using a microplate spectrophotometer.

To assess the proliferation of AGS and HGC27 cells, a cell cloning assay was performed. Cells were seeded into a six-well plate and incubated under standard conditions for 10-14 days. Subsequently, the supernatant was removed, and cells were fixed with 4% paraformaldehyde (30072418, China National Pharmaceutical Group Chemical Reagent Co., Ltd). Following fixation, cells were stained with crystal violet, and the number of colonies formed was manually counted.

To assess the migration of AGS and HGC27 cells, a cell scratch assay was performed. Cells were seeded in a 24-well plate and incubated for 2-3 days until reaching approximately 90% confluency. A wound was created on the cell monolayer using a 200 µL pipette tip. Cell migration into the wound area was observed and evaluated after 48 hours.

Use cell lysis buffer (containing 20mM Tris (pH 7.5), 150mM NaCl, 1% Triton X-100, and other components) to lyse cell samples. Protein concentration in the collected lysates was determined using the BCA quantification method. Subsequently, 25-40 µg of protein was separated on an SDS-PAGE gel, followed by transfer to a polyvinylidene fluoride (PVDF) membrane. The membrane was incubated overnight at 4°C with primary antibodies against TXNIP and GAPDH. Afterward, the membrane was incubated with a secondary antibody at room temperature for 2 hours. Protein bands were visualized using ECL chemiluminescence.

The work in this study was primarily performed by R software (version 4.3.1). Univariate and multivariate Cox regression analysis were employed to evaluate the independent prognostic significance of variables. For survival analysis, the Kaplan-Meier method was employed to plot the survival curves of different risk groups, and the Log-rank test was utilized to evaluate the significance of survival differences between groups. Regarding the comparison of continuous data, the independent Student’s t-test was carried out for normally distributed data between two groups, while the Wilcoxon test was applied for non - normally distributed continuous variables. The Spearman’s correlation coefficient was computed to evaluate the associations between two variables. Statistical significance was determined by setting the threshold at p < 0.05 for all analyses.

Firstly, a general workflow was constructed to outline the entire analysis process ( Figure 1 ). Supplementary table 1 displayed a fundamental table containing baseline data of certain participants in this study. These participants have available survival status and survival time (≥30 days), and their T, N, M and clinical stage are clearly defined. As illustrated in the Venn plot ( Figure 2A ), 685 overlapping genes were obtained by intersecting 2269 autophagy-related genes and 4136 senescence-related genes from the GeneCards database with TCGA-STAD genes. Through differential expression analysis of these overlapping genes, 161 autophagy- and senescence-related DEGs (DEASRGs) were filtered comparing normal and tumor groups ( Figure 2B ). To provide a clearer understanding of the functional properties of DEASRGs in STAD, GO enrichment analysis was performed. The GO terms for biological processes and molecular functions revealed that the DEASRGs were mainly involved in the regulation of mitotic cell cycle phase transition, regulation of response to DNA damage stimulus, and ubiquitin-like protein ligase binding ( Figure 2C ). Additionally, KEGG pathway analysis was conducted to investigate the possible mechanistic pathways associated with these DEGs, including the cell cycle, cellular senescence, and PI3K-Akt signaling pathway ( Figure 2D ).

Through univariate Cox regression analysis, we identified 29 genes from the initial 161 DEASRGs as potential prognostic factors for STAD patients ( Figure 2E ). We investigated CNV alterations in these 29 genes, revealing predominantly copy number gains. However, COL3A1, HMGB2, DNMT1, CASP2, EZH2, DCN, TP53, LMNB2, TTF2, LMNB1, TNFRSF10B and UHRF1 exhibited a greater frequency of CNV losses ( Figure 2F ). The chromosomal location of CNV alterations for these 29 genes is depicted in Figure 2G . Finally, LASSO regression analysis further reduced the gene set to 14 ( Figures 3A, B ).

After calculating risk scores, patients were categorized into high- and low-risk groups using the median risk score as a threshold. KM survival analysis demonstrated a statistically better OS for the individuals classified as low risk ( Figures 3C, D ). The risk score exhibited a negative correlation with patient survival, indicating that a higher risk score predicts lower overall survival and higher mortality ( Figures 3E, F ). Similar analyses were conducted on the validation cohort ( Figures 3G, H ). Figures 3I, J presents gene expression profiles of the prognostic model genes as heatmaps. Our model demonstrated superior prognostic accuracy for gastric cancer patients compared to traditional clinical indicators in TCGA cohort ( Figure 3K ). Time-dependent ROC curves with AUCs of 0.645, 0.673, and 0.721 at 1, 3, and 5 years further validated the model’s efficiency ( Figure 3L ). In the GSE66229 dataset, our model outperformed most clinical indicators ( Figure 3M ), with AUCs of 0.627, 0.648, and 0.633 at 1, 3, and 5 years ( Figure 3N ). After that, we employed GSEA method to investigate disparities in biological functionality between patients classified as high and low risk ( Supplementary Figure S1 ). Interestingly, the findings indicated that the biological functions of the high-risk group were intricately linked to the composition and specific activities of the extracellular matrix, while the low-risk group was primarily enriched in cell cycle, certain activities related to DNA and mitochondria. Univariate and multivariate Cox analyses confirmed the risk score as a significantly independent predictor of gastric cancer in the TCGA-STAD and GSE66229 datasets ( Supplementary Figures S2A-D ). Next, we examined the association between risk scores and clinical factors. The findings suggested that the younger populations (<65) had considerably higher risk scores ( Supplementary Figure S2E ). In addition, the differences in risk scores among gender, T stage, N stage, M stage, and clinical stage were not statistically significant ( Supplementary Figures S2F-J ). Moreover, Supplementary Figure S3 illustrated the difference in survival rates among high- and low-risk patients across several clinical subgroups.

Based on the expression of DEASRGs, we employed the Consensus Cluster algorithm to identify three distinct patient subtypes within the TCGA-STAD cohort, designated as clusters 1, 2, and 3 ( Figure 4A ). A PCA plot visualized the transcriptional differences among the three clusters in a three-dimensional space ( Figure 4B ). Survival analyses revealed a significant survival advantage for patients in cluster 1 ( Figure 4C ). A Sankey diagram illustrated patient transitions among gene clusters, risk groups, and survival status, demonstrating higher survival rates in cluster 1 and low-risk group ( Figure 4D ). The risk score in cluster 1 exhibited a statistically significant decrease compared to the other two clusters ( Figure 4E ). We subsequently assessed inter-cluster variations in the immunological microenvironment ( Figures 4F-H ). Notably, cluster 2 exhibited considerably higher expression levels of most ICs than the other clusters ( Figure 4I ), suggesting a suboptimal response to immunotherapy. A heatmap comparing immune cell infiltration patterns across clusters using algorithms from multiple platforms is presented in Figure 4J . C2 exhibited the highest overall immune cell infiltration, consistent with the findings from the ESTIMATE algorithm ( Figure 4G ). However, elevated expression of M2 macrophages, myeloid-derived suppressor cells, and tumor-associated fibroblasts was detected in C2 across multiple platforms, suggesting that this subgroup resides in an immunosuppressive microenvironment.

To enhance the clinical applicability and predictive accuracy of our signature, a nomogram was constructed incorporating risk score and other clinical indicators ( Figures 5A, E ). Calibration curves exhibited robust concordance between the expected and observed survival probability at 1, 3, and 5 years ( Figures 5B, F ), indicating high nomogram accuracy and reliability. DCA curve revealed that the nomogram exhibited larger net benefit compared to the nomogram without prognostic signature (nomogram excluding ASRG) and other factors ( Figures 5C, G ). ROC curve demonstrated superior predictive accuracy of the nomogram compared to other factors, such as nomogram excluding ASRG, risk score, gender, age, and TNM stage, with AUC values of 0.691 in the training set and 0.826 in the validation set ( Figures 5D, H ). The above results indicated that the incorporation of prognostic signature contributed to enhancing the superiority of the nomogram.

Resistance to therapeutic medications is a common challenge in cancer therapy, often leading to poor drug efficacy and worse clinical outcomes in gastric cancer. To enhance therapeutic benefits, we figured out whether the ASRG signature could accurately predict drug sensitivity in the training cohort. By utilizing the “oncoPredict” package, we estimated IC₅₀ values for 198 drugs in all patients. The high-risk patients showed markedly elevated IC₅₀ values for Oxaliplatin, Paclitaxel, Cisplatin, Docetaxel, 5-Fluorouracil, and Afatinib, which were positively correlated with risk scores. This suggested that individuals with lower risk scores might exhibit a more favorable response to therapies containing these medications ( Figures 6A-F ). Gemcitabine, Camptothecin, KRAS (G12C) Inhibitor, Dabrafenib, and Sorafenib also exhibited increased IC₅₀ values in the high-risk group ( Figure 6G ). Conversely, SB505124, JQ1, IGF1R, JAK, and NU7441 had higher IC₅₀ values in the low-risk patients, implying a poorer response to these drugs ( Figure 6H ).

The TME, composed of diverse immune cells, cancer-associated fibroblasts (CAFs), endocrine cells, extracellular matrix (ECM) components, and other elements, significantly influences tumorigenesis. Disrupting the tumor immune tolerance feedback loop by targeting the TME is a promising strategy to enhance cancer therapy (16). To examine the correlation between our signature and immune infiltration, we employed the CIBERSORT algorithm to determine the composition of tumor-infiltrating immune cells in STAD ( Figure 7A ). Comparative analysis of immune cell distribution between high- and low-risk groups revealed significant differences. Plasma cells were notably reduced in the high-risk patients, whereas naïve B cells, activated NK cells monocytes, resting dendritic cells, and resting mast cells were increased ( Figure 7B ). Moreover, elevated expression of multiple ICs in high-risk patients suggested increased susceptibility to immune evasion ( Figure 7C ).

Subsequently, we examined the correlation between microsatellite status and risk score. Figure 7D indicated that patients with MSI-H, known for increased immunotherapy sensitivity, exhibited lower risk scores. As expected, the MSI-H prevalence was considerably higher in low-risk patients (30%) ( Figure 7E ). To predict immune system evasion, we performed TIDE analysis. As shown in Figure 7F , individuals at high risk expressed elevated TIDE scores, indicating a greater risk of immunological escape and reduced immunotherapy responsiveness. The ESTIMATE algorithm revealed significantly elevated stromal, immune, and estimate scores in the high-risk group, positively correlated with the risk score. Conversely, a negative correlation between tumor purity and the risk score was observed ( Figures 7G, H ).

Human tumors exhibit varying levels of somatic mutations collectively termed TMB, which has been linked to immunotherapy efficacy (17, 18). To investigate the association between risk score and gene mutation, we analyzed simple nucleotide variation data from TCGA. Figures 8A, B present the top 20 genes with the highest frequency of mutations in two groups. TTN, TP53, MUC16, ARID1A, and LRP1B emerged as the five most prominent mutated genes. TMB analysis revealed an inverse relationship between TMB and risk score ( Figures 8C, D ). Spearman correlation analysis further differentiated clusters based on TMB and risk score ( Figure 8E ). Significantly, the high TMB had superior survival rates compared to the group with a low TMB ( Figure 8F ). An integrated survival analysis, including both TMB and risk groups ( Figure 8G ), demonstrated that GC patients with high TMB and low risk scores presented the most favorable outcome.

To investigate the interaction of model genes, we constructed a PPI network. Next, the results of the STRING database were exported into Cytoscape for further analysis to obtain the hub gene ( Supplementary Figure S4A ). Then, by intersecting the top 7 hub genes determined by MCC, Stress, Degree and Closeness algorithms in cytohub plug-in, we identified 6 core genes ( Supplementary Figure S4B-F ). Moreover, the impact of IRAK1, SERPINE1, KIT, CXCL1, CD36, TXNIP on prognosis of STAD was analyzed by GEPIA online tools ( Supplementary Figure S5 ). The findings indicated that lower expression levels of SERPINE1, KIT, CD36 and TXNIP were associated with longer OS, while differences in IRAK1 and CXCL1 expression levels did not have a statistically significant effect on prognosis.

We selected scRNA-seq data from GSE112302 dataset for further analysis of the model gene. To ensure the reliability of the single-cell data, we applied a filter to exclude genes expressed in fewer than three cells and cells expressing fewer than 50 genes ( Supplementary Figure S6A ). The correlation of sequencing depth with mitochondrial content and gene number was shown in the Supplementary Figure S6B . Subsequently, the data was standardized and the top 1,500 genes with significant intercellular coefficients of variation were extracted for further analysis ( Supplementary Figure S6C ). We then employed PCA analysis to reduce the dimensionality of the data ( Supplementary Figure S7A ). Supplementary Figures S7B-C illustrates the characteristic genes of the initial four principal components in the PCA analysis. We chose the initial 14 PCA components with a significance level of p<0.05 for subsequent analysis ( Supplementary Figure S7D ). The t-Distributed Stochastic Neighbor Embedding (tSNE) algorithm was applied to classify the cells into six distinct clusters, illustrating the global distribution of the single-cell transcriptomes ( Supplementary Figure S8A ). Each cluster represents a distinct cell population. Supplementary Figure S8B depicted, in the form of a heat map, the top 10 genes exhibiting the most substantial variances within each cluster. The distribution and expression of key model genes are visualized in Supplementary Figures S8C, D . Furthermore, Supplementary Figure S8E visualizes the expression of prognostic genes identified through the PPI network, including SERPINE1, KIT, CD36, and TXNIP, across the clusters: IRAK1 was significantly expressed in Clusters 4 and 5; SERPINE1 and KIT were not highly expressed in any of the clusters; CXCL1 showed high expression in Clusters 1 and 5; CD36 was predominantly expressed in Cluster 5; and TXNIP was most highly expressed in Clusters 0 and 2. Cell type annotation ( Supplementary Figure S8F ) reveals that the six clusters can be broadly classified into two major cell types: Clusters 0-4 primarily represent epithelial cells, while Cluster 5 is primarily composed of monocytes. Based on these findings, we conclude that IRAK1 and CXCL1 are expressed in both epithelial cells and monocytes, CD36 is predominantly expressed in monocytes, and TXNIP is mainly expressed in epithelial cells.

We initially examined TXNIP protein expression levels in gastric cancer cell lines SGC7901, AGS, and HGC27, as well as in normal human gastric mucosal epithelial cells GES-1. TXNIP protein expression was significantly higher in AGS and HGC27 cells compared to GES-1 cells ( Figure 9A ). To investigate the biological role of TXNIP in gastric cancer, we employed lentiviral transduction to knock down TXNIP gene expression in AGS and HGC27 cells. Based on knockdown efficiency, shRNA3 was selected for subsequent experiments ( Figure 9B ). MTT assays revealed that TXNIP knockdown significantly inhibited the growth of AGS and HGC27 cells, with a marked reduction in cell viability on days three and five ( Figures 9C, D ). Furthermore, colony formation and migration assays demonstrated that TXNIP knockdown suppressed the clonal formation and migratory capacity of AGS and HGC27 cells ( Figures 9E-H ). Overall, these results indicate that TXNIP exerts an oncogenic role in gastric cancer.