RNA sequencing data (in fragments per kilobase per million (FPKM)) of TCGA-STAD (stomach adenocarcinoma) cohort were curated from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/). Thereafter, the FPKM format was converted to the transcripts per kilobase million (TPM) format for further analysis. Clinical features of GC patients were also harvested from TCGA project. The GSE15459 dataset was downloaded from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/), which was used as the external validation set.

The gene sets of TNFα-signaling were curated from the Molecular Signatures Database (MSigDB; http://www.broadinstitute.org/msigdb) (Liberzon et al., 2015). Single-sample gene set enrichment analyses (ssGSEA) derived from the gene set variation analysis (GSVA) package were presented for quantifying the activities of TNFα signaling across GC specimens (Hänzelmann et al., 2013). The ssGSEA ranked the mRNA expression in each specimen and used empirical cumulative distribution function of genes in the signature and the remaining genes to produce an enrichment score. The ssGSEA score was normalized through the Z-score.

The WGCNA package was adopted for performing co-expression analysis (Langfelder and Horvath, 2008). The expression profiling of the first 5,000 genes according to SD was included for the WGCNA. The soft-thresholding power ß was set as 3 with the pickSoftThreshold function. Additionally, the scale-free R ² = 0.85 calculated with the softConnectivity function was set as the soft-thresholding parameter for ensuring a scale-free topology network and producing a TOM matrix. Thereafter, co-expression modules were clustered. Pearson correlation analysis was carried out between the co-expression modules and TNFα score. Moreover, the correlations between module membership and gene significance were plotted. Genes in the co-expression module that presented the strongest correlation strength with the TNFα score were deemed as TNFα-derived genes.

Univariate cox regression analysis was conducted for screening prognostic TNFα-derived gene signatures (p < 0.05). Thereafter, this study input the aforementioned genes into the Least Absolute Shrinkage and Selection Operator (LASSO) analyses with the glmnet package (Engebretsen and Bohlin, 2019). Characteristic genes were screened in accordance with the optimal λ value. The TNFα-derived risk score was determined following the formula: risk score = ∑ X i * coef i, in which coef i was the coefficient, and X i was the mRNA expression of each characteristic gene. GC patients were randomly separated into training and testing sets with a 1:1 ratio. With the median value of the risk score, patients were divided into high- and low-risk groups in each dataset. Kaplan–Meier curves of overall survival (OS), disease-free survival (DFS), disease-specific survival (DSS), and progression-free survival (PFS) were conducted between high- and low-risk groups. Survival differences were estimated with log-rank tests. Time-independent receiver operating characteristic (ROC) curves were presented for evaluation of the efficacy of the risk score in predicting GC OS outcomes. Uni- and multivariate cox regression models were constructed for screening independent prognostic factors of GC.

The TNFα-derived gene signature and clinicopathological characteristics (age, T stage, N stage, M stage, and pathological stage) were input into the nomogram model in TCGA-STAD dataset. Calibration curves, ROC curves at 5-, 6- and 7-year survival, and decision curve analyses (DCA) were presented for evaluating whether this nomogram was useful as an ideal model.

GSEA was presented for comparing activated hallmark gene sets between high- and low-risk groups in TCGA-STAD cohort (Subramanian et al., 2005). For each analysis, 1,000 gene set permutations were carried out. The hallmark gene sets curated from the MSigDB project were utilized as the reference set. Additionally, the ssGSEA score of hallmark gene sets was calculated across GC tissues.

Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm (Yoshihara et al., 2013) was utilized for evaluations of immune and stromal scores in accordance with mRNA expression signatures. Immune and stromal scores represented the tumor immune and stromal infiltrations within a bulk tumor. Thereafter, the ESTIMATE score was defined by combining immune and stromal scores within tumor tissues. Through ssGSEA, the abundance of immune cells was scored within tumor tissues in accordance with mRNA expression profiles.

The T-cell dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was calculated for predicting the clinical responses to immune checkpoint inhibitors (Jiang et al., 2018). The immunophenoscore (IPS) was determined for the prediction of the responses to CTLA-4 or PD-1 inhibitors in accordance with the marker genes of MHC-relevant signatures, checkpoint molecules, immunomodulators, effector cells, and suppressor cells (Charoentong et al., 2017). All steps within the cancer immunity cycle that reflected the anticancer immune response were quantified through the ssGSEA algorithm (Chen and Mellman, 2013). The half-maximal inhibitory concentration (IC50) values of chemotherapeutic agents from the Cancer Genome Project (Geeleher et al., 2014b) were estimated utilizing the pRRophetic package (Geeleher et al., 2014a).

Somatic mutation profiling (mutation annotation format (MAF) files) of TCGA-STAD was analyzed with MuTect2 and visualized with the maftools package (Mayakonda et al., 2018). The tumor mutational burden (TMB) was determined through non-synonymous somatic mutations utilizing 38 Mb as the estimate of the exome size (Chalmers et al., 2017).

Human gastric mucosal cells GES-1 and human GC cells HGC-27, MKN-28, and AGS were purchased from the Chinese Academy of Sciences. All cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM; Hyclone, United States) containing 10% fetal bovine serum (Gibco, United States), 100 U/ml penicillin sodium, and 100 μg/ml streptomycin (Hyclone, United States ). The cell culture flask was placed in a 5% CO₂ incubator at 37°C.

