Gene expression data, somatic mutation data, copy number variation (CNV) data, and corresponding clinicopathological information of gastric cancer patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The GSE84337 dataset was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). After excluding patients with missing survival time, 804 samples were included in this study, 371 from the TCGA-STAD dataset and 433 from the GSE84437 dataset. To eliminate batch effects of different datasets, we converted fragments per kilobase million (FPKM) values of the TCGA-STAD dataset to transcripts per kilobase million (TPM) and merged two datasets using the ComBat algorithm of the R package sva (Leek et al., 2012). 54 TSRGs were obtained from the MSigDB (HALLMARK_TGF_BETA_SIGNALING) (Supplementary Table S1) (Yu et al., 2022).

We performed differential expression analysis of 54 TSRGs in gastric cancer samples and normal samples using R package limma with the adjusted p < 0.05 and | log2 FC)|>1 (Ritchie et al., 2015). The protein-protein interaction network of 54 TSRGs was constructed in the STRING database (https://string-db.org/). R package maftools was utilized to map the somatic mutation waterfall of 54 TSRGs in gastric cancer patients (Mayakonda et al., 2018). Lastly, we calculated the CNV gain or loss percentage of 54 TSRGs in gastric cancer patients and analyzed the chromosomal location using the R package RCircos (Zhang et al., 2013).

We first extracted the expression of 54 TSRGs in 804 samples and then performed consensus unsupervised clustering analysis based on 54 TSRGs expression levels using the R package ConsensusClusterPlus (Wilkerson and Hayes, 2010). PCA was performed to visualize the distribution between the two different TGF-β subgroups. To explore the clinical significance of different TGF-β subgroups, we performed Kaplan–Meier survival analysis using the R package survival and survminer (Wang et al., 2020). In addition, we mapped the expression heat map of 54 TSRGs using the R package pheatmap in conjunction with the clinicopathological features of the patients.

To explore the differences in TME between TGF-β subgroups, we first analyzed the stromal score, immune score, and ESTIMATE score between two subgroups using the ESTIMATE algorithm. We analyzed the differences in the expression of critical immune checkpoints such as PD-1, PD-L1, and CTLA-4 between the two subgroups. Subsequently, we calculated the infiltration level of 22 immune cells in each sample using the CIBERSORT algorithm (Newman et al., 2015) and analyzed the abundance of immune cell infiltrates between the two subgroups using the single sample gene set enrichment analysis (ssGSEA) algorithm (Zeng et al., 2022). In addition, gene set variation analysis (GSVA) was performed with the hallmark gene set (h.all.v7.5.1.symbols) to investigate the differences in TGF-β subgroups in signaling pathways (Hänzelmann et al., 2013).

To identify DEGs in the distinct TGF-β subgroups, R package limma was utilized with |log2-fold change (FC)| ≥ 1 and adjusted p < 0.05. Based on the DEGs, we performed gene ontology (GO) enrichment analysis and kyoto encyclopedia of genes and genomes (KEGG) signaling pathway analysis. We performed a clustering analysis based on the expression of DEGs and performed a Kaplan–Meier survival analysis among gene subgroups. In addition, we combined TGF-β subgroups, gene subgroups, and clinicopathological features of patients to map the expression heat map of DEGs.

To quantify the degree of risk for each patient, we constructed a risk model based on DEGs. First, we performed univariate regression analysis to screen DEGs associated with the prognosis of gastric cancer patients. Second, we randomly divided the patients into training and testing sets in a 1∶1 ratio (Qing et al., 2022). The training set is used to construct the risk model, and the testing set and the entire set are used to validate the risk model. Third, the LASSO Cox regression analysis was executed in the training set to reduce overfitting genes with 10-fold cross-validation and 1000 repeated times (Tibshirani, 1997). Finally, we performed a multivariate regression analysis using the genes screened by the LASSO regression analysis and calculated the risk score for each patient according to expression levels and regression coefficients of genes. The formula was as follows: Risk score = β _gene1 × exp _gene1 + β _gene2 × exp _gene2 + … + β _genen × exp _genen (Qing et al., 2022). Patients were divided into high- and low-risk groups based on the median risk score. Furthermore, we analyzed the relationship between the TGF-β cluster, gene cluster, risk score, and survival status using the R package ggalluvial and the differences in risk scores between distinct subgroups (Zeng et al., 2022). In the training and validation sets, we performed Kaplan-Meier survival analysis with the R package survminer and survival (Wang et al., 2020) and ROC curve analysis with the R package timeROC (Zeng et al., 2022), respectively.

To explore the performance power of the risk score among different subgroups of clinicopathological characteristics, we first analyzed the correlation between risk scores and clinicopathological characteristics using the Student’s t-test. In addition, Kaplan–Meier survival analysis was performed in different subgroups stratified by age (≤65 years or >65 years), sex (female or male), T stage (T1-2 or T3-4), and N stage (N0 or N1-3).

Univariate and multivariate Cox regression analyses were performed to explore whether the risk score could be an independent prognostic factor for gastric cancer patients in the training, testing, and entire set, respectively. Age, gender, tumor size (T), lymph node metastasis (N), and risk score were included for analysis. In addition, we constructed a nomogram integrated the risk score and clinicopathological factors to predict the survival of gastric cancer patients at 1-, 3-, and 5-year using R package rms in the training set, testing set, and entire set, respectively (Zeng et al., 2022). Calibration curves were plotted to determine the performance of the nomograms in predicting OS.

To explore the differences in the tumor immune microenvironment between high- and low-risk groups of gastric cancer patients based on the risk model, we first analyzed the stromal score, immune score, and ESTIMATE score between the two groups using the ESTIMATE algorithm. Then, we analyzed the Spearman correlation between the risk score and immune cells using seven methods, including the XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT algorithms (Zeng et al., 2022). We further analyzed the Spearman correlation between the expression of 5 genes in the model and immune cells. In addition, the ssGSEA was subjected to calculate the infiltrating immune cells’ scores and assess the activity of immune-related pathways between high- and low-risk groups using the R package gsva (Hänzelmann et al., 2013). Finally, we analyzed the expression levels of immune checkpoint-related genes between high- and low-risk groups.

TMB(Rizzo et al., 2021), MSS(Rizzo et al., 2021), IPS(Wu et al., 2021), and TIDE (Zeng et al., 2022) scores were considered markers to predict immunotherapy response. First, we downloaded the mutation data of gastric cancer patients in MAF format from the TCGA database and annotated them using the R package maftools (Mayakonda et al., 2018), and subsequently analyzed the correlation between the risk score and TMB as well as the mutated genes common to patients in high- and low-risk groups. Second, we downloaded IPS and MSS data from the TCIA database (http://tcia.at/) for gastric cancer patients and analyzed the differences between patients in high- and low-risk groups. Finally, we analyzed the response rate of gastric cancer patients to immunotherapy based on the TIDE website (http://tide.dfci.harvard.edu/).

Next, we used the R package pRRophetic to calculate the half inhibitory centration (IC50) of antitumor drugs for each patient and analyzed the differences in sensitivity to antitumor drugs between patients in high- and low-risk groups (Geeleher et al., 2014).

R software (version 4.1.2) and related R packages were utilized for statistical analyses. The Wilcoxon test was used to compare clinicopathological characteristics, immune status, TMB, IPS, TIDE scores, and IC50 values between different groups. Kaplan-Meier curves were used to compare OS between different groups. Univariate and multivariate Cox regression analyses were used to analyze independent prognostic factors. ROC curves and nomograms were used to evaluate the predictive power of the risk model. p < 0.05 was considered statistically significant. *p < 0.05; **p < 0.01; ***p < 0.001.

The design idea of this study is shown in Supplementary Figure S1. We first performed differential expression analysis of 54 TSRGs in gastric cancer tissues and normal gastric tissues. We obtained 43 DEGs with the adjusted p < 0.05 and | log2 FC)|>1, of which JUNB, ID1, CDKN1C, ID3, and BCAR3 were lowly expressed in gastric cancer tissues, and the remaining DEGs were highly expressed in gastric cancer tissues (Figure 1A). Protein-protein interaction network analysis based on the String database revealed a close linkage between most genes (Figure 1B). Next, we explored the somatic mutation levels and the frequency of CNVs alteration in 54 TSRGs in gastric cancer patients. The waterfall plot in Figure 1C shows that 197 (45.5%) of the 433 samples had TSRG mutations. Among them, APC (11%) had the highest mutation frequency, followed by CDH1 (8%) and NCOR2 (6%). Missense mutations are the most common form of mutation in TSRGs. We also investigated the frequency of CNVs alterations of TSRGs and found that FURIN, SKIL, and ARID4B had the most significant copy number increase, while HIPK2, ID3, and BMPR1A had the most significant copy number deletion (Figure 1D). Figure 1E shows the site of CNVs of TSRGs on chromosomes.

To understand the expression pattern of TSRGs involved in tumorigenesis, data from 804 gastric cancer samples from TCGA-STAD and GSE84437 datasets were enrolled in our study for further analysis (Supplementary Table S2). To explore the characteristics of 54 TSRGs expression profiles in gastric cancer, we performed unsupervised clustering analysis to identify gastric cancer subtypes based on 54 TSRGs expression levels. The results showed that K = 2 was the most appropriate cluster, and 804 gastric cancer patients were classified into TGF-β cluster A (n = 443) and TGF-β cluster B (n = 361) (Figures 2A–C and Supplementary Table S3). The PCA results further demonstrate the excellent grouping effect (Figure 2D). Kaplan-Meier survival analysis showed a more significant survival advantage for TGF-β cluster B (p < 0.001, Figure 2E). In addition, we combined TGF-β subgroups and clinicopathological features of gastric cancer patients to map 54 TSRGs expression heatmaps and found that 54 TSRGs were expressed at higher levels in TGF-β cluster A compared to TGF-β cluster B (Figure 2F).

To explore the correlation between TSRGs and TME in gastric cancer, we first performed an ESTIMATE analysis. The results showed that patients in TGF-β cluster A had a higher stromal score, immune score, and ESTIMATE score (Figures 3A–C), suggesting that gastric cancer patients in the TGF-β cluster A have higher immune activity and lower tumor purity. Then, expression analysis of three crucial immune checkpoint genes (PD1, PD-L1, and CTLA4) showed higher expression levels of PD1, PD-L1, and CTLA4 in gastric patients in the TGF-β cluster A compared to patients in TGF-β cluster B (Figures 3D–F). We further analyzed the level of infiltration of 23 immune cells in patients with two distinct TGF-β clusters using the CIBERSORT algorithm. As shown in Figure 3G, the infiltration levels of activated B cell, activated dendritic cell, CD56 bright natural killer cell, eosinophil, gamma delta T cell, immature B cell, immature dendritic cell, MDSC, macrophage, mast cell, natural killer T cell, natural killer cell, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell, type 1 T helper cell, and type 2 T helper cell were higher in the TGF-β cluster A than those in the TGF-β cluster B, while activated CD4 T cell and neutrophil had significantly lower infiltration in TGF-β cluster A than those in the TGF-β cluster B. In addition, GSVA enrichment analysis revealed multiple tumor-associated signaling pathways enriched in TGF-β cluster A, including KRAS, IL2/STAT5, inflammatory response, hypoxia, apoptosis, and wnt/β-catenin signaling pathways (Figure 3H).

To further explore the potential biological functions of the TGF-β clusters, we obtained 202 TGF-β clusters-related DEGs (Supplementary Table S4) using R package limma and performed functional enrichment analysis. These TGF-β cluster-related DEGs are mainly enriched in biological processes associated with the extracellular matrix (Figure 4A). KEGG analysis showed that DEGs were associated with metastasis and tumor-related signaling pathways (Figure 4B), suggesting that TSRGs play an essential role in tumorigenesis and metastasis. Then, 202 TGF-β cluster-related DEGs were subjected to univariate Cox regression analysis to screen for genes associated with OS in gastric cancer. We obtained 199 genes related to the prognosis of gastric cancer patients at p < 0.05 (Supplementary Table S5). To further explore the potential mechanisms of prognosis-related DEGs in gastric cancer, based on the expression level of 199 prognostic genes, unsupervised consensus clustering analysis was utilized to classify gastric cancer patients into three different gene clusters, namely gene cluster A, gene cluster B, and gene cluster C (Supplementary Table S6). Kaplan-Meier survival analysis showed that patients in gene cluster A had the worst OS, whereas patients in gene cluster C showed a superior OS (Figure 4C). In addition, we combined the TGF-β cluster, gene cluster, and clinicopathological features of gastric cancer patients to map heat maps and found significant expression differences among gene clusters (Figure 4D). The three gene clusters showed significance in TSRGs expression, as expected from the TGF-β clusters (Figure 4E).

To quantify the risk of each gastric cancer patient, we constructed a prognostic risk model based on TGF-β cluster-related prognostic DEGs. First, the R package caret was used to randomize patients into a training set (n = 402) (Supplementary Table S7). And a testing set (n = 402) (Supplementary Table S8) at a ratio of 1:1. Second, in the training set, LASSO and multivariate Cox regression analyses were used to construct an appropriate risk model. Based on the minimum partial likelihood deviance, 12 potential candidate genes were screened by LASSO regression analysis (Figures 5A,B; Supplementary Table S9). Subsequent multivariate Cox regression of 12 prognosis-related genes yielded five genes used to construct the risk model, namely SRPX2, SGCE, DES, MMP7, and KRT17. We calculated the risk score for each patient based on the formula. Risk score= (0.1586×expression of SRPX2) + (0.1438×expression of SGCE) + (0.0728×expression of DES) + (0.0554×expression of MMP7) + (0.0754×expression of KRT17) (Figure 5C). The Sankey diagram showed the correlation between risk score and TGF-β clusters, gene clusters, and survival status (Figure 5D). In addition, we observed an obvious difference in the risk score of the TGF-β clusters and gene clusters (Figures 5E,F). The previous survival analysis showed shorter OS in the TGF-β cluster A and gene cluster A groups, and our model showed the highest risk scores in TGF-β cluster A and gene cluster A groups, which further demonstrated the excellent performance of our risk model.

Next, we divided gastric cancer patients into high- and low-risk groups based on the median risk score. The risk score curve and survival status scatter plots show that the number of deaths in gastric cancer patients increases as the risk score increases (Figures 6A,B). Kaplan-Meier survival analysis showed that patients in the high-risk group had worse OS than those in the low-risk group (Figure 6C). The risk score’s 1-, 3-, and 5-year AUC values were 0.612, 0.668, and 0.694, respectively (Figure 6D). Meanwhile, we did the same analysis in two validation sets (the testing set and the entire set), respectively, and we obtained similar results (Figures 6E–L). In the IMvigor210 cohort, patients in the high-risk group had significantly lower survival than the low-risk group (Supplementary Figure S2A), which further validates the accuracy of our constructed prognostic model. Taken together, our established risk model has an excellent performance in predicting the survival outcome of gastric cancer patients.

To explore the correlation between the risk score and available clinicopathological characteristics, we first analyzed differences in risk scores across clinical subgroups. The subgroups were divided by age (≤65 years or >65 years), sex (female or male), T stage (T1-2 or T3-4), and N stage (N0 or N1-3). The results showed that the risk scores were not statistically different across age and gender subgroups (Figures 7A,B), while patients in the T3-4 and N1-3 subgroups had higher risk scores (Figures 7C,D). In addition, we performed Kaplan-Meier survival analysis for different subgroups. We found that in the age ≤65 years (Figure 7E), age >65 years (Figure 7F), female (Figure 7G), male (Figure 7H), T3-4 (Figure 7J), N1-3 (Figure 7L) subgroups of gastric cancer patients, the OS of patients in the high-risk group was significantly lower than that of low-risk patients, while no significant differences were seen for T1-2 (Figure 7I), N0 (Figure 7K) subgroups.

To explore whether the risk score is an independent prognostic factor for patients with gastric cancer, we performed univariate and multivariate Cox regression analyses in the training set and two validation sets (testing set and entire set) in combination with clinicopathological characteristics. In the training set, univariate Cox regression analysis displayed that age (HR = 1.025, 1.011–1.039, p < 0.001), T stage (HR = 1.233, 1.011–1.505, p = 0.039), N stage (HR = 1.472, 1.256–1.725, p < 0.001), and risk score (HR = 2.122, 1.675–2.690, p < 0.001) predicted worse OS (Figure 8A). Multivariate Cox regression analysis showed that the age (HR = 1.028, 1.018–1.038, p < 0.001), N stage (HR = 1.392, 1.181–1.604, p < 0.001) and risk score (HR = 2.005, 1.562–2.574, p < 0.001) were independent prognostic factors in gastric cancer patients (Figure 8B). In the testing set, univariate Cox regression analysis displayed that age (HR = 1.027, 1.012–1.042, p < 0.001), T stage (HR = 1.276, 1.051–1.550, p = 0.014), N stage (HR = 1.633, 1.387–1.923, p < 0.001), and risk score (HR = 1.722, 1.337–2.217, p < 0.001) predicted worse OS (Figure 8C). Multivariate Cox regression analysis showed that the age (HR = 1.033, 1.018–1.048, p < 0.001), N stage (HR = 1.576, 1.334–1.863, p < 0.001) and risk score (HR = 1.674, 1.293–2.166, p < 0.001) were independent prognostic factors in gastric cancer patients (Figure 8D). In the entire set, univariate Cox regression analysis displayed that age (HR = 1.026, 1.016–1.036, p < 0.001), T stage (HR = 1.255, 1.093–1.442, p = 0.001), N stage (HR = 1.549, 1.383–1.735, p < 0.001), and risk score (HR = 1.922, 1.617–2.285, p < 0.001) predicted worse OS (Figure 8E). Multivariate Cox regression analysis showed that the age (HR = 1.028, 1.018–1.038, p < 0.001), N stage (HR = 1.475, 1.312–1.659, p < 0.001) and risk score (HR = 1.819, 1.519–2.179, p < 0.001) were independent prognostic factors in gastric cancer patients (Figure 8F). Taken together, the risk score is an independent prognostic factor for patients with gastric cancer.

Given the close correlation between the risk score and prognosis of gastric cancer patients, we integrated gender, age, T stage, N stage, and risk score to plot a nomogram to predict the 1-, 3-, and 5-year survival rates in the training set and two validation sets (testing set and entire set) (Figures 9A,C,E). Furthermore, the 1-, 3-, and 5-year calibration curves showed great accuracy between the nomogram-predicted OS and the actual observed OS (Figures 9B,D,F).

To explore the differences in tumor immune microenvironment between high- and low-risk groups of gastric cancer patients, we first performed ESTIMATE analysis. The results showed that gastric patients in the high-risk group had a higher stromal score, immune score, and ESTIMATE score (Figure 10A). Subsequently, seven algorithms were used to assess the correlation between the level of immune cell infiltration and the risk score. As shown in Figure 10B, the risk score was positively correlated with myeloid dendritic cell, CD4⁺ T cell, CD8⁺ T cell, cancer-associated fibroblast, hematopoietic stem cell, neutrophil, and macrophage M2, while negatively correlated with T cell CD4⁺ memory activated, T cell follicular helper, NK cell resting, and mast cell resting (Supplementary Table S10). We also performed a correlation analysis between the five genes in our prognostic model and the immune cells. We found that DES, KRT17, SGCE, and SRPX2 were significantly correlated with most immune cells, while MMP7 only correlated with macrophages M1 and eosinophils (Figure 10C). In addition, we further explored the difference of 16 immune cells and 13 immune-related pathways between the two subgroups using ssGSEA. We found that B cells, DCs, iDCs, macrophages, mast cells, neutrophils, TIL, CCR, HLA, parainflammation, type I IFN response, and type II IFN response were more enriched in the high-risk group, while the Th1 cells, Th2 cells, APC co inhibition, and MHC class I is less enriched in the high-risk group (Figures 10D,E). Finally, we analyzed the expression levels of immune checkpoint-related genes between two subgroups. Figure 10F showed that 24 immune checkpoint-related genes were differentially expressed in the high- and low-risk groups.

TMB and MSI are considered biomarkers of tumor immunotherapy response rate (Rizzo et al., 2021), and patients with high TMB and MSI-H benefit from immunotherapy and have more prolonged survival. Therefore, we first analyzed the correlation between the TMB and risk score. The results showed a negative correlation between the TMB and risk score (Figure 11A), and the TMB of gastric cancer patients in the low-risk group was significantly higher than that of gastric cancer patients in the high-risk group (Figure 11B). Kaplan-Meier survival analysis showed that the risk score diminished the prognostic advantage of patients with gastric cancer in the high TMB group (Figure 11C). We further analyzed the somatic mutations in the high- and low-risk groups of gastric cancer patients. The results showed that the most common form of mutation was missense mutation, and the top five mutated genes were TTN, TP53, MUC16, ARID1A, and LRP1B, and the frequency of mutations was higher in the low-risk group (Figure 11D), which was consistent with the results of the above study. In addition, we analyzed the correlation between the MSI and risk scores and showed that patients in the low-risk group had a higher proportion of MSI-H and that patients with MSI-H had lower risk scores (Figure 11E).

IPS and TIDE scores are novel tumor immunotherapy response rate biomarkers that better assess the efficacy of anti-PD1 and anti-CTLA4 therapies. A high IPS score represents higher immunogenicity, and a high TIDE score represents a greater likelihood of tumor immune escape (Wu et al., 2021; Zeng et al., 2022); therefore, the higher the IPS and the lower the TIDE score, the better the patient’s outcome to immunotherapy. Our results showed that gastric cancer patients in the low-risk group had higher IPS (Ips_ctla4_neg_pd1_neg, ips_ctla4_pos_pd1_neg, and ips_ctla4_pos_pd1_pos scores) than those in the high-risk group, but there was no statistically significant difference between the two groups in the ips_ctla4_neg_pd1_pos score (Figures 11F–I). Furthermore, the TIDE, dysfunction, and exclusion scores of gastric cancer patients in the low-risk group were lower than those in the high-risk group (Figures 11J–L). In addition, analysis of immunotherapy response based on the IMvigor210 cohort showed that patients in the immunotherapy-responsive group (complete response (CR)/partial response (PR) group) had significantly lower risk scores than the immunotherapy non-responsive group (stable disease (SD)/progressive disease (PD) group) (Supplementary Figure S2B). The above results suggest that patients with gastric cancer in the low-risk group may be may be more sensitive to immunotherapy.

To explore the potential role of our established risk model for clinical treatment, we analyzed the differences in IC50 of common antitumor drugs between high- and low-risk groups. We found that gastric cancer patients in the low-risk group were more sensitive to ATRA, cytarabine, gefitinib, gemcitabine, methotrexate, metformin, paclitaxel, rapamycin, sorafenib, tipifarnib, and vorinostat than those in the high-risk group, while gastric cancer patients in the low-risk group were less sensitive to axitinib, bleomycin, bortezomib, docetaxel, doxorubicin, erlotinib, imatinib, lapatinib, and pazopanib than those in the high-risk group (Figure 12). The above results suggest that our prognostic model can be an essential indicator of antitumor drugs for patients with gastric cancer.