In silico development and validation of a novel glucose and lipid metabolism-related gene signature in gastric cancer.

Background: Abnormal glucose and lipid metabolism plays a critical role in gastric carcinogenesis and development. Hence, we presented a systematic analysis of glucose and lipid metabolism-related genes to explore their function and prognostic value in gastric cancer (GC).
Methods: The consensus clustering algorithm was used to identify the molecular subtypes based on glucose and lipid metabolism-related genes. Subsequently, cox regression analysis and lasso regression analysis were utilized to establish a risk prediction model. A clinical nomogram was constructed to assist prognosis assessment. In addition, ESTIMATE and single-sample gene set enrichment analysis (ssGSEA) algorithms were performed to evaluate the immune infiltration of the metabolic model, and GSEA was used for enrichment analysis of the metabolic signature. Finally, we explored the association between the risk model and anti-cancer therapy for the purpose of clinical application for GC treatment.
Results: GC samples were divided into 2 subtypes based on glucose and lipid metabolism-related genes, patients in cluster 2 had a better overall survival (OS) than those in cluster 1. Fifty-two genes were identified by univariable regression analysis. Finally, a 13-gene metabolic signature (CACNA1H, CHST1, IGFBP3, NASP, STC1, VCAN, NUP205, NUP43, PGM2L1, CAV1, ELOVL4, PRKAA2, TNFAIP8L3) was successfully constructed that demonstrated good performance in different datasets, as well as an independent hazardous factor for prognosis. In addition, the nomogram constructed with the clinical variables showed higher predictive efficacy for predicting the 1-, 3-, and 5-year OS. The 13-gene metabolic signature was significantly associated with immune scores and immune cell infiltration in high-risk group. Moreover, GSEA analysis revealed that cancer- and immune-related pathways were enriched in the high-risk group. Finally, our results indicated that there might exist an immunosuppressive status in the high-risk groups. Conclusions: This study demonstrated that glucose and lipid metabolism-related genes were significantly associated with prognosis. Meanwhile, it will provide novel insights into exploring the immunoregulation roles of these genes. 2022 Translational Cancer Research. All rights reserved.

Introduction

Gastric cancer (GC) is among the most common and lethal digestive malignancy worldwide, according for 768,793 deaths in 2020 (1). Most patients with GC are frequently detected at later stages due to atypical and unobvious symptoms. Despite great progress in GC treatment, the 5-year survival rate for patients with advanced GC still need to be improved (2,3). The high recurrence and metastasis rate lead to shorter survival (4). Therefore, there is an urgent need to explore the mechanism of tumorigenesis and identify prognostic biomarkers of GC.

Emerging evidence has confirmed that there is a systemic link between tumor growth and metabolic pathways. Aberrant metabolic reprogramming, especially glycolysis and lipid metabolism, is currently considered a hallmark feature of cancers, including GC (5-7). Glucose metabolic reprogramming maintains the acquisition of energy allowing cancer cells to survive during different disease states (8,9). Glucose metabolism is also highly correlated to lipid metabolism (10). The lipid metabolic reprogramming refers to the sufficient supply of energy and the synthesis of structural and functional lipids, thereby satisfying the demands of increased membrane biogenesis in cancer cells (11,12). The dysfunctional metabolic programming based on glucose and lipid metabolism has also been proven to contribute to tumor metastasis, chemotherapy drug resistance, and immune escape (13). However, glucose and lipid metabolism have a high level of complexity in its regulation. No such risk model has been established based on glucose and lipid metabolism related genes. Moreover, the potential mechanism in relation to disease is still unknown. Several coding or/and noncoding RNAs have been shown to mediate a wide range of biological processes in cancer cells, including glucose and lipid metabolism (14,15). These genes exert direct or indirect effects on glucose and lipid metabolism both in vivo and in vitro.