Cells were washed twice lasting 2 min through PBS and resuspended by using radioimmunoprecipitation assay buffer at 4°C. The protein content was evaluated utilizing a BCA kit (Beyotime, China), in accordance with the manufacturer’s protocols. Then, 200 µl protein lysates were separated via 10% SDS-PAGE and transferred onto the polyvinylidene difluoride (PVDF) membrane. Thereafter, the membrane was incubated with TBST (TBS with 1% Tween 20) containing 5% BSA lasting 1 h at room temperature and subsequently incubated with primary antibodies targeting AKR1B1 (1/1000; ab192865; Abcam, United States), CPVL (1/1000; ab180147; Abcam, United States), CTSL (1/1000; ab200738; Abcam, United States), and GAPDH (1/10000; ab8245; Abcam, United States) overnight at 4°C. The membrane was washed by TBST lasting 5 min at room temperature, followed by incubation with horseradish peroxidase-conjugated goat anti-rabbit secondary antibodies (1/2000; ab7090; Abcam, United States) at 37°C lasting 1 h. Through an electrochemiluminescence detection kit (Bio-Rad, United States), the protein bands were developed, and the protein expression was tested with an X-ray film. The bands were quantified by ImageJ software.

Total RNA was extracted from cells utilizing RNeasy kits (Beyotime, China) and reverse transcribed with one-step RT-PCR kits (Beyotime, China) at 37°C lasting 30 min, in accordance with the manufacturer’s protocols. qPCR was conducted utilizing SYBR Green RT-PCR kits (Takara, China). The thermocycling conditions were as follows: 95°C lasting 5 min; 40 cycles of 95°C lasting 40 s, 60°C lasting 30 s, and 72°C lasting 30 s. The following primers were used for PCR: AKR1B1: 5′-TTT​TCC​CAT​TGG​ATG​AGT​CGG-3′ (forward), 5′-CCT​GGA​GAT​GGT​TGA​AGT​TGG-3′ (reverse); CPVL: 5′-TGG​AAG​GTG​ATT​GTT​TCG​CTG-3′ (forward), 5′-GTC​TCC​CTT​AGG​TGG​CAT​GGA-3′ (reverse); CTSL: 5′- CTT​TTG​CCT​GGG​AAT​TGC​CTC -3′ (forward), 5′-CAT​CGC​CTT​CCA​CTT​GGT​C-3′ (reverse); and GAPDH: 5′-GGA​GCG​AGA​TCC​CTC​CAA​AAT-3′ (forward), 5′-GGC​TGT​TGT​CAT​ACT​TCT​CAT​GG-3′ (reverse). The fold change in mRNA expressions was determined with the 2^−ΔΔCq method.

All analyses were executed through R (version 4.0.1) and GraphPad Prism (version 8.0.1) software. With Student’s or Wilcoxon test, comparisons between groups were conducted. Pearson’s or Spearman’s correlation test was utilized to evaluate the associations between variables. p < 0.05 was indicative of statistical significance.

