Identification and Development of Subtypes with Poor Prognosis in Gastric Cancer Based on Both Hypoxia and Immune Cell Infiltration

Introduction

Gastric cancer is the fourth most common incident cancer and the third most common cause of cancer death.¹ In addition, previous literature suggests that gastric cancer is a heterogeneous disease featuring many distinct histological and molecular subtypes.² As an example, Lauren classification, proposed in 1965, has been recognized by both clinicians and pathologists and continues to be used. However, the limitations of Lauren classification as a prognostic predictor are also gradually unveiled. Currently, whether intestinal type gastric cancer patients have better prognosis in comparison with diffuse-type patients remains controversial.^3,4 Even more, Lauren classification could not be used as an independent prognostic factor on large population-based studies, calling into question their reliability.^5,6

The exploration of gastric cancer subtypes used for diagnosis and treatment has never stopped in recent years. In 2014, according to key DNA defects and molecular abnormalities, The Cancer Genome Atlas Consortium (TCGA) group described four subtypes for gastric cancer including Epstein–Barr virus (EBV)-positive, microsatellite instability (MSI), gene stable (GS), and chromosome instability (CIN) types.⁷ Moreover, Oh et al indicated two different transcriptomic subtypes in gastric cancer, MP and EP, where the MP subtype was closely linked to significantly poor survival and resistance to chemotherapy.⁸ So an in-depth exploration and development of tumor subtypes of gastric cancer is imperative to improve therapeutic effects and prognosis.

Hypoxia is an important and common feature of solid tumors.⁹ The occurrence of this condition is generally considered to be caused by the unrestricted proliferation of the tumor and the hypertonicity and tortuosity of the vascular network in the tumor.^10,11 The disturbance of immune cell infiltration in tumors caused by hypoxia has been widely discussed, such as myeloid-derived suppressor cells (MDSCs) and tumor-associated macrophages (TAMs) which are important components in TME. Tumor-associated macrophages (TAMs) derived from Macrophages were found to preferentially appear in hypoxic areas of tumors.¹² In tumor-bearing mice, a quick increase of PD-L1 is observed in MDSCs, macrophages, dendritic cells, and tumor cells under hypoxia, and this phenomenon depends on an oxygen-regulated subunit HIF-1α. At the same time, Noman et al also found that T cell activation mediated by MDSCs can be significantly enhanced when PD-L1 is blocked under hypoxic conditions, which is evidence for the simultaneous inhibition of PD-L1 and HIF-1α can be utilized as novel immunotherapy.¹³

In previous literatures, it has received heightened attention that the heterogeneity of hypoxia or tumor immune cell infiltration can be used to guide the therapeutic strategy and predict the prognosis in different kinds of cancers.^14–16 However, most of the studies were based on only hypoxia or tumor immune cell infiltration. The research on identification of tumor subtypes by combining the characteristic of both hypoxia and immune cell infiltration is regrettably scarce.

In this study, we tried to utilize a large sample size of gastric cancer transcriptome data to describe the characteristics of tumor hypoxia and immune-related cell infiltration. Combining these two characteristics, we attempt to develop hypoxia-immune classifiers to identify gastric cancer subtypes with different biological features and outcome.

Materials and Methods

Retrieval, Downloading, and Processing of Gene Expression Profiles in Patients with Gastric Cancer

Cohorts with patient survival information or relatively complete clinical information were retrieved from the Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) database. We selected several qualified cohorts and downloaded all their expression profiles and corresponding clinical information (including GSE15459, GSE34942, GSE57303, GSE62254/ACRG, GSE84437, GSE26253, and TCGA-STAD) (Supplementary Table 1).^17,18 The above data were all downloaded from the website https://portal.gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/. In the TCGA database, we also obtained somatic mutation data of patients with gastric cancer.

We converted high-throughput sequencing data (FPKM) from TCGA-STAD to transcripts per kilobase million (TPM) values.¹⁹ R software and the “ComBat” algorithm are used to remove batch effects and then we merge five of the cohorts (GSE15459, GSE34942, GSE57303, GSE62254 (ACRG), GSE84437) into one GC (gastric cancer) metadata for further analysis (n=1059).²⁰ We also downloaded the transcriptome data and clinical information of GSE13911, including the patient’s microsatellite status information.

We download GSE114083 as another validated data set, which contains mRNA and lncRNA data of gastric cancer cell lines treated with hypoxia (n=4) and normoxia (n=4).

Unsupervised Clustering for Hypoxic and Immune-Related Characteristic Gene Sets

We extracted a total of 23 different immune cell marker gene sets from Charoentong P’s research,²¹ and we searched and selected 13 gene sets that up-regulated with Hypoxic in the Molecular Signatures Database.^22,23 These 36 gene sets in total are shown in Supplementary Table 2.

We capitalized on the ssGSEA (single-sample gene-set enrichment analysis) algorithm in the “GSVA” R packages to measure the enrichment scores of these 36 gene sets in each patient of the GC metadata.²⁴ These enrichment scores were considered to represent the infiltration of immune cells and the intensity of hypoxia in tumors. By integrating the enrichment scores of these 36 gene sets, we put the Consensus ClusterPlus package into service to complete unsupervised clustering analysis and identified different hypoxic-immune patterns in gastric cancer patients. The above algorithm ran repeatedly 1000 times to ensure the stability of the results.²⁵

ESTIMATE algorithm was used to predict tumor purity, and the presence of infiltrating stromal/immune cells in tumor tissues in different clusters, and output stromal score, immune score, and estimate score (Supplementary Table 3).²⁶ In addition to the ssgsea algorithm, we also used CIBERSORT analysis to calculate the 22 types of immune cell proportions in the tumor, and select samples with P<0.05 for the next step to verify (n=974) the accuracy of the hypoxic-immune patterns and further analysis.²⁷

Determination of Differential Genes Between Hypoxic-Immune Patterns and Construction of Hypoxic-Immune Gene Signature