To enhance our understanding of glucose and lipid metabolism in GC, it is necessary to perform an integrative analysis of the potential roles of glucose and lipid metabolism related genes. In this study, we identified the metabolic subtypes based on the expression of genes related to glycolysis and lipid metabolism. A GC prognostic signature was established and verified based on glucose and lipid metabolism-related genes. In addition, we constructed a nomogram based on the prognostic signature and clinical features to improve the markers’ clinical utility. Finally, we also clarified the relationship between prognostic signature and other in patients with GC, including the signaling pathway and tumor microenvironment. This study might provide research clues for the accurate prevention and treatment of GC. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-168/rc).

Methods

Data source

The RNA sequencing data and the corresponding clinical information of GC patients were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) dataset (Table S1), respectively. Subsequently, the Fragments per Kilobase Million (FPKM) data were normalized to transcripts per kilobase million (TPM) values. The expression values in GSE84437 and GSE62254 dataset [the Asian Cancer Research Group (ACRG) cohort] were performed to log2-transformed. Four samples without clinical follow-up information were removed. TCGA-Stomach Adenocarcinoma (STAD) (371 samples) and GSE84437 (433 samples) data sets were merged to produce the final training set (804 samples), and GSE62254 (300 samples) served as the final validation set. The R package SVA was utilized to perform batch effect removal on the different data sets. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The flowchart is shown in .

Figure 1

The flowchart of this study. ACRG, Asian Cancer Research Group; GEO, Gene Expression Omnibus; GSEA, gene set enrichment analysis; ICI, immune checkpoint inhibitor; STAD, Stomach Adenocarcinoma; TCGA, The Cancer Genome Atlas.

Glucose and lipid metabolism gene set

Four glucose metabolism datasets [Kyoto Encyclopedia of Genes and Genomes (KEGG) glycolysis gluconeogenesis, Hallmark glycolysis, Reactome regulation of glycolysis, and Reactome glycolysis] and four lipid metabolism datasets (KEGG glycerophospholipid metabolism, hallmark fatty acid metabolism, Reactome metabolism of lipids, and Reactome phospholipid metabolism) were obtained from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb) and chosen as reference gene sets. A total of 1,118 glucose and lipid metabolism related genes were selected for the study (Tables S2,S3).

Using the consensus clustering algorithm to identify molecular typing based on glucose and lipid metabolism-related genes

Cluster analysis algorithms were gradually utilized to explore hidden groupings, and exhibited a good performance. In this study, a consensus clustering analysis was utilized to cluster samples of training set using the R package ConsensusClusterPlus. The parameters were shown below: number of repetitions =1,000, number of clusters from 2 to10, pFeature =1, and pItem =0.8. According to the inflection point, as illustrated in , the optimal number of clusters was further determined. To evaluate the prognostic implication of glucose and lipid metabolism-related genes, the Kaplan-Meier curve was selected to compare the overall survival (OS) of the different subgroups.

Figure 2

Cluster analysis based on glucose and lipid metabolism-related gene and OS in the cluster 1/2 subgroups. (A-C) The clustering results when the number of classifications is k=2, 3, and 4. (D) Consensus clustering CDF with k valued 2 to 9. (E) Relative alteration in the area under the CDF curve with k valued 2 to 9. (F) Kaplan-Meier survival curves for clusters C1 and C2 in the training set (P=0.003). C1 and C2 represent different metabolic statuses. GC in C1 is marked with red, and C2 is marked with blue. CDF, cumulative distribution function; GC, gastric cancer; OS, overall survival.

Lasso cox regression analysis

Lasso regression analysis was conducted to perform variable selection, which could compress variables to overfit the risk mode (16,17). Certain coefficients were compressed via setting others to zero to determine the optimal factor (16,17). Ten-fold cross-validation were used running lasso regression to avoid overfitting. Subsequently, these factors were utilized to establish an optimized model based on multivariate Cox regression analysis. Patients with GC were segregated into high- and low-risk groups using the median risk score value. Kaplan-Meier method was used to assess the difference in survival outcomes between two groups. The receiver operating characteristic (ROC) curve was plotted to evaluate the performance of the risk model. The risk model was evaluated using a validation set, which was consistent with its coefficient from training set.

Construction of a nomogram