Through the ssGSEA method derived from the GSVA package, we quantified the activities of TNFα signaling across GC tissues. In accordance with the median value of z-scores of TNFα signaling, we separated GC patients into high and low z-score groups. Kaplan–Meier curves demonstrated that GC patients with high z-scores displayed more undesirable OS outcomes than those with low z-scores (Figure 1A). This indicated that TNFα signaling might be linked to GC prognosis, which was consistent with previous research (Ju et al., 2020). This study employed the WGCNA approach to further identify TNFα signaling-derived genes in GC. First, the top 5,000 genes according to SD were included for co-expression analyses. The genes with similar expression patterns would be clustered into one module. Hierarchical clustering analyses indicated that there was no outlier among GC samples (Figure 1B). For constructing an appropriate scale-free topological overlap matrix, we calculated the scale independence and mean connectivity under diverse soft thresholds. Consequently, when soft thresholding was set as 3, the scale-free R² was 0.853, indicating the constructed co-expression network met the scale-free topology criterion (Figure 1C). Thereafter, GC samples were clustered into 16 co-expression modules (Figure 1D). To determine the correlation between co-expression modules and the TNFα score as a phenotype, we carried out a Pearson correlation analysis. Our results uncovered that the “tan” module displayed the strongest correlation to the TNFα score (R = 0.51 and p < 0.0001; Figure 1E). Moreover, we compared the gene significance of each module with the TNFα score. In particular, we noted that the “tan” module presented the highest gene significance with the TNFα score (Figure 1F), indicating that the genes in the “tan” module were prominently associated with TNFα signaling. Herein, the 80 genes in the “tan” module were considered TNFα-derived genes, which are listed in Supplementary Table S1.

To determine prognosis-related TNFα-derived genes, we conducted univariate cox regression analyses. Our results showed that 15 TNFα-derived genes displayed significant associations with GC prognosis (p < 0.05; Table 1). The aforementioned genes were input into LASSO analyses. With the optimal λ (0.0494), three genes (AKR1B1, CPVL, and CTSL) were retained following LASSO regularization (Figures 2A,B). Figure 2C showed the prognostic significance of AKR1B1, CPVL, and CTSL in GC. The risk score of each GC specimen was quantified utilizing the established formula: risk score = 0.00453439137355748 * AKR1B1 expression +0.0023541802071365 * CPVL expression +0.00307599022458496 * CTSL expression. With the increase in the risk score, the expressions of AKR1B1, CPVL, and CTSL were gradually increased in all GC patients (Figure 2D). In accordance with the median value of the risk score, GC patients were separated into high- and low-risk groups (Figure 2E). We noted there were more patients with the dead and recurred or progressed status in the high-risk group (Figures 2F,G). Thereafter, GC patients were randomly separated into two parts (1:1) for training and testing sets. Table 2 summarized the clinical characteristics of GC patients from training, testing, and entire sets. Our data demonstrated that high-risk patients presented more undesirable OS than low-risk patients in training (Figure 2H), testing (Figure 2I), and entire sets (Figure 2J). ROCs at 5-, 6- and 7-year survival confirmed that the TNFα-derived risk score was accurately and sensitively predictive of GC prognosis in training (Figure 2K), testing (Figure 2L), and entire sets (Figure 2M).

Time-independent ROCs revealed that the TNFα-derived risk score displayed a prominent advantage in predicting GC prognosis (Figure 3A). Univariate cox regression analyses showed that the TNFα-derived risk score was indicative of an undesirable prognosis of GC (Figure 3B). Furthermore, multivariate cox regression analyses uncovered that the TNFα-derived risk score acted as an independent risk factor of GC outcomes (Figure 3C). Compared with G1/2, a higher risk score was detected in G3/4 patients (Figure 3D). Additionally, we noted that compared with stage I, there was a prominently increased risk score in stages II, III, and IV (Figure 3E). In comparison to the T1 stage, a significantly higher risk score was investigated in T2, T3, and T4 stages (Figure 3F). Also, N1 and N3 patients presented an increased risk score compared to those with N0 (Figure 3G). The aforementioned findings demonstrated that the TNFα-derived risk score was in relation to GC progression. Further survival analyses suggested that high-risk patients indicated poorer DFS (Figure 3H), DSS (Figure 3I), and PFS (Figure 3J) than low-risk patients. The clinical applicability of this model was further validated in the GSE15459 dataset. Consistently, the high-risk score predicted poorer OS than the low-risk score (Figure 3K). Additionally, ROCs at 3-, 4- and 5-year survival demonstrated that this model enabled the prediction of GC prognosis accurately and sensitively (Figure 3L).

A prognostic nomogram was established by integrating the TNFα-derived gene signature, age, T stage, N stage, M stage, and pathological stage, which might be predictive of BC patients’ survival outcomes through a quantitative scoring method (Figure 4A). In accordance with the nomogram, each patient would obtain a total point from each prognostic indicator. Calibration curves demonstrated that the predictive accuracy of this nomogram was similar to the actual OS outcomes (Figure 4B). With the median value of the nomogram score, GC patients were clustered into high- and low-risk groups. It was found that high-risk patients were indicative of more undesirable OS outcomes than low-risk patients (Figure 4C). ROCs at 5-, 6-, and 7-year OS demonstrated that this nomogram displayed excellent efficacy in the prediction of OS outcomes (Figure 4D). Additionally, decision curve analyses demonstrated that the nomogram had a remarkable advantage of the TNFα-derived gene signature alone and displayed a high potential for clinical utility (Figure 4E).