Based on differences in immune cell infiltration and hypoxia pathway enrichment, we divide gastric cancer patients into three hypoxic-immune patterns. To obtain differentially expressed genes (DEGs) related to different hypoxic-immune types, limma R package was utilized in every two patterns. We set the screening criteria to “FDR < 0.05”.²⁸ A web-based portal, Metascape, was utilized to conduct biological functional enrichment analysis of these DEGs and we considered p < 0.01 statistically significant.²⁹

Univariate Cox proportional hazard regression analysis was utilized to extract prognosis-related genes in ACRG samples, with p < 0.05 as the criteria.

With the above prognostic-related genes, we performed unsupervised clustering in the ARCG data set and we classified patients into several clusters for the next analysis. We used R package, glmnet, to perform the least absolute shrinkage and selection operator (Lasso) Cox regression analysis, and selected the potential prognostic genes to construct a prognostic model. According to the median score of this model (HIscore), gastric cancer patients were divided into high-HIscore groups and low-HIscore groups.

The Coef_i means the coefficients, and x_i is the TPM of each gene included in the model.

Identification of candidate small-molecule drugs for gastric cancer patients with high HIscore.

The CMAP database (https://portals.broadinstitute.org/cmap/), composed of 7056 gene expression profiles induced by 1309 small molecules, is widely applied to explore the potential unknown roles of existing drugs on diseases.³⁰

The R package “limma” was utilized to identify DEGs in the high and low HIscore samples in ACRG and TCGA cohorts, respectively. The criteria for DEG screening were as follows: false discovery rate (FDR) <0.01. The top quarter DEGs with largest |log2 (fold change)| were extracted from the up-regulated and down-regulated DEGs, respectively. We take the intersection of The top quarter up-regulated and down-regulated DEGs in ACRG and TCGA, respectively, and the result obtained is used as the input file of the CMAP database. Then, we obtained the potential small drug molecules, and P < 0.05 was regarded as the cut-off criteria.

Acquisition of Biological Processes

We obtained 21 biological processes from previous studies and the KEGG database,³¹ shown in Supplementary Table 4, including EMT1, EMT2, angiogenesis, antigen processing and presentation, antigen processing machinery, CD8 T effector, immune checkpoint, base excision repair, cell cycle, DNA damage repair, DNA replication, mismatch repair, ECM–receptor interaction, JAK-STAT signaling pathway, MAPK signaling pathway, NF-kappa B signaling pathway, nucleotide excision repair, Pan-F-TBRS, PI3K-Akt signaling pathway, TGF-beta signaling pathway, Wnt signaling pathway. We not only calculated and compared their enrichment scores in different patient clusters but also obtained the correlation between them and hypoxic-immune scores (HIscore) through Pearson’s correlation analysis (Supplementary Table 5).³²

Tumor Mutation Burden, Microsatellite Status, Genomic Instability

We downloaded the somatic mutation data of gastric cancer patients in TCGA and calculated the number of mutation events/million bases as tumor mutation burden (TMB). Non-coding alterations were excluded from TMB calculation.³³ Microsatellite instability information was downloaded from The Cancer Immunome Atlas (https://tcia.at/) including MSS (microsatellite stability, n=229), MSI-H (microsatellite instability-high, n=59), MSI-L (microsatellite instability-low, n=49). In the TCGA cohort, we also developed a framework for evaluating genome stability. Sorting the patients according to the accumulated number of somatic mutation genes from high to low, the first 25% were considered to be the group genome unstable (GU), and the last 25% were considered to be the group genome stable (GS).³⁴

Gene Set Enrichment Analysis

GSEA was performed between hypoxia and normoxia gastric cancer cell lines in GSE114083. The number of random sample permutations was set at 1000, and NOM p-value <0.05, FDR q-value <0.25, and | NES | > 1 were set as the significance threshold.³⁵

Statistical Analysis

Correlations were assessed using Pearson’s coefficients. A Kruskal–Wallis test was used in the comparison between three or more groups and a Wilcox test was used in the comparison between two groups. We used the median value of HIscore in ACRG as the cutoff value, and with the cut-off, we divided the gastric cancer patients into the corresponding data set into high and low HIscore groups. Kaplan–Meier method was adopted to generate survival curves, Log rank tests were utilized to determine whether there is a significant difference in prognosis between high and low groups.³⁶ The hazard ratios (HR) of genes related to patient classification are obtained from univariate Cox regression analysis. Whether it is an independent prognostic factor is judged by multivariable Cox regression analysis.³⁷ Multivariable Cox regression analysis was only used for gastric cancer patient data sets with detailed clinical information (ACRG and TCGA-STAD). The results of multivariable Cox regression analysis were visualized by the forest plot. The R package timeROC was used to generate the receiver operating characteristic, and the AUC was also calculated to evaluate the specificity and sensitivity of the HIscore. A combined nomogram was established as a quantitative tool for predicting the likelihood to die of each patient using the “regplot” R package. The concordance index (C-index) was calculated to assess the consistency between model prediction and actual clinical outcomes of patients. The calibration plot was established to evaluate the accuracy of the prediction for 1-, 3-, and 5-year overall survival using this nomogram by the “rms” R package.

The somatic mutation landscape between high and low HIscore groups was obtained from TCGA-STAD, and the R package maftools were used for visualization.³⁸ In all the above statistical processes, p < 0.05 means statistical significance. We utilized R software for data processing and analysis, version 4.0.4.

Results

Hypoxic-Immune Infiltration Patterns

We merged GSE15459, GSE34942, GSE57303, GSE62254, and GSE84437, and then we obtained a large sample size meta-cohort with complete patient survival information. With the R package Consensus ClusterPlus, we divided patients in meta-cohort into different subtypes based on the enrichment scores of 23 immune cell-related gene sets and 13 hypoxia-related gene sets. Based on such unsupervised two-dimensional hierarchical clustering, we identified three distinct HIclusters (hypoxic-immune cluster) with different immune infiltration and hypoxic states (Figure S1A), including 416 cases in HIcluster A, 372 cases in HIcluster B, and 271 cases in HIcluster C (Figure 1A). Kaplan–Meier analysis pointed out that there are significant differences in OS among the three subgroups (Figure S1B). The median survival time of patients in HIcluster C was 10.30 years, and the median survival time of HIcluster A and HIcluster B was 7.13 years and 4.80 years, respectively. HIcluster C has the best 5-year survival rate of 51.7%, while the 5-year survival rates of HIcluster A and HIcluster B are 45.5% and 39.8%, respectively. This implies that gastric cancer patients with low hypoxia and high immune characteristics may have a better prognosis.

Figure 1 Gastric cancer subtypes with different immune cell infiltration and hypoxia status. (A) Unsupervised clustering of 36 gene sets related to hypoxia or Immune cell infiltration in GC metacohort including the 5 independent stomach cancer cohorts. Each column represented patients and each row represented the enrichment scores for these gene sets. The clinical information collated serves as the patient’s annotation. Differences in stromal scores (B), immune score (C) and estimate score (D) between different HIcluster. The Kruskal–Wallis test was used to compare the statistical difference between HIcluster. (E) The component differences of immune cells among the three HIcluster analyzed by CIBERSORT. Differences in PD-L1 (F) and CTLA4 expression (G) among three HIclusters in GC metadata. The Kruskal–Wallis test was used to compare the statistical difference. (***p < 0.001.).

The expression level of gene HIF1A in HIcluster B was significantly higher than HIcluster A and HIcluster C (Figure S1C). We scaled and compared the enrichment scores of 23 immune-related gene sets and 13 hypoxia-related gene sets in different HIclusters involved in unsupervised clustering (Figure S1D and E). The enrichment scores of all hypoxia-related gene sets in HIcluster B were significantly higher than HIcluster A and HIcluster C. The trend of the enrichment scores of CD56dim natural killer cell and neutrophil in HIclusters was similar to that of hypoxia-related gene sets, both of which were that HIcluster B is significantly higher than other groups. The enrichment scores of other immune cell infiltration-related gene sets were the highest in HIcluster C, the lowest in HIcluster A, and the middle level in HIcluster B. The stromal score obtained by the ESTIMATE algorithm was significantly lower in HIcluster A than the other two groups, while there was no difference in stromal scores between HIcluster B and C (Figure 1B). The immune score and estimate score were the lowest in HIcluster A and the highest in HIcluster C (Figure 1C and D). The differences between the three groups are statistically significant. By CIBERSORT analysis, 22 types of immune cell proportions in different clusters are shown in Figure 1E. The median expression of immune checkpoint genes PDL1 and CLTA4 was the lowest in HIcluster A, and that of CLTA4 was the highest in HIcluster C. There was no significant difference in the median expression of PDL1 between HIcluster B and C (Figure 1F and G).

Combining the above results, HIcluster A was classified as low hypoxic/low immune subtype, HIcluster B was classified as high hypoxic/high immune subtype and HIcluster C was classified as low hypoxic/high immune subtype.

Construction of Hypoxic-Immune Signatures and Functional Annotation

Because the ACRG cohort has more comprehensive clinical information on the 300 gastric cancer patients, we extracted the ACRG cohort from the meta-cohort to study the clinical and biological characteristics in the ACRG cohort that may be closely linked to different hypoxic-immune infiltration patterns. Similar to the clustering results of meta-cohort, ACRG cohort also found three different hypoxic-immune infiltration patterns (Figure S2A and B). Through the limma R package, we obtained 5607 differentially expressed genes between HIcluster A and B, 7770 differentially expressed genes between HIcluster A and C,2776 differentially expressed genes between HIcluster B and C. Taking the intersection of these genes, we obtained a total of 880 DEGs that were significantly differentially expressed in the three clusters (Figure S2C). A web-based portal, Metascape, was utilized to conduct biological functional enrichment analysis of the DEGs. When P < 0.01 is considered statistically significant, we obtained a total of 1907 enriched biological processes (Supplementary Table 6). The top 20 enrichment results with the smallest p-value are shown in Figure S2D. Many results were related to activation of an immune response, immune cell recruitment, such as GO:0046649 (lymphocyte activation), GO:0002366 (leukocyte activation involved in immune response), and GO:0050900 (leukocyte migration). The largest number of DEGs (133) was enriched in lymphocyte activation, and its p-value was the smallest (Log10 P=−60.86).

With univariate Cox regression analysis, we selected 293 OS-related genes from the ACRG cohort of 880 DEGs (p<0.05), and they are shown in Supplementary Table 7.

Based on these 293 genes, unsupervised clustering analyses were performed in the ARCG cohort to classify patients into different subtypes. Similar to the results of clustering based on immune infiltration and hypoxia gene sets, we obtained three different genomic phenotypes (HIgenecluster A–C) (Figure 2A and B; Supplementary Table 8). It implied that three different hypoxia-immune infiltration patterns do exist. The median survival time of HIgenecluster B–C was 7.453 and 1.970 years, and the 5-year survival rates of HIgenecluster A–C were 56%, 58.5%, 28.8%, respectively (Figure 2C).

Figure 2 Identify different genomic subtypes in ACRG cohort. (A) Consensus matrices of the ACRG cohort for k = 3. (B) Unsupervised clustering of overlapping hypoxia-immune-related genes in ACRG cohorts to classify patients into different genomic subtypes (genecluster A–C). The geneclusters, HIclusters, molecular subtypes, stage, gender, survival status and age were utilized as patient annotations. (C) Survival analyses for the three geneclusters in ACRG cohort including 132 cases in gencluster A, 116 cases in genecluster B, and 52 cases in genecluster C. Kaplan–Meier curves with Log rank p value< 0.001 showed a significant survival difference among three geneclusters. Genecluster C showed significantly worse OS than the other two genclusters.

In the Lauren classification of gastric cancer obtained in the ACRG cohort, we found that the diffuse subtype occupies the largest proportion in genecluster C, and the intestinal subtype occupies the largest proportion in genecluster A. In geneclusters A–C, the proportion of diffuse subtypes gradually increased, and the proportion of intestinal subtypes gradually decreased (Table 1, Figure 3A). Interestingly, tumor patients with EMT molecular subtypes have an absolute advantage in HIgenecluster C, while genecluster C does not include MSI molecular subtype, tumor patients. In genecluster, the proportion of patients with molecular subtypes MSS/TP53+ or MSS/TP53- is less than that of geneclusters A and B. The molecular subtypes of gene clusters A and B are the same (Table 1, Figure 3B). In addition, we also compared the proportion of patients with EMT subtypes in different tumor stages and T stages with HIgeneclusters, and found that they were much smaller than that in the HIgenecluster C group (Figure 3C and D).

Table 1 Clinicopathological Characteristics in Relation to Different HIgeneclusters in ARCG Cohorts

Figure 3 The proportion of Lauren-classifications (A) and ACRG molecular subtypes (B) in the three geneclusters. The proportion of ACRG molecular subtypes in different T stages (C) and tumor stages (D). (E) Difference in the expression of 21 known biological processes related to tumor progression among three geneclusters. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and black dots showed outliers. The asterisks represented the statistical p value (****P < 0.0001).

Characteristics of Transcriptome Traits in Different Hypoxic-Immune-Related Genecluster (HIgenecluster)

We have collated and obtained 21 biological processes from previous studies and the KEGG database. The sources of these biological processes, the genes involved in these processes, and the enrichment scores of each sample in ACRG were shown in (Figure 3E).

The processes EMT1 and EMT2 have the highest enrichment scores in HIgenecluster C, which was mutually corroborated the fact that the EMT subtype (molecular subtypes in ACRG cohort) in HIgenecluster C accounts for a large proportion. And the enrichment scores of EMT1 and EMT2 in HIgeneclusterB were higher than that of HIgenecluster A. Immune-related biological processes, such as Antigen processing and presentation, CD8 T effector, etc were enriched more in HIgenecluster B. Surprisingly, for some pathways that related to genome stability, such as Nucleotide excision repair, base excision repair, DNA damage repair, DNA replication Mismatch repair, etc, all revealed a low enrichment in HIgenecluster C. Other signal pathways closely related to the occurrence and development of tumors also had significant differences between different geneclusters.

Based on the enrichment results of these biological processes in different HIgeneclusters, we classified HIgenecluster B as an immune-activation group, and HIgenecluster C classified as EMT-Genomic stability group. Combined with classification and previous research, the poor prognosis of gastric cancer patients in HIgenecluster C has been explained to a certain extent.

Construction of the Hypoxic-Immune-Related Signature

Based on the previously obtained 293 Hypoxic-immune infiltrations phenotype-related DEGs, unsupervised clustering allowed us to identify a part of patients with poor prognosis (HIgenecluster C). In order to dig deeper into their prognostic value, we constructed a hypoxic-immune-related signature. With ACRG cohort as the training set, the Lasso Cox Regression Analysis was utilized to construct the hypoxic-immune-related signature, and 13 hypoxic-immune-related signatures were built for predicting the overall survival of cular sugastric carcinoma patients (Figure S3A and D). Based on the TPM value and related coefficients of these 13 signatures, we obtained the risk scores (HI score) of each patient.

HIscore = (BEX4*0.060280828757125) + (BTN3A2*0.0480478080183311) + (CCDC50*0.221803823749138) + (CLIP4*0.0185094741736464) + (COL4A1*0.0245555436398619+CPE*0.0137755467630004) + (HLA.DRB6*-0.0326550077240013) + (JAZF1*0.0282045392726293) + (LMO2*0.0299067833201967)+(PSMB10* - 0.120416260002512) + (RASGRP2*0.0162557186287523) + (MXD1*-0.007241171408135) + (ZNF326* - 0.114149705889599).

According to the median value (1.61406348) of HI-score in the ACRG cohort, patients were divided into high HIscore group (n=150) and low HIscore group (n=150). Kaplan–Meier survival curves depicted the specifically remarkable survival advantage in the high HIscore group (P<0.0001) (Figure S3B). The median survival time of the high HIscore group (2.3556 years) was worse than that of the low HIscore group (>7.133 years). The 5-year survival rates of the low HIscore group were better than that of the high HIscore group (66.6% VS 36.5%). ROC curves indicated that HI-score ability to predict OS is extraordinary in the ACRG cohort (3-year AUC = 0.72, 5-year AUC = 0.73, 8-year AUC = 0.74; Figure S3C).

We utilized the same formula to obtain the HIscore of each gastric cancer patient in the TCGA cohort and made the median value of the HIscore (1.61406348) in the ACRG cohort as the cut-off to classify the patients into a high HIscore group (n=70) and a low HIscore group (n=267). ACRG cohort was utilized as a validation set to test the robustness and stability of hypoxic-immune-related signatures. The K–M survival curves attested to the higher HIscore were related to the better OS (Figure S3E). The median survival time of the high HIscore group (1.663 years) was worse than that of the low HIscore group (3.1589 years). The 3-year survival rates of the low HIscore group were better than that of the high HIscore group (53.1% VS 31.9%). The ROC analysis indicated that the AUC of predicting 5-year OS in the TCGA cohort was 0.7 (Figure S3F).

Including the factors of age, sex, TNM status, tumor stage, and grade (only TCGA cohort), multivariate Cox regression analysis was utilized to prove that HIscore can be an independent prognostic biomarker for gastric cancer patients in both ACRG (HR 3.710 (2.659–5.178); Figure S3G) or TCGA (HR 2.026 (1.163–3.529); Figure S3H).

To establish a quantitative tool for predicting the prognosis of gastric cancer patients in the ACRG cohort, a nomogram was generated by integrating clinicopathological factors and HIscore based on the multivariable Cox proportional hazards model (Figure S3I). The point scale in the nomogram was used to generate point to these variables, and the risk of death of each patient was evaluated by accumulating total points of all variables. The HIscore had the most excellent weight among all these variables, which was consistent with the result of the previous multivariable Cox regression analysis. The calibration plots confirmed the significant consistency between predicted and observed actual clinical outcomes of gastric cancer patients (Figure S3J).

Characteristics of Transcriptome Traits in High and Low HIscore Groups

We evaluated the 23 TME immune cells infiltrating tumors in the high- and low-risk groups, and found that most of the TME immune cells have significant differences between the two groups. activated CD4 T cell, activated dendritic cell and neutrophil had higher enrichment scores in the low HIscore group, while other immune cell enrichment scores had higher enrichment scores in the high HIscore group (Figure 4A). The expression levels of some genes related to tumor hypoxia status were also significantly differentially expressed in high and low HIscore groups. Among the hypoxia genes with statistically significant differences in the Kruskal–Wallis test, except for VKORC1 and HIF1A, which is highly expressed in the high HIscore group, the others were highly expressed in the low HIscore group (Figure 4B).

Figure 4 Immune cell infiltration characteristics and transcriptome traits in different group with high or low HIscore. (A) The abundance of each TME infiltrating cell in high and low HIscore groups. (B) Difference in the expression of 18 hypoxia gene signatures between high and low HIscore groups. The asterisks represented the statistical p value (Wilcox test, *P < 0.05; **P < 0.01; ***P < 0.001). (C) Differences in HIscore among three geneclusters in ACRG cohort. The Kruskal–Wallis test was used to compare the statistical difference between three geneclusters (P < 0.0001). (D) Differences in HIscore between different ACRG molecular subtypes. The Kruskal–Wallis test was used to compare the statistical difference between four ACRG molecular subtypes (p < 0.0001). (E) Correlations between HIscore and the known gene signatures in ACRG cohort using Pearson's analysis. Negative correlation was marked with blue and positive correlation with red. (F) Difference in the expression of 21 known biological processes related to tumor progression between high and low HIscore groups. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and black dots showed outliers. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001).

The vast majority of patients with EMT subtypes were classified as HIgenecluster C, and all patients in HIgenecluster C carried high HIscore. Kruskal–Wallis test indicates that HIscore has a significant difference between HIgenecluster A–C. HIgenecluster A had the lowest median HIscore and HIgenecluster C had the highest median HIscore (Figure 4C), which indicates that HIscore may serve as a signature of EMT-Genomic stability. More clearly in ACRG cohort, the median HIscore of EMT mole type was significantly higher than other.

The median HIscore of MSI molecular subtypes was higher than that of MSS subtypes, and whether TP53 is mutated will not significantly affect the HIscore of MSS subtypes (Figure 4D).

To legitimately decipher the characteristics of HIscore, we calculated the correlation between HIscore and the signatures of some biological processes (Figure 4E; Supplementary Table 5), and the enrichment scores of these biological processes between the high and low HIscore groups were shown in Figure 4F. Pan-F-TBRS was the pathway was the most positively related to HIscore (R = 0.706, P < 0.001), and TGF-beta signaling pathway is also closely related to HIscore (R = 0.629, P < 0.001). Echoing the results of the previous alluvial diagram, two signatures related to epithelial-mesenchymal transition also showed a strong positive correlation with HIscore (EMT1, R = 0.593, P < 0.001; EMT2, R = 0.550, P < 0.001) and their median enrichment scores were higher in the high HIscore group. It was striking that the biological process signatures related to the maintenance of genome stability are strongly negatively correlated with HIscore, such as base excision repair (R = −0.739, P < 0.001), DNA replication (R = −0.659, P < 0.001), cell cycle (R = −0.630, P < 0.001), DNA damage repair (R = −0.577, P < 0.001), nucleotide excision repair (R =−0.568, P < 0.001) and mismatch repair (R = −0.564, P < 0.001). The median enrichment scores of genome stability-related biological processes are lower in the high HIscore group.

The Role of HIscore in Assessing Genome Instability

In the TCGA cohort, based on the prognostic model we constructed, we obtained the HIscore of all gastric cancer patients and utilized 1.61406348 (the median value of HI-score in ACRG cohort) as the cut-off to divide the patients into high and low HIscore groups. Pearson correlation analysis indicated that a high HIscore could be closely linked to high TMB (R =−0.46, P < 0.001; Figure 5A). We divided the patients into high or low TMB groups based on the median value of TMB and the high TMB group had a lower median HIscore (Figure 5B). We found that both high and low HIscore groups accounted for a large proportion of patients with MSS, but the low HIscore group accounted for the higher proportion of MSI-L (16% VS 9%) and MSI-H patients (21% VS 6%) than that in high HIscore group (Figure 5D). Kruskal–Wallis test revealed that patients with different Micro satellite status in the TCGA cohort had significantly different median HIscore (Figure 5E). MSS patients had the highest median HIscore, and MSI-H patients had the lowest median HIscore. Taking the upper quartile and lower quartile of the total number of somatic mutation genes as cut-off, we obtained the genome unstable (GU) group (n=81) and the group genome stable (GS) group (n=87) in the TCGA cohort, and we extracted these 168 patients and combined their HIscore. In the low HIscore group, genome unstable patients accounted for a larger proportion (Figure 5G). Not only that the median HIscore in the genome unstable group was lower than that in the genome stable group (P<0.0001, Figure 5H). We also evaluated the expression levels of all genes that constitute the prognostic model in different TMB, Micro Satellite, and genome stability states (Figure 5C, F, and I). In the TCGA cohort, the high TMB, the MSI, and the genome unstable (GU) group all showed low expressions of BEX4, CCDC50, CLIP4, CPE, JAZF1, LMO2, and RASGRP2, while PSMB10 was highly expressed in these groups.

Figure 5 The close relationship between HIscore and TMB, microsatellite status, genomic instability. (A) Correlation between the TMB and HIscore in the TCGA cohort. (Pearson correlation analysis, R=−0.46, P < 0.0001) (B) Difference in the HIscore of patients between high and low TMB group. (C) Difference in the expression of 13 members of the hypoxia-immune signature between high and low TMB group. (D) The proportion of MSS, MSI-L and MSI-H patients in the high and low HIscore groups. (E) Difference in the HIscore of patients among MSS, MSI-L and MSI-H groups. (F) Difference in the expression of 13 members of the hypoxia-immune signature among MSS, MSI-L and MSI-H groups. (G) The proportion of genome unstable and genome stable patients in the high and low HIscore groups. (GU: genome unstable patient; GS: genome stable patient) (H) Difference in the HIscore of patients between GU and GS group. (I) Difference in the expression of 13 members of the hypoxia-immune signature between GU and GS group. The Kruskal–Wallis test was used to compare the statistical difference (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

In order to verify our inference that HIscore is related to genome instability, we also evaluated the expression of these common differentially expressed genes in another GSE13911 containing Micro satellite information of gastric cancer patients. Because the Micro satellite status in GSE13911 only included MSS (n=20) and MSI (n=20), the difference of these genes between the MSI-L and the MSI-H group could not be observed. Wilcox test revealed that except for LMO2, PSMB10, and RASGRP2, the median gene expression levels of other genes in the MSI group were significantly lower than those in the MSS group (P<0.05), which was consistent with the phenomenon we found in the TCGA cohort (Figure S4). The median expression of LMO2 was lower in the MSI group, although the P-value was not less than 0.05 (Wilcox test, P=0.07).

Identification of Candidate Small Molecule Drugs for Gastric Cancer Patients with High HIscore

In total, we identified 546 up-regulated genes and 337 down-regulated genes in high HIscore groups of both ACRG and TCGA cohort (Supplementary Table 9). By using the CMAP database we obtained total 96 potential small drug molecules, and P < 0.05 was regarded as the cut-off criteria (Table 2). The enrichment scores (− 1 to + 1) that could assess the similarity between genes and drugs were calculated. Specifically, enrichment score <0 revealed molecules had potential antagonistic effects to biological processes closely linked to high HIscore, indicating they could reverse the poor prognosis and could serve as therapeutic drugs, especially for gastric cancer patients with high HIscore. In our research, the enrichment scores of adiphenine (−0.939), isoflupredone (−0.929), vinblastine (−0.916), and viomycin (−0.907) are all lower than −0.9. They are the four small drug molecules with the enrichment score closest to −1. This showed that they have great potential in reversing the poor prognosis and serving as therapeutic drugs.

Table 2 Identification of Candidate Small-Molecule Drugs for Gastric Cancer Patients with High HIscore Using the CMAP Database

Verify the Accuracy of the Hypoxia Pathway Gene Set Included in the Study

To determine the correctness of the 13 hypoxia-related gene sets involved in unsupervised clustering in this study, we performed GSEA between hypoxia and normoxia gastric cancer cell lines in GSE114083. The genes of gastric cancer cell lines in the hypoxic group were significantly enriched in 13 hypoxia-related gene sets we used (Figure S5 and Supplementary Table 10).

Discussion

Previous researches have forced us to get to the deeper issue of the biological and molecular subtypes of gastric cancer. However, relevant classifiers are still in great demand for guiding clinical treatment or predicting prognosis. Hypoxia, a kind of tumor microenvironmental stress related to poor prognosis of patients, is a hallmark of solid tumors,³⁹ and immune cell infiltration characterizations are also a research hotspot in many solid tumors.⁴⁰ Most classifications focus on a single biological and molecular feature and in this study, a new classifier was developed by combining the hypoxic state and immune cell infiltration of gastric cancer.^41,42

Due to the complexity of immune activity and hypoxia in the tumor micro environment, few transcriptome-based biomarkers could evaluate hypoxia and immune status, so we took advantage of ssGSEA to calculate the enrichment scores of related immune and hypoxia gene sets as our evaluation index.^43,44 Of course, the gene sets of these hypoxia-related processes have been verified in GSE114083 and the results show that they can represent tumor hypoxia traits in gastric cancer cell lines. In the GC (gastric cancer) meta-cohort after removing the batch effect, the unsupervised clustering algorithm revealed that there are three markedly different gastric cancer subtypes. Among the three subtypes, HIcluster A was classified as low hypoxic/low immune, HIcluster B was classified as high hypoxic/high immune, and HIcluster C was classified as low hypoxic/high immune. We utilized the expression of HIF1A to verify the correctness of our judgment of the three types of hypoxia, while the ImmuneScore, the expression of PDL1, and CLTA4 were used to verify the correctness of the immune cell infiltration type.

Firstly, we found that the median OS of HIcluster with high hypoxia levels was smaller than that with low hypoxia levels. This phenomenon was in line with our expectation. Hypoxia leads to increased glycolysis and lactate production.⁴⁵ Glycolysis can reduce the production of reactive oxygen species (ROS), and the low concentrations of ROS can contribute to tumor cell survival.⁴⁶ In addition, lactate, generated during glycolysis, decreases the tumor environmental pH and an acidic pH distinctly impedes the function of normal immune cells, such as T-cell and tumor-infiltrating lymphocytes.⁴⁷ Furthermore, the lactate-mediated enhancement of tumor cell motility is reported not only in single-cell motion but also in enforced bulk migration.⁴⁸ Hyaluronan, synthesized by tumor-associated fibroblasts (TAF), increases around the carcinoma regions at high lactate concentrations and encourages the growth and motility of cancer cells.⁴⁹ Moreover, hypoxia induces EMT and the formation of angiogenesis and lymphangiogenesis, which will allow tumor cells to escape the adverse microenvironment and disseminate into secondary sites.^50,51 Taken together, these findings may at least partially explain the facilitation of tumor progression and the worse prognosis in the high hypoxia levels cluster as compared to the low hypoxia levels cluster.

Even for the same tumor subtype in which immune cells highly infiltrated (HIcluster B and C), different hypoxia status still caused different infiltration patterns. For example, activated CD8+ T cells, activated B cells and immature B cells have a higher infiltration in HIcluster C. This phenomenon is consistent with the research of Akio Ohta et al, which indicated that lower tissue oxygenation decrease the extent of T cell activation.⁵² Previous studies have also pointed out that the narrow maturation process of dendritic cells under lower local O2 tension may also reduce T cell activation.⁵³ HIFs (HIF-1, HIF-2, HIF-3) are involved in most biological processes caused by hypoxia.⁵⁴ Burrows et al generated a mouse with B cell-specific deficiency of VHL (an important protein that causes the dissociation of HIFs and indicated that HIFs activation results in immature B cell developmental arrest and peripheral B cell loss.⁵⁵ The infiltration pattern of B cell especially activated B cell in HIcluster B and C also conforms to this rule.

We compared the median survival time between HIcluster B and C. Both of them are highly immune cells infiltrated, however, HIcluster C had a better prognosis. We speculate that it is not only due to hypoxic effects but also closely related to the different infiltration characteristics of activated CD8+ T cells and B cells in tumors. Abnormal angiogenesis, desmoplasia, and inflammation will all be promoted due to the oxygen supply obstacles in the tumor that proliferates unrestrictedly, which can cause tumor progression and drug resistance.^56,57 Furthermore, activated HIF1a up-regulation of the activity of Snail and Twist, two transcription factors that decrease E-cadherin expression and facilitate EMT.⁵⁸ Previous studies have already shown that pancreatic cancer and bladder urothelial carcinoma patients with high densities of CD8+ T-cells have a better prognosis.⁵⁹ Kim et al constructed B-cell depleted mice of which the tumor size is larger and the growth rate is higher than the control group. They further researched and found that the B-cell-specific gene expression is a prognostic biomarker of head-neck carcinoma.⁶⁰ Therefore, we speculate that hypoxia, CD8+ T cells, and B cells may also play a similar role in gastric cancer, resulting in the poor prognosis in HIcluster B. Immune checkpoint inhibitors therapy has developed rapidly in recent years, such as that target CTLA4 and PD1 pathway, but only a fraction of patients can benefit from it.⁶¹ The median expression of PDL1 and CLTA4 was the lowest in HIcluster A, and that of CLTA4 was the highest in HIcluster C. There was no significant difference in the median expression of PDL1 between HIcluster B and C. Plentiful researches in different cancer types have indicated a positive correlation between PDL1 expression and ICI response and OS.^62,63 The new classifier we created has the potential in guiding ICI treatment for patients with gastric cancer.

Furthermore, We explored the differences in transcriptome traits between different HIclusters, and performed the enrichment analysis with the differentially expressed genes obtained. These biological processes are mainly concentrated on the activation and differentiation of immune cells, leukocyte migration, which is consistent with the characteristics of our classifier. In order to amplify the ability of the classifier, we created to classify the poor prognosis patients with gastric cancer. We selected genes with the prognostic value among the differential genes to perform unsupervised clustering again. We were surprised to isolate a genecluster C with a significantly poor prognosis. We evaluated the interaction between other gastric cancer classification methods in ACRG cohort and genecluster, and we found that EMT subtype and diffuse subtype have an absolute advantage in genecluster C. Diffuse-type gastric cancer, a subtype with a poor prognosis, lacks effective and specialized treatments and accounts for 30% of the total number of gastric cancer patients.^64,65 As the early stage of tumor metastasis, EMT includes not only the dissolution of cell–cell junctions but also loss of apico-basolateral polarity. Throughout the EMT, down-regulation of cellular adhesion proteins and up-regulation of mesenchymal markers occur on gastric cancer cells.⁶⁶ Clinically, EMT is closely linked to a poor prognosis.⁶⁷ Therefore, we suppose that The epithelial-mesenchymal transition (EMT) is a pivotal characteristic of tumors in genecluster C, and it is also the basis for our classifier to classify patients with worse prognosis. In order to verify our ideas, we measured the expression levels of EMT-related genes in different geneclusters, and found that the EMT-related gene set (EMT1, EMT2) obtained in the previous literature got the highest enrichments score in genecluster C. Similar to EMT1 and EMT2, other pathways and biological processes closely related to the EMT processes are also enriched in genecluster C, such as ECM–receptor interaction,⁶⁸ MAPK signaling pathway,⁶⁹ PI3K-Akt signaling pathway,⁷⁰ TGF-beta signaling pathway and Wnt signaling pathway.^71,72 This supports our conjecture.

The pathways related to repairing DNA damage in genecluster A and B are more activated than genecluster C, such as nucleotide excision repair, base excision repair, and mismatch repair.^73,74 This implies that geneclusters A and B (especially genecluster A) will undergo more processes of detecting and removing bases, repairing modified bases, and removal of bulky adducts in DNA. Each of these DNA repair pathways resects a damaged region and inserts new bases to fill the gap. It implies that the stability of genomes in geneclusters A and B may have already been challenged.^75,76 Cell cycle and DNA replication have a negative correlation with genome stability, so their lower enrichment in genecluster C supports our guess that genecluster C belongs to the subtype of genome stability.^77,78

Furthermore, we utilized lasso Cox regression analysis to screen out 13 differentially expressed genes among HIclusters to construct a hypoxic-immune-related signature (HIscore). According to the median value of HIscore, patients are divided into high- and low-risk groups. Gastric cancer patients in the high-HIscore group have a poorer prognosis. We utilized TCGA cohort as the validation set and found that HIscore is still robust in predicting prognosis. Multivariate Cox regression analysis indicated that whether in ACRG or TCGA cohort, HIscore was an independent risk factor for overall survival (ACRG (HR 3.710 (2.659–5.178))); TCGA (HR 2.026 (1.163–3.529)). The accuracy of HIscore in ACRG cohort was evaluated by the AUC of receiver operating characteristic analysis and the AUC was 0.72 at 3 years, 0.73 at 5 years, and 0.74 at 8 years, respectively. The AUC was 0.70 at 5 years in the TCGA cohort. This suggests that HIscore is clinically significant in predicting the long-term prognosis of patients. In addition, we found that whether the patient undergoes adjuvant chemotherapy or not, HIscore’s ability to predict the prognosis will not be diminished, which could make HIscore more widely used in clinical practice. In addition, we developed a nomogram to describe the risk of death of each gastric cancer patient and the calibration plots confirmed the significant consistency between predicted and observed actual clinical outcomes. This nomogram is established by the combination of hypoxic-immune-related signature and clinical characteristics including TNM stage and Lauren classification. So that it showed considerable property for predicting the survival and can be more clinically applicable.

The vast majority of patients with EMT subtypes were classified as HIgenecluster C, and patients in HIgenecluster C carried higher HIscore. In ACRG cohort, the median HIscore of the EMT molecular subtype was significantly higher than other subtypes. This suggests that HIscore may be an indicator of the characteristics of genecluster C (EMT-genome stability type). Pearson correlation analysis showed that EMT-related biological processes are strongly positively correlated with HIscore, and DNA damage repair processes are strongly negatively correlated with HIscore. At the same time, the enrichment scores of EMT-related biological processes are higher in the high HIscore group, and the enrichment score of the DNA damage repair process was lower in the high HIscore group. Furthermore, we compared the HIscores between different TMB and microsatellite status, which is closely linked to genomic instability.^79,80 The results are consistent with our expectations that gastric cancer patients with high TMB and MSI have lower HIscore, and Pearson correlation analysis showed that HIscore is negatively correlated with TMB. We also ranked the total number of mutated genes from most to least. The first 25% were identified as genomically unstable patients, and the last 25% were identified as genomically stable patients.³⁴ In this new framework for evaluating the extent of genomic instability in gastric cancer, high HIscore is still a good indicator of genome stability.

We focused on the members that constitute hypoxic-immune-related signature (HIscore) and found that most of their expression levels are significantly different between genomic stability and instability group. BEX4 protein mainly localizes at centrosomes, contractile rings, and midbodies, and BEX4 contributes to the stabilization and elongation rates of microtubules. Moreover, BEX4 plays an important role in apoptosis and the adaptation to aneuploidy induced by spindle damage or a mitotic checkpoint defect.⁸¹ CLIP4, one of the CLIP-170 family, connects microtubule with organelle by CAP-Gly which contributed to cell polarity maintenance, intracellular transport, cell migration, and oncogenesis.^82,83 Multiple studies indicated that microtubule is altered in cancer cell divisions and linked to chromosomal instability, aneuploidy, and drug resistance.⁸⁴ Parvin et al found that LMO2 expression is sufficient to increase genomic instability in DLBCL and T-ALL cell lines, leading to the accumulation of chromosomal rearrangements and phosphorylated H2AX (g-H2AX), a marker of DNA double-strand breaks (DSBs).⁸⁵ Therefore, these genes may be the cornerstones of HIscore to predict the genome stability of gastric carcinoma.

In the TCGA cohort, BEX4, CCDC50, CLIP4, CPE, JAZF1, LMO2, and RASGRP2 all showed lower expressions level in the high TMB group, the MSI group, and the genome unstable (GU) group, while PSMB10 was highly expressed in all these groups. We chose another independent cohort (GSE13911) from the GEO database as the validation set. The expression trend of the majority in the validation set is similar to that in the TCGA cohort. The P-value of LMO2 is 0.07, which is close to 0.05, and we think it makes partial sense that the expression of LMO2 was lower in the MSI group.⁸⁶ The P-value of LMO2, PSMB10, and RASGRP2 is greater than 0.05, which may be caused by insufficient sample size in the verification set.⁸⁷ The results of the validation set further confirm that the HIscore we constructed is closely related to genome stability.

In short, in clinical practice, HIscore is not only an independent prognostic biomarker but also a forecast for clinicopathological features of patients with gastric cancer, such as EMT and genomic stability. This provides new insights into personalized treatment for some traditional chemotherapy drugs, such as cisplatin and gemcitabine that target DNA damage.^88,89 At the same time, because of the significant differences in immune cell infiltration and TMB between patients with high and low HIscore, it may also have the potential to screen out patients who will benefit from immunotherapy.⁹⁰ In addition, based on the fact that patients with high HIscore have a worse prognosis, we used the CMAP database to screen 96 small-molecule drugs that had potential antagonistic effects to biological processes closely linked to high HIscore. Specially, in our study, adiphenine, isoflupredone, vinblastine and viomycin have great potential in reversing the poor prognosis and serving as therapeutic drugs.

Our research provides new ideas and materials for the personalized clinical treatment plans for patients with gastric cancer, although some limitations of this study should be acknowledged. Firstly, our study only included a bioinformatics analysis, lacking the validation of solid clinical specimens. Additionally, the research was conducted with a retrospective design rather than a prospective one. However, our risk model was validated by multiple data cohorts; therefore, the result is still reliable and acceptable. Thus, future studies with prospective clinical trials and mechanistic exploration are warranted to further validate the present result.

Conclusion

In conclusion, our study evaluated the hypoxic state and immune cell infiltration in gastric tumors and identified different gastric cancer subtypes through the interaction of them. Subsequently, we obtained three subtypes with different prognosis: low hypoxic/low immune subtype, high hypoxic/high immune subtype, and low hypoxic/high immune subtype. Moreover, based on the DEGs between the above three subtypes, we can clearly identified another three sub-classification of patients with gastric cancer through unsupervised clustering including immune-activation subtype and EMT-Genomic stability subtype.

We extensively explored the differences in clinical and transcriptome traits among those cancer subtypes. Then, we constructed hypoxic-immune-related signature (HIscore) by LASSO regression. HIscore is an independent prognostic biomarker that could not be ignored to link to tumor EMT and genomic stability. Based on HIscore, we used the CMAP database to explore small-molecule drugs that may have the potential in serving as therapeutic drugs such as adiphenine, isoflupredone, vinblastine and viomycin. In addition, a nomogram was established by using the combination of HIscore and clinical characteristics to describe the risk of death of each gastric cancer patient.

This study proposes new insights on the molecular classification of gastric cancer and may contribute to the advancement of personalized treatment for gastric cancer. However, further experimental studies are needed to verify our proposal.