Univariate and multivariate Cox regression analyses were performed to choose potential risk factors. The clinical variable was added to the nomogram model for further analysis. The nomogram was utilized to predict the 1-, 3-, and 5-year OS rate with the R package rms. The calibration and ROC curve were applied to estimate the performance of the nomogram model.

Gene set enrichment analysis (GSEA)

To elucidate relevant biological significance of risk score, GSEA analysis was utilized to expound on the noteworthy differences between the high- and low-risk groups. We performed 1,000 repetitions of gene set permutations, and utilized the normalized enrichment score (NES) as evaluation metrics. The thresholds were set as false discovery rate (FDR) <0.05.

Evaluation of immune scores and immune infiltration

To evaluate the immune status of different samples, we performed ESTIMATE and single-sample GSEA (ssGSEA) algorithms to explore differences between the different groups. The ESTIMATE algorithms were used to analyze the immunological characteristics of the high- and low-risk groups. In addition, we performed ssGSEA to analyze the immune infiltration statuses and relevant immune-related pathways based on the R package gsva.

Exploration of the potential susceptibility to anti-cancer therapy

To explore the potential effect of anti-cancer therapy, we also evaluated the correlation between the expression of key immune checkpoint genes. Besides, we performed correlation analysis based on Pearson method using R software.

Statistical analysis

All the statistical analyses were conducted with R software (Version 3.6.3). P<0.05 or FDR <0.05 was served as statistically significant.

Results

Cluster analysis based on glucose and lipid metabolism-related genes

The consensus clustering algorithm was used to cluster the training set samples (804 samples) through the R package ConsensusClusterPlus. Due to higher inter-group correlation when k=4 not k=2, the optimal number of clusters selected was 2 (). Namely, the GC samples in the training set could be clustered into subtypes, named cluster 1 and cluster 2, based on glucose and lipid metabolism-related genes. Besides, the survival analysis was further utilized to evaluate the prognostic value of subcluster. As shown in , patients in cluster 2 showed more favorable OS rates than patients of cluster 1 (P=0.003) (). The result showed that glucose and lipid metabolism-related genes are associated with prognosis in GC.

Identification of glucose and lipid metabolism-related prognostic genes

Next, univariate Cox analysis was employed to identify the potential prognostic value of these glucose and lipid metabolism-related genes in the training set. P value less than 0.01 was selected as the threshold for filtering, and 52 genes were significantly associated with OS. These genes were selected for further analysis.

Construction and validation of the prognostic risk model in different datasets

The R package glmnet was utilized to perform Lasso-penalized Cox regression analysis for 52 genes with potential prognostic value based on training set. The optimal lambda value was associated with the determinants of penalized term. As shown in Figure S1, we obtained the theoretically optimal model when lambda =13. Therefore, 13 genes were eventually selected for detailed model building. The coefficient of risk score formula was shown in . To assess the model’s robustness, the Asian Cancer Research Group (ACRG) cohort was introduced as the independent validation datasets. Kaplan-Meier curve was further drawn to show the survival status of the GC samples (). This study suggested that high RiskScore sample has a worse prognosis. The same coefficients as that from the training set were utilized to obtain the risk score of each patient according to the expression level. The survival status of the GC samples in the high-risk group was still worse than that of in the low-risk group (). Subsequently, we evaluated 1-, 3-, and 5-year predictive capability based on ROC. As shown in , the model has a good area under curve (AUC) value, and the AUC values calculated from TCGA for 1, 3, and 5 years were 0.611, 0.619, and 0.641, respectively. Besides, the validation results indicated that the AUC value of the glucose and lipid metabolism-related prognostic signatures was 0.635 in 1 year, 0.644 in 3 years, and 0.634 in 5 years (). Moreover, we calculated the risk score based on the expression level of these prognostic genes, and the samples were divided into the high- and low-risk groups based on the median risk score (). A high proportion of deaths with GC were still shown in the high-risk group, which was consistent with the training set (). The above results showed that our model has good robustness, and could support a stable potential for a predictive prognostic signature.

Table 1

The gene coefficients of risk model