It was found that the epithelial-mesenchymal transition, UV response up, and Notch signaling presented enhanced activities in high-risk specimens, in accordance with GSEA results (Figure 5A). Meanwhile, protein secretion, G2M checkpoint, and mitotic spindle exhibited reduced activities in low-risk specimens. Moreover, we quantified the activities of ssGSEA gene sets in each GC specimen (Figure 5B). Compared with the low-risk group, graft rejection, angiogenesis, apical junction, complement, epithelial-mesenchymal transition, IL6-JAK-STAT3 signaling, inflammatory response, interferon-gamma response, and KRAS signaling showed remarkedly enhanced activities in the high-risk group (Figures 5C,D). Oppositely, late estrogen response, glycolysis, heme metabolism, MYC targets V2, p53 pathway, protein secretion, unfolded protein response, and UV response up had prominently reduced activities in the high-risk group.

Furthermore, analyses were presented for investigation of the difference in responses to small molecular agents between groups. Our study noted that CHIR.99021 and CI.1040 displayed higher IC50 values in high- than low-risk patients (Figure 6A). Additionally, high-risk patients showed reduced IC50 values of pazopanib, VX.702, PF.562271, FTI.277, TW.37, bosutinib, AZD8055, docetaxel, AZD6482, rapamycin, and DMOG in comparison to low-risk patients. The aforementioned data suggested that low-risk patients presented higher sensitivity to CHIR.99021 and CI.1040, while high-risk patients displayed enhanced responses to pazopanib, VX.702, PF.562271, FTI.277, TW.37, bosutinib, AZD8055, docetaxel, AZD6482, rapamycin, and DMOG. We also compared the differences in somatic mutations between high- and low-risk groups. The first ten mutated genes included TTN, TP53, MUC16, LRP1B, SYNE1, CSMD3, ARID1A, FLG, PCLO, and FAT4. Higher mutational frequencies of the aforementioned genes were observed in low- than high-risk groups (Figure 6B). Both in high- and low-risk groups, missense mutation was the most frequent mutational type (Figures 6C–E). In particular, there was a significant difference in TTN mutation between groups (Figure 6F).

Through the ESTIMATE algorithm, we estimated the infiltration levels of immune and stromal cells. As a result, the high-risk score was in relation to increased immune and stromal scores, as well as the ESTIMATE score (Figures 7A–C). The abundance levels of immune cells were quantified within GC tissues by the ssGSEA method. There were enhanced abundance levels of activated B cells, activated CD4 T cells, activated CD8 T cells, activated dendritic cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, eosinophils, gamma delta T cells, immature B cells, immature dendritic cells, macrophages, mast cells, MDSCs, memory B cells, natural killer T cells, neutrophils, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, type 1 helper cells, and type 2 helper cells in high- compared with low-risk specimens (Figures 7D–F).

Further analyses uncovered that immune checkpoint molecules containing HAVCR2, CD209, LAG3, SIRPA, TNFRSF4, CD274, CD28, CD27, CD96, TIGIT, and ICOS displayed enhanced expressions in high- compared with low-risk groups (Figure 8A). Additionally, the TNFα-derived risk score was positively associated with most immune checkpoint molecules (Figure 8B). We calculated the TMB score across GC tissues, with a median value of 2.1/MB (Figure 8C). A higher TMB score was investigated in low- than high-risk patients (Figure 8D). Moreover, we noted that high-risk patients presented elevated TIDE scores (Figure 8E). Nevertheless, no prominent difference in the IPS score was noted between groups (Figure 8F). The activities of all steps within the cancer immunity cycle were estimated in GC tissues (Figures 8G,H). In particular, there were reduced activities of the release of cancer cell antigens and enhanced activities of cancer antigen presentation in high- than low-risk groups (Figure 8I). The aforementioned data were indicative that the TNFα-derived genomic model might be applied as a predictor of immune responses in GC.

We noted that AKR1B1, CPVL, and CTSL within the TNFα-derived gene signature presented remarkably increased expressions in GC than in normal tissues (Figures 9A–C). Their expressions were further verified in human gastric mucosal cells GES-1 and human GC cells HGC-27, MKN-28, and AGS. Our data confirmed the significant upregulation of AKR1B1, CPVL, and CTSL in HGC-27, MKN-28, and AGS cells than GES-1 cells (Figures 9D–J).