Gene	Coefficient
CACNA1H	0.0073492701064191
CHST1	0.0481537527082623
IGFBP3	0.0350611211681865
NASP	−0.0014241982646828
STC1	0.0830987998258511
VCAN	0.0209038567717716
NUP205	−0.076027779878957
NUP43	−0.108261343313484
PGM2L1	0.10183699405434
CAV1	0.0559587258963077
ELOVL4	0.0506717462447253
PRKAA2	0.0994553688728891
TNFAIP8L3	0.00733755211794015

Figure 3

Prognostic value of the risk model in the training and validation set. (A,B) Kaplan-Meier survival analysis of GC patients between the high- and low-risk group based on the risk score formula in the training and validation set. (C,D) ROC curves of the 1-, 3-, and 5-year OS based on the risk signature in the training and validation set. (E,F) The median value of risk scores with survival and statuses of GC patients depends on the risk signatures. GC in high-risk group is marked with red, and in low-risk group is marked with blue. AUC, area under the curve; GC, gastric cancer; OS, overall survival; ROC, receiver operating characteristic curve.

The relationship between risk score and clinical variables

Univariate and multivariate Cox regression analyses were selected for evaluation of the independent prognostic predictor successively. We demonstrated that the age [P<0.001, hazard ratio (HR) =1.026, 95% confidence interval (CI): 1.016–1.036], T (P=0.001, HR =1.255, 95% CI: 1.093−1.442), N (P<0.001, HR =1.549, 95% CI: 1.383–1.735), and risk score (P<0.001, HR =1.241, 95% CI: 1.162−1.326) were significantly correlated with OS in the training set (). After adjusting for confounding factors, the age (P<0.001, HR =1.029, 95% CI: 1.019−1.039), N (P<0.001, HR =1.467, 95% CI: 1.305−1.649), and risk score (P<0.001, HR =1.223, 95% CI: 1.143−1.309) were confirmed to be independent prognostic factors by multivariate Cox regression analysis (). Similarly, univariate and multivariate Cox regression analyses showed the age, N, and risk score were independent prognostic factors as those in external dataset ().

Figure 4

Identification of the risk score as an independent prognostic factor by Cox regression analysis. (A,B) Training set. (C,D) Validation set.

Establishment of a nomogram model to further improve the predictive power

To better improve the predictive potential of GC patients, a nomogram was then constructed to facilitate prediction of OS. The clinical variable, including age, gender, T, and N, was added to the nomogram model. As shown in , the calibration curve for the nomogram is close to the 45° line, which indicated a good performance of the nomogram. Figure S2 showed that the nomogram model presents good accuracy in predicting OS, with AUC values of 0.715, 0.736, and 0.743 at 1, 3, and 5 years, respectively, in the training set. Besides, the nomogram model was shown with the AUC values at 1, 3 and, 5 years were 0.791, 0.761, and 0.723, respectively.

Figure 5

Development and validation of a prognostic nomogram based on risk signature. (A,B) Construction of a nomogram for 1-, 3-, and 5-year OS prediction along with risk score and clinical variables. (C,D) Calibration curves of OS for the GC patients. **, P<0.01; ***, P<0.001. GC, gastric cancer; OS, overall survival; Pr, predicted.

Evaluation the degree of correlation between risk score and immune score

ESTIMATE algorithm was used to calculate each GC sample’s immune and matrix scores based on R package estimate. As shown in , there were significant differences in the StromalScore, ImmuneScore, and ESTIMATEScore between the two groups in the training set (). Besides, these significant differences were still existed in the independent validation ACRG cohort (). In both cohort, patients in the high-risk group had a higher StromalScore, ImmuneScore, and ESTIMATEScore. All these data demonstrated that the risk score constructed by glucose and lipid metabolism-related genes might correlate with immune response.

Figure 6

Correlation between the different risk signature based on ESTIMATE algorithm. The high-risk group shows a significantly elevated immune score (A), Stromal score (B) and ESTIMATE score (C) in the training set. The results in the validation set were consistent with that of in the training set (D-F). GC in high-risk group is marked with red, and in low-risk group is marked with blue. GC, gastric cancer.

Exploration of relation between the risk score and pathways

To elucidate the potential links between the different samples’ risk scores and underlying biological significance, we performed GSEA analysis on the two risk-score groups. The results showed that several cancer-related pathways and immune-related pathways were significantly enriched in the high-risk group. Among them, the top five pathways in the high-risk group, include KEGG_CALCIUM_SIGNALING_PATHWAY, KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION, KEGG_PATHWAYS_IN_CANCER and KEGG_AMINOACYL_TRNA_BIOSYNTHESIS, KEGG_GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM, KEGG_HOMOLOGOUS_RECOMBINATION, KEGG_MISMATCH_REPAIR, KEGG_ONE_CARBON_POOL_BY_FOLATE in the low-risk group ().

Figure 7

Gene set enrichment analysis based on risk scores. (A) The result of high-risk group in the training set. (B) The result of low-risk group in the training set. (C) The result of high-risk group in the validation set. (D) The result of low-risk group in the validation set. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Three pathways (KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION) in the high-risk group and two pathways (KEGG_AMINOACYL_TRNA_BIOSYNTHESIS, KEGG_GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM) in the low-risk group were validated by the ACRG cohort, respectively (). These results suggested that these pathways’ imbalance was related closely to GC development.

Correlation between the risk score and tumor immune infiltration

To further explore the correlation between the risk score and tumor immune infiltration, ssGSEA algorithm was applied to find the differences in infiltration of primary immune cells in GC samples. We revealed that several types of tumor-infiltrating immune cells were enriched in the high-risk group, including dendritic cells (DCs), macrophages, mast cells, neutrophils, and tumor-infiltrating lymphocyte (TIL). These results demonstrated good consistency between the training and validation set (). Besides, we could differentiate between the high- and low-risk groups in both cohorts based on contents of the antigen presentation process, including chemokine-chemokine receptor (CCR), major histocompatibility complex (MHC) class I, parainflammation, and type II interferon (IFN) response (). Interestingly, the risk score was negatively correlated to MHC class I, suggesting there was a potential trend in the high-risk group toward enhanced immunosuppression compared with the low-risk group. The result was further confirmed by the findings of immune checkpoint inhibitors (ICIs).

Figure 8

The relationship between the risk score and immune status. (A,B) Box plot of differences in immune cell infiltration between two groups in the training and validation set, respectively. (C,D) Box plot of differences in immune function between two groups in the training and validation set, respectively. *, P<0.05; **, P<0.01; ***, P<0.001. aDCs, activated dendritic cells; APC, antigen-presenting cell; CCR, chemokine-chemokine receptor; DCs, dendritic cells; HLA, human leukocyte antigen; iDCs, inhibited dendritic cells; IFN, interferon; MHC, major histocompatibility complex; NK, natural killer; ns, not significant; Th, T helper; TIL, tumor-infiltrating lymphocyte.

Evaluation of the association between the risk model and anti-cancer therapy

To evaluate the potential clinical application for GC treatment, we attempted to investigate the potential susceptibility to ICIs. We firstly assessed the correlation between risk scores and several genes affecting ICIs. Our study revealed that the risk score was related to the expressions of FAP, TAGLN, LOXL2, whereas it was negatively associated with POLE2, FEN1, MCM6, MSH6, MSH2 (). Furthermore, we investigated the potential susceptibility to ICIs based on the expression level of the immune regulator. The results indicated that the expressions of CD200, CD44, TNFRSF4, NRP1, CD276, CD48, and CD28 were elevated in the high-risk groups, whereas the expression of HHLA2 was downregulated (). Our results indicated that there might exist an immunosuppressive status in the high-risk groups.

Figure 9

Analysis of the relationship of the ICIs and/or DNA in different group. (A) Correlation analysis between risk score and ICIs and/or DNA in the training set (asterisk indicates a statistically significant). (B) Correlation analysis between risk score and ICIs and/or DNA in the validation set (asterisk indicates a statistically significant). (C) Comparison of the ICIs between the high and low-risk group in the training set. (D) Comparison of the ICIs between the high and low-risk group in the validation set. *, P<0.05; **, P<0.01; ***, P<0.001. ICIs, immune checkpoint inhibitors.

Discussion

Several macromolecules, such as proteins, nucleic acids, and lipids, have been proven to contribute to cell growth and proliferation via the production of metabolic intermediates (18,19). These agents act as the precursors necessary for life, which divert biological media from biomass production to improve organismal or cell viability (18,19). Accumulating evidence indicates that abnormal metabolic reprogramming is critically involved in the oncogenic signaling pathways, and has emerged as a hallmark of cancer (20-22). The results suggest metabolic markers might be possible therapeutic targets. As increasingly high-throughput techniques emerge, the novel technical means provide valuable molecular information to understand the potential roles of glucose and lipid metabolism. Previous studies demonstrated that a variety of biological characteristics could alter metabolic phenotype during a period of the transformation of normal cells into malignant cells, including enhanced glucose uptake and reactivation of de novo lipid biosynthesis (20-22). In addition, glucose provides energy and precursors for biosynthetic processes, and lipids function as an important component of some pathways and cell-cell communication (23-25). There is tight connection between oncogenic transformation and abnormal metabolic reprogramming. However, no specific glucose and lipid metabolism related model has been developed for GC to predict prognosis and treatment response until now. Therefore, we constructed a risk model based on glucose and lipid metabolism gene signatures that could act as prognostic and therapeutic biomarker in GC.

Previous study confirmed that GC has been identified as a glycolytic-enhanced malignancy, and metabolomic studies have shown that there exists a predominance of lipid storage over consumption in gastric epithelial neoplasms (26-29). In this study, we have applied the consensus clustering algorithm to the data set, and identified two subtypes of GC first based on genes related to glucose and lipid metabolism. Further, univariable regression analysis was performed to screen out OS-related gene sets, and 52 genes were obtained to display the cancer-related functions. Finally, the 13-gene risk signature was constructed by lasso penalized regression and multivariate Cox analysis, and demonstrated that the signature had high reliability and stable predictability. The application of Nomogram improved the predictive accuracy of the model notably, which could provide support to clinical decision-making.

The gene signature was constructed with CACNA1H, CHST1, IGFBP3, NASP, STC1, VCAN, NUP205, NUP43, PGM2L1, CAV1, ELOVL4, PRKAA2, and TNFAIP8L3, which was correlated significantly with the OS and tumor immune microenvironment (TIME). Wang et al. reported that stanniocalcin-1 (STC1) positively regulated Bcl-2 to mediate gastric carcinogenesis, metastasis and chemoresistance under hypoxia, indicating that the carcinogenic and therapeutic roles of STC1 in GC (30). Wang et al. found that Caveolin‑1 could enhance RANKLinduced multi-cancer cells migration, including stomach, lung, renal and breast cancer cells (31). Rao et al. indicated that the signaling pathway STK11-PRKAA2-ULK1 was linked to Gastrin-mediated autophagy activation in gastric adenocarcinoma cells (32). Jimenez et al. demonstrated that CACNA1H mediates the antiproliferative effects and promotes cancer stem cell suppression in hepatocellular carcinoma following exposure to tumor-specific electromagnetic fields (33). Tu et al. found that cancer-related gene vasohibin-2 was critically involved in the proliferation of breast cancer cells by activating insulinlike growth factors, including insulin like growth factor binding protein 3 (IGFBP3) and insulin like growth factor binding protein 6 (IGFBP6) (34). Kang et al. indicated that nuclear autoantigenic sperm protein (NASP) knockdown could facilitate transcription release and failure of replication initiation in hepatocellular carcinoma, which was achieved by abolished the supply of histone H3 and enhanced chromatin accessibility (35). Several bioinformatic analyses were performed to VCAN in GC without further exploration of the mechanisms using in vivo/in vitro experiments (36-38). Tian et al. found that NUP43 might predict prognosis in luminal A and in HER2+ breast tumors based on bioinformatics analysis (39). However, the potential roles of CHST1, NUP205, PGM2L1, ELOVL4, and TNFAIP8L3 remain unclear. Their roles need to be classified based on further validation.

Increasing evidence has confirmed that the TIME was strongly associated with the disease development and progression. In this study, the relationship between the risk scores and TIME in GC was under comprehensive evaluation. We demonstrated that GC patients with high-risk scores had higher StromalScore, ImmuneScore, and ESTIMATEScore than those with low-risk scores. Previous study has found that GC patients with high stromal scores had a better prognosis, which was consistent with the present study (40). In-depth analysis of immune infiltrates indicated that the risk signature was correlated positively with the expression of DCs, Macrophages, mast cells, neutrophils, and TIL. Additionally, the risk score was negatively correlated to MHC class I. MHC class I plays the role of a gatekeeper, which promotes pathogens and transformed cells recognition via displaying peptides to CD8+ T cells (41). For the tumor to arise and progress, the evolutionary mechanism will be dominant to avoid elimination by CD8+ T cells and hence the immune response (41). Down-regulation of the MHC class I expression at the cell surface was considered one of the hallmarks of immune escape, which significantly impaired the ability of CD8+ T cells to recognize the cancer cells (42). The present study indicated that the down-regulation of the MHC class I expression was more evident in our high-risk group, suggesting there might exist the reduction of tumor-associated antigenic epitopes to facilitate immune evasion. In addition, among the immune checkpoint markers, the expressions of CD200, CD44, TNFRSF4, NRP1, CD276, CD48, and CD28 were highly expressed in the high-risk groups, whereas the expression of HHLA2 was downregulated. The result indicated that GC patients grouped by this risk signature showed differences in responsiveness to immunotherapy. CD200 and CD44 are stem cell-specific markers that play an important role in immunosuppression, which are highly expressed in several cancers and are correlated with unfavorable prognosis (43,44). Blocking CD200 or CD44 inhibit immune activation, and contribute to improve the efficacy of immunotherapy (45,46). TNFRSF4, also known as OX40L, is expressed on regulatory T cells, which promotes immune escape of leukemia stem cells (47). The high expression of NRP1 is associated with the poor prognosis of several cancers, which may assist cancer immune escape and enhance tumor progression via mediating Treg cell infiltration (48). CD276, also known as B7-H3, has been detected in several malignancies and associated with tumor progression and poor outcome (49,50). Agonistic targeting of CD276 can enhance immune surveillance and promote the anti-tumor response (51,52). CD28 is a T cells costimulatory factors, which transduces a positive signal to enhance the proliferation of T cells (53). Besides, CD28 ligation can also promote the production of various cytokines, such as IL-2, IL-4 and IL-10 (53). However, the role of CD28 in Treg cell differentiation seems to be determined by some other conditions, such as the level of T cell receptor (TCR) engagement and cytokine environment. Although HHLA2 functions as the inhibitor of CD4+ and CD8+ T cells and is helpful for cancer immunotherapy (54), the potential roles in GC deserve further study.

Although we performed this study based on a large sample of data sets, there were still several limitations about this study. First, this is a retrospective analysis, and all the data were obtained from public databases. Our risk signature must be further validated by prospective cohort. Secondly, high-risk factors for GC, such as unhealthy diets (particularly diets high in salt) and infections with Helicobacter pylori, were not included in this study. These factors might also affect GC patients’ prognosis. Thirdly, the AUC value ranging from 0.6 to 0.7 is not very satisfactory. We constructed the nomogram model to improve it, and showed good performance. However, there still needs more accurate model later through developing novel algorithms. Finally, although GSEA is an accepted method for functional analysis, the potential mechanisms mediated by the glucose and lipid metabolism-related risk genes must be further explored comprehensively and in-depth.

Conclusions

In this study, we constructed a robust 13-gene signature prognostic risk model in different datasets that was better for prediction of prognosis and treatment response. It will provide novel insights into exploring the immunoregulation roles of these genes.

The article’s supplementary files as