Ferroptosis Patterns Correlate with Immune Microenvironment Characterization in Gastric Cancer.

OBJECTIVE: We aimed to build a ferroptosis-based classifier to characterize the molecular features of gastric cancers (GC) and investigate the relationship between different ferroptosis patterns and GC tumor microenvironment (TME).
METHODS: Based on the genomic and clinical information from TCGA portal and GEO database, non-negative matrix factorization (NMF) was used to identify ferroptosis subtypes in GC patients. In order to estimate the ferroptosis levels, we established ferroptosis subtype score (FSS) to quantify ferroptosis patterns and ferroptosis potential index (FPI) by principal component analysis (PCA). The correlations of different ferroptosis patterns with TME cell-infiltrating characteristics (including immune cell infiltration, immune checkpoints expression levels, tumor mutational burden (TMB) and immunotherapy response) were systematically analyzed.
RESULTS: Two ferroptosis subtypes, C1 (with lower FSS) and C2 (with higher FSS), were determined. C2 displayed a significantly lower FPI than C1. Besides, C2 was associated with diffuse subtype while C1 with intestinal subtype. As for TME characteristics, C2 was in accordance with the immune-excluded phenotype as it showed more active immune and stromal activities but lower TMB, less probability of immunotherapy response and poorer prognosis. C1 was linked to immune-inflamed phenotype as it had lower stromal activities but increased neoantigen load, enhanced response to immunotherapy and relatively better prognosis.
CONCLUSION: The systematic assessment of ferroptosis patterns and ferroptosis levels presented in our study implied that ferroptosis serves as an important factor in the formation of TME, which may expand the understanding of TME and provide a novel perspective for the development of targeted immunotherapeutic strategies for GC patients.

Introduction

Gastric cancer (GC) is the fifth most common malignancy with a third leading reason of cancer-related mortality.1 A high rate of recurrence, metastasis, resistance to conventional chemotherapy and toxic side effects are the primary death causes of GC patients.2

Over recent years, immunotherapy has bloomed into the new therapeutic frontier of cancer treatment and emerged evident survival benefits in a variety of malignancies.3,4 Unfortunately, only a small fraction of patients can benefit from it in clinic.5 Therefore, it is urgently needed to discover promising biomarkers with the potential to identify GC subgroups that respond to immunotherapy.

Ferroptosis, a newly discovered form of regulating cell death, is driven by the accumulation of lipid peroxidation and lethal reactive oxygen species (ROS).6 Several studies have indicated a dual role of ferroptosis in antitumor immunity. For example, CD8+T cells responding to immune checkpoint blockade secrete interferon gamma to suppress tumors through the induction of ferroptosis.7,8 In contrast, tumor cells undergoing ferroptosis might act as arachidonic acid (AA) donors for the transcellular biosynthesis of eicosanoids and participate in the production of immunomodulatory AA metabolites, exerting an impact on antitumor immunity.9 Besides, the increased intratumor production of prostaglandin E2 (PGE2), a principal mediator of inflammation associated with cells undergoing ferroptosis,10 assists tumor escaping from immune surveillance.11,12 Collectively, the dual role of ferroptosis in antitumor immunity suggests the complexity of their interaction, which indispensably requires further exploration.

In this study, based on the genomic and clinical information from TCGA portal and GEO database, we identified two ferroptosis subtypes with distinct prognosis and ferroptosis levels in GC and revealed a remarkable correlation between the ferroptosis patterns and TME cell-infiltrating characteristics. To evaluate the ferroptosis level and quantify the ferroptosis patterns, we established the ferroptosis potential index (FPI) and ferroptosis subtype score (FSS), both of which were demonstrated as precise predictors for the immunotherapy response and clinical outcome.

Materials and Methods

The overview of study design is displayed in Figure 1.

Figure 1

Overview of study design.

Data Collection and Preprocessing

Gene expression data and corresponding clinical information of GSE62254,13 GSE15459,14 GSE3494215 and GSE5730316 datasets were retrieved from the Gene Expression Omnibus database (GEO) and consolidated as a GC meta-GEO cohort comprising 626 GC patients, 9 of which were removed due to lack of follow-up information or 0-day follow-up time. Above four GEO datasets were performed on the same microarray platform of [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. The batch effects between different datasets were reduced by the “combat” algorithm of the “sva” package within R.17 The gene expression and corresponding clinical information of stomach adenocarcinoma (STAD) samples of TCGA cohort were downloaded from the UCSC Xena browser ().18 After removing samples without follow-up information or with a 0-day follow-up duration, 350 STAD samples were collected. The clinical characteristics of the patients included are detailed in Table 1. In addition, FerrDb ()19 collected 259 ferroptosis-related genes including driver, suppressor and marker. The confidence levels of genes involved in ferroptosis were assigned to 4 degrees including validated, screened, predicted and deduced. The species involved included human, mice, rat and drosophila. To ensure the accuracy and stability of the model, 121 human-related and validated ferroptosis-related genes were obtained and provided in .

Table 1

Clinical Characteristic of the GC Patient Used in This Study

	TCGA	Meta-GEO
No. of patient	350	617
Age
<60	112	198
≥60	238	417
Unknown	0	2
Gender
Female	124	206
Male	226	411
Lauren type
Intestinal	73	295
Diffuse	59	261
Mixed	1	59
Unknown	217	2
Grade
G1	9	NA
G2	125	NA
G3	207	NA
GX	9	NA
Stage
I	46	72
II	110	149
III	145	228
IV	35	166
Unknown	14	2
Survival status
OS day (median)	475	956
Ending
Survival	204	307
Death	146	310

Consensus Clustering for Ferroptosis-Related Genes

Non-negative matrix factorization (NMF) involves in factorization into matrices with nonnegative entries that can identify the original component of the composite data and decompose the complex data intuitively.20 NMF clustering of 121 ferroptosis-related gene expression profiles was applied to categorize the patients in the meta-GEO and TCGA cohorts into distinct subtypes, individually. Using “NMF” package of R, the samples were decomposed into clusters for k = 2–6 against 50 random initializations of both the actual and a permuted gene expression matrix. The k value reaching the largest difference between cophenetic correlation coefficients calculated from the actual and permutated data was chosen. Principal component analysis (PCA) projected high-dimensional data into low-dimensional space and constructs a new linear combination with maximum variance to realize data dimension reduction. Subsequent PCA was built around the ferroptosis-related gene expressions to visualize and validate the subtype distribution.

Dimension Reduction and Generation of FPI

To evaluate the ferroptosis level, we adapted a method similar to Liu’s research21 to construct the ferroptosis potential index (FPI). We downloaded the gene expression profiling of HepG2 cells treated with ferroptosis inducer erastin and inhibitor ferrostatin (GSE104462) which was performed using the GPL16791 Illumina HiSeq 2500 (Homo sapiens) from GEO database.22 Spearman correlation analysis was employed to select the genes positively and negatively correlated with the ferroptosis level with |r| ≥ 0.9 and p < 0.05, which were designated FPI gene signatures A and B, respectively. Principal component 1 was chosen as the signature scores for FPI gene signatures through PCA. The FPI of each sample was calculated by the following formula: FPI=∑PC1A-∑PC1B. To verify its effectiveness and robustness, FPI of different subgroups in GSE104462 was compared. Also, an external cohort, GSE112384, was downloaded and analyzed as another validation cohort. GSE112384 was performed on GPL10558 platform Illumina HumanHT-12 v4.0 expression beadchip, consisted of gene expression profiles from IMR32 and SKNSH cell lines treated with withaferin A (WA, a natural ferroptosis inducing agent in neuroblastoma) and DMSO (control).23

Dimension Reduction and Generation of FSS

After the decomposition of GC samples, the differentially expressed genes (DEGs) between ferroptosis subtypes were identified in the meta-GEO cohort with the significance cutoff criteria set as adjusted p < 0.05 and |logFC| ≥ 0.585, via the “limma” package of R. The functions of DEGs were annotated based on the Reactome Knowledgebase () using the FunRich software (version 3.1.3). DEGs with positive and negative correlation to clusters were designated FS gene signatures A and B, respectively. The Boruta algorithm was conducted to further reduce noise and redundant genes.24 To quantify the ferroptosis pattern, we applied a method similar to Gene expression grade index25 to define the FSS score of each sample: FSS=∑PC1A-∑PC1B.

Difference of TME Between Ferroptosis Subtypes

To explore the difference of TME between ferroptosis subtypes, we used the R “estimate” package26 to calculate the stromal score, immune score and ESTIMATE score of each included sample and compare them grouped by subtypes. The correlations between immune-stromal score and FSS or FPI were also evaluated using correlation test. To assess the relative immune cell infiltration levels, the single sample gene set enrichment analysis (ssGSEA) programme27 was performed according to the geneset provided in the study of Charoentong,28(). The level of immune cell infiltration and gene expression of 9 immune checkpoints were compared between subtypes.

Collection of Somatic Alteration Data

With the goal of determining the tumor mutation burden (TMB) of GC patients in different ferroptosis subtypes, their non-synonymous mutations were computed according to the corresponding mutation data. The correlations between TMB and FSS or FPI were also analyzed using correlation test. Then, we compared the somatic alterations of GC driver genes which were identified within R “maftool” package29 between ferroptosis subtypes. The driver genes of top 20 mutation frequency were elected to represent the overall mutation level.

Collection of Immunotherapy Gene Expression Data

Two independent immunotherapeutic cohorts, IMvigor210 and GSE78220, were used to estimate the predictive value of the FSS and FPI. The gene expression data and clinical information of IMvigor210 cohort, which contains urothelial cancer cases received intervention of atezolizumab, an anti-PD-L1 antibody, was downloaded from under the Creative Commons 3.0 license.30 As for GSE78220, it contains metastatic melanoma cases treated with pembrolizumab, an anti-PD-1 antibody, and the expression data as well as clinical information were available from GEO ().31 Totally, 298 urothelial cancer cases and 26 metastatic melanoma cases received immunotherapy with complete clinical features were used to investigate whether the FSS and FPI could predict patients’ response to immunotherapy via comparison and prognostic analysis.

Statistical Analysis

All statistical analyses were conducted using GraphPad Prism version 8.0 or R programming. Comparison of continuous variables between two subgroups was performed using Student’s t test (normally distributed) or Mann–Whitney U-test (non-normally distributed). The prognostic analysis was conducted using the Kaplan–Meier method with the significance of differences identified utilizing Log rank tests. For FPI and FSS, the “survminer” software package in R was used to determine the optimal cut-off value. The correlation between ferroptosis subtypes and the clinical characteristics was estimated by Chi-square analysis. Correlations between two variables were examined via Spearman or Pearson correlation analysis in appropriate situations. A two-tailed p < 0.05 was considered statistically significant.

Results

Generation of FPI

With the aim of understanding the role of ferroptosis during gastric carcinogenesis, we built a FPI model to estimate the ferroptosis level based on the gene expression of a GEO dataset (GSE104462) of liver cancer cell lines treated with ferroptosis inducer erastin, or inhibitor ferrostatin. PCA was performed to define the FPI of each sample as positive principal component 1 minus negative principal component 1. As presented in , in comparison with controls, FPI was markedly elevated in the erastin group (p = 0.0023) while significantly decreased in the ferrostatin group (p = 0.0491). To validate this result, we calculated the FPI of the neuroblastoma cells, including IMR32 and SKNSH cell lines, from another GEO dataset (GSE112384), which showed a similar significant increase of FPI in the ferroptosis inducer (WA) than control group (DMSO) (). Collectively, these outcomes implied that the FPI may be a representative of the ferroptosis level based upon transcriptome data.

Identification of Two Ferroptosis Subtypes in GC by NMF

Based on the expression profiles of 121 ferroptosis-related genes in 617 GC samples from the meta-GEO cohort and 350 STAD samples from the TCGA cohort, we applied NMF consensus clustering to describe the patterns of ferroptosis in GC patients. The numbers of sample clusters decomposed by NMF depend on a desired rank k. For both primary and validation cohorts, the optimal k value turned out to be 2, that is, assigning them into two ferroptosis subtypes (meta-GEO: C1, n = 292 and C2, n = 325; TCGA: C1, n = 228 and C2, n = 122). The consensus matrix heatmap with sharp and clear edges indicated the stableness and robustness for our clustering (Figure 2A). Besides, the scatter diagram plotted by PCA also confirmed the two-dimensional distribution patterns (Figure 2B). According to the survival analysis (Figure 2C), the meta-GEO GC patients of C1 showed longer survival time than that of C2 (log rank test, p < 0.0001). Similar difference was also observed in TCGA cohort, with C1 showing better overall survival (OS) than C2 (log rank test, p = 0.0133).

Figure 2

Distribution of Ferroptosis Subtypes using NMF Consensus Clustering. (A) NMF clustering based on 121 ferroptosis-related genes decomposes the samples in meta-GEO and TCGA cohorts. The corresponding cophenetic correlation coefficient of k value between 2 and 6 is shown. (B) Scatter diagrams plotted by PCA sustained decomposition of two clusters. (C) Kaplan–Meier curves show OS for patients in C1 and C2. (D) The violin plots show FPI between C1 and C2. (E) Kaplan–Meier curves show OS for patients stratified by high FPI and low FPI.

The FPI of GC patients were also computed and drew in violin plots (Figure 2D). Compared to C1, the FPI of C2 was elevated in both primary and validation cohorts (meta-GEO: Mann–Whitney test, p < 0.0001; TCGA: Mann–Whitney test, p < 0.0001). Using the “survival” package in R, we grouped the GC patients into high FPI and low FPI subclasses according to the optimal cut-off value (Figure 2E). The Kaplan–Meier curves revealed that patients in low FPI showed better OS than patients in high FPI (meta-GEO: Log rank test, p < 0.0001; TCGA: Log rank test, p = 0.0089).

Generation of FSS

To quantify the ferroptosis pattern, FSS was constructed. Firstly, 1297 DEGs were identified between two subtypes in the meta-GEO primary cohort (). Eight hundred and forty-three genes were positively correlated with the ferroptosis subtypes (FS signature genes A) while the left 454 genes were negatively correlated (FS signature genes B). After utilizing Boruta algorithm to further reduce noise and redundant genes, 91 of FS signature genes A and 46 of FS signature genes B were reserved. Figure 3A delineates the transcriptomic profiling of the reserving signature genes in the meta-GEO cohort. Applying the screened signature genes obtained from the meta-GEO cohort, GC samples in TCGA cohort could be divided into C1 and C2 as well (Figure 3B). Functional enrichment analysis revealed that FS signature genes A were associated with epithelial–mesenchymal transition while FS signature genes B showed opposite enriched results (mesenchymal–epithelial transition) (Figure 3C). Then, FSS was constructed based on the expression profiles of FS signature genes A and B obtained from the meta-GEO cohort. FSS of each patient in TCGA cohort was also calculated using signatures the same with the meta-GEO cohort. As expected, a higher FSS was observed in C2 both in meta-GEO and TCGA cohorts (Mann–Whitney test, p < 0.0001, Figure 3D). Based on the optimal cut-off value, FSS had an ability to distinguish the prognosis of patients with GC (meta-GEO: log rank test, p < 0.0001; TCGA: log rank test, p = 0.0011, Figure 3E). FSS also exhibited a positive correlation with FPI (meta-GEO: Spearman coefficient: R = 0.67, p < 2.2e-16; TCGA: Spearman coefficient: R = 0.20, p = 0.00017, Figure 3F).

Figure 3

Generation of FSS. (A) Heatmap of DEGs between C1 and C2 to assign patients in meta-GEO cohort into two subgroups: FS gene signature A and B. (B) Heatmap shows FS gene signatures obtained from the meta-GEO cohort can classify patients in TCGA cohort into C1 and C2. (C) Reactome enrichment analysis of the FS signature genes A and B. (D) The violin plots show FSS between C1 and C2. (E) Kaplan–Meier curves show OS for patients stratified by high FSS and low FSS. (F) Scatterplots exhibit the positive correlation between FPI and FSS.

According to the ESTIMATE algorithm, compared to C1, C2 had a significantly higher level of immune score, stromal score and ESTIMATE score in the meta-GEO cohort (Figure 4A) and TCGA cohort (Figure 4B). We also unveiled that both FSS and FPI were positively correlated with immune score, stromal score and estimate score with all p-value < 0.01 (Figures 4C and D). Combining with the ssGSEA outcomes in meta-GEO (Figure 4E) and TCGA (Figure 4F) cohorts, C2 had a higher abundance in 21 immune cell populations (activated B cell, CD8 T cell, central memory CD4 T cell, effector memory CD4 T cell and CD8 T cell, gamma delta T cell, immature B cell, regulatory T cell, T follicular helper cell, type 1 and 17 helper cell, CD56dim natural killer, eosinophil, immature dendritic cell, macrophage, mast cell, MDSC, monocyte, natural killer cell, natural killer T cell and plasmacytoid dendritic cell). For immune checkpoint genes, LAG3, CTLA4, HAVCR2, CD80, IDO1, TIGIT, PDL1 and CD86 were remarkedly overexpressed in C2 while only PDCD1 showed inconsistent expression trend in these two GC cohorts (Figures 4G and H).

Figure 4

The Correlation between Ferroptosis Subtypes and TME. (A) Boxplot shows the differences of immune score, stromal score and estimate score between C1 and C2 in meta-GEO cohort. (B) Boxplot shows the differences of immune score, stromal score and estimate score between C1 and C2 in TCGA cohort. (C) Correlation of FSS and FPI with the Immune score, stromal score and estimate score in meta-GEO cohort. (D) Correlation of FSS and FPI with the Immune score, stromal score and estimate score in TCGA cohort. (E) Violin plot shows the abundance of 21 immune cell populations distinguished by C1 and C2 in meta-GEO cohort. (F) Violin plot shows the abundance of 21 immune cell populations distinguished by C1 and C2 in TCGA cohort. The green box represents C1, red one represents C2. (G) Boxplot shows the expression levels of 9 immune checkpoint genes between C1 and C2 in meta-GEO cohort. (H) Boxplot shows the expression levels of 9 immune checkpoint genes between C1 and C2 in TCGA cohort. *p < 0.05; **p  < 0.01; ***p  < 0.001.

The Correlation Between FSS, FPI and Somatic Variants

Over past decades, efforts to develop more effective treatment approaches facilitate the wide application of immunotherapy (immune checkpoint inhibition) in many malignant tumors, including GC. Accumulating evidence have illustrated that high tumor mutational burden (TMB), also known as non-synonymous variants, was closely related to tumor CD8+ T-cell infiltration and immune response such as activity of anti-PD-1 therapies.32 Although the TMB was considered with significance in predicting clinical response to immune cell infiltration, the intrinsic relationship between TMB and ferroptosis needs to be elucidated. First, the TMB of GC samples in C1 and C2 was compared. As demonstrated in Figure 5A, a higher TMB was observed in the TCGA-GC patients of C1 than that of C2 (Mann–Whitney test, p < 0.0001). Figures 5B and C exhibited that TMB were higher in low FSS or FPI than high FSS or FPI (p = 0.0006, p < 0.0001, respectively, Mann–Whitney test). After categorizing the patients into low TMB (n = 300) and high TMB (n = 43) subclasses with the optimal cut-off value, we evaluated their survival time using Kaplan–Meier curves, which indicated that prognosis of patients with low TMB was better than that with high TMB (Log rank test, p = 0.0021; Figure 5D). Then, we analyzed the correlation between TMB and FSS or FPI. As shown in Figures 5E and F, TMB exhibited a significant negative correlation with both FSS (Spearman coefficient: R = −0.29, p = 4.2e-08) and FPI (Spearman coefficient: R = −0.45, p < 2.2e-16). Moreover, the distribution of somatic alterations in GC driven genes were estimated. We listed the top 20 GC driver genes with highest mutation frequency, which showed that the overall mutation rate of driver genes in C1 was higher than C2 (Figure 5G).

Figure 5

The Correlation between Somatic Variants and FSS, FPI. (A) TMB difference between C1 and C2 in the TCGA cohort. (B) TMB difference between high FSS and low FSS in the TCGA cohort. (C) TMB difference between high FPI and low FPI in the TCGA cohort. (D) Kaplan–Meier curves show OS for patients stratified by high TMB and low TMB. (E) Scatterplots depict the negative correlations between TMB and FSS. (F) Scatterplots depict the negative correlations between TMB and FPI. (G) The oncoPrint grouped by C1 on the left (green) and C2 on the right (red), with individual patients represented as columns.

The Correlation Between Clinical Characteristics and Ferroptosis Subtypes

The relationship between clinical features and subtypes in the meta-GEO and TCGA cohorts was analyzed by Chi-square test. Table 2 illustrates that there was no significant difference between subtypes in gender both in the meta-GEO and TCGA cohorts. In terms of age and tumor grade, C2 tended to be diagnosed at a younger age and have a higher tumor grade in the meta-GEO cohort (age: p = 0.0077; tumor grade: p = 0.0003). As for Lauren type and stage, C2 was inclined to have a higher proportion of diffuse-type and associate with more advanced-stage in the TCGA cohort (Lauren type: p = 9.045e-7; stage: p = 0.0052).

Table 2

Baseline Characteristic of the Patient in Different Ferroptosis Subtypes

Characteristic	TCGA			Meta-GEO
	C1	C2	P value	C1	C2	P value
Gender			0.3757			0.6702
Female	77	47		95	111
Male	151	75		197	214
Age			0.0077			0.2257
<60	60	49		87	111
≥60	168	73		205	212
Lauren type			0.4436			9.045e-7
Intestinal	47	26		170	125
Diffuse	33	26		92	169
Mixed	1	0		28	31
Tumor grade			0.0003			NA
G1+G2	103	31		NA	NA
G3	120	87		NA	NA
Stage			0.3958			0.0052
I+II	108	52		121	100
III+IV	120	70		170	225

Ferroptosis Patterns and Ferroptosis Potential Level in the Role of Anti-PD-1/L1 Immunotherapy

The characteristics of FSS and FPI in speculating the therapeutic benefit of anti-PD-1 and anti-PD-L1 were further investigated. For FSS, both in anti-PD-L1 cohort (IMvigor210) and anti-PD-1 cohort (GSE78220), patients with low level exhibited significantly clinical benefits and a markedly prolonged survival (IMvigor210: Log rank test, p = 0.0336, Figure 6A; GSE78220: Log rank test, p = 0.0293, Figure 6D). Better therapeutic benefit in patients with low FSS compared to those with high FSS were confirmed (Figure 6B, C, E and F). Similar outcomes were observed in the FPI (Figure 6G–L). Collectively, FSS and FPI could be correlated with response to anti-PD-1/L1 immunotherapy.

Figure 6

Ferroptosis Patterns and Ferroptosis Potential Level in the Role of Anti-PD-1/L1 Immunotherapy. (A) Kaplan–Meier curves of low and high FSS patient groups in the anti-PD-L1 immunotherapy cohort (IMvigor210). (B) The proportion of patients with response to PD-L1 blockade immunotherapy in low or high FSS groups; SD: stable disease, PD: progressive disease; CR, complete response; PR, partial response. (C) Distribution of FSS in different anti-PD-L1 clinical response groups. (D) Kaplan–Meier curves of low and high FSS patient groups in the anti-PD-1 immunotherapy cohort (GSE78220). (E) The proportion of patients with response to anti-PD-1 immunotherapy in low or high FSS groups. (F) Differences in FSS among different anti-PD-1 clinical response groups. (G) Kaplan–Meier curves of low and high FPI patient groups in the anti-PD-L1 immunotherapy cohort (IMvigor210). (H) The proportion of patients with response to anti-PD-L1 immunotherapy in low or high FPI groups. (I) Distribution of FSS in different anti-PD-L1 clinical response groups. (J) Kaplan–Meier curves of low and high FPI patient groups in the anti-PD-1 immunotherapy cohort (GSE78220). (K) The proportion of patients with response to PD-1 blockade immunotherapy in low or high FPI groups. (L) Differences in FPI among different anti-PD-1 clinical response groups. The statistical significance of differences in all survival analysis was determined via the Log rank test.

Discussion

Although immunotherapy has opened a new frontier in GC treatment, only a part of patients could obtain clinical benefits, which demanded novel biomarkers enabling to predict the response to immunotherapy. Ferroptosis has been reported closely related to tumorigenesis and immunity but with a dual role. Therefore, identifying distinct ferroptosis patterns is helpful to construct an effective immunotherapy strategy.

In the current study, two ferroptosis subtypes (C1 and C2) with different prognosis and ferroptosis level in GC were identified with a relationship to the characteristics of TME. Compared with C1, C2 showed higher immune score and stromal score. Besides, C2 had a higher abundance of immune cell infiltration and was significantly enriched in epithelial–mesenchymal transition, which suggested a correlation with active immune and stromal activities. It was reported that immune infiltration in the TME is generally characterized into three patterns:33 1) Immune-inflamed: dense infiltration of immune cell within the tumor ensures the high probability of response to immunotherapy, 2) Immune-excluded: a large number of immune cells infiltrated around the tumor but not penetrating inward lowers the probability of response, 3) Immune-desert: lacking of T cells in both tumor cells and peripheral matrix leads to immune tolerance or immune ignorance. According to the exhibitions that C2 with a higher abundance of immune cell infiltration but lower probability in response to immunotherapy and associated with stromal activities, it was reasonable to speculate that C2 was in accordance with the immune-excluded phenotype. As a successful immunotherapy depends on the ability of innate and adaptive immune cells to penetrate into the tumor parenchyma and eradicate cancer cells,34 in the background of immune-excluded phenotype, activated immune cells cannot effectively penetrate into the tumor parenchyma to interact with tumor cells, resulting in no intended targets for immunotherapy. Conclusively, although C2 was described with a more active immunity and a higher expression of immune checkpoints, it was not surprising C2 tended to carry a poorer prognosis. Previously, we developed a ferroptosis-related model to predict the prognosis of GC and found that the high-risk group also exhibited higher Immune Score, Stromal Score, and ESTIMATE Score, suggesting higher levels in immune cells infiltration and stromal component in the TME of high-risk group.35 Additionally, the high-risk group was also speculated in accordance with the immune-excluded subtype. Therefore, we inferred that C2 patients in this study have a high-risk score in the ferroptosis-related model. As for C1, it was reasonable to deduce as immune-inflamed phenotype because of lower stromal activities but increased neoantigen load, enhanced response to immunotherapy and relatively better prognosis.

In the meta-GEO cohort, we observed that C2 was significantly associated with diffuse-type while C1 with intestinal-type. The TME characteristic of diffuse-type GC is similar to immune-excluded phenotype as it is accompanied with extensive stromal fibrosis and abundant extracellular matrix,36,37 which also impede the accessibility of T-cells.38 As for intestinal-type, it is considered as the result of Correa pathway under chronic inflammation stress,39 which increases the infiltration of immunosuppressive cells and results T cell exhaustion.40

There also existed a negative correlation between FSS and TMB in GC. Previous studies have indicated that gene mutations were associated with cancer immunotherapy response.41–43 It was also reported that tumors with higher levels of TMB are more likely to produce new immunogenic neo-antigens, making tumors more susceptible to immunotherapy.44 Our analysis demonstrated that the FSS integrated not only the characteristics of stromal and immune TME but also tumor mutation load, which could serve as an accurate predictor for patient response to immunotherapy and thereby direct the individual immunotherapy of GC.

To evaluate the ferroptosis level, we constructed a FPI model. Similar to FSS, a higher FPI was observed in C2 and carried with a poorer prognosis. Besides, the correlations of FPI with the characteristics TME and TMB resemble to FSS. A higher FPI was associated with more active immune and stromal activities but a lower probability of immunotherapy response. Liu et al21 systematically analyzed the aberrances and functional implications of ferroptosis in cancer and found that epithelial–mesenchymal transition (stromal-related pathway) and IL-6/JAK/STAT3 (immune-related pathway) signaling were enriched in the high-FPI group. Besides, they also demonstrated that higher FPI was correlated with a worse survival in glioblastoma multiforme (GBM), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), and lung adenocarcinoma (LUAD). The unconventional phenomena that a high ferroptosis level was associated with worse status in GC, GBM, KIRC, KIRP, LIHC and LUAD indicated that ferroptosis plays a dual role in tumor immune response. Also, regulating the ferroptosis level in TME with different strategies may benefit patients and improve prognosis.

In our study, C2 corresponding to diffuse-type GC had a higher level of FPI compared with C1 corresponding to intestinal-type GC. It was reported that low-grade inflammation could increase the cell population size due to a decreased rate of DNA damage-independent apoptosis, whereas higher grade inflammation could suppress cell growth through an increasing rate of DNA damage-induced apoptosis.45 These may be due to different levels of ROS corresponding to different grades of inflammation. At low levels, ROS could exert a protective effect through activating protective signaling pathways against inflammation and increase cell proliferation as well as survival but at high levels ROS cause DNA damage and trigger apoptosis.46 Ferroptosis, an iron-dependent form of regulated cell death, is also induced by the accumulation of lethal lipid ROS.6 Diffuse-type GC characterized with highly active inflammation47,48 may have a higher level of ROS, which induce ferroptosis in a higher proportion of cells, and therefore a higher FPI was observed. For intestinal-type GC, it undergoes chronic inflammatory process49,50 and is accompanied with a relatively lower level of ROS. As a result, the proportion of cells that undergo ferroptosis is smaller.

There were some limitations that must be emphasized. Firstly, our findings mainly come from retrospective data and prospective clinical validation in a larger GC cohort receiving immunotherapy is demanded henceforth. Secondly, the mechanism of ferroptosis shaping TME characteristics in GC was unknown and should be investigated experimentally in the future.

In conclusion, we analyzed the ferroptosis patterns in GC and its relationship with TME. The difference of ferroptosis patterns possessing significantly different ferroptosis levels was associated with the heterogeneity and complexity of individual TME. Therefore, the systematic assessment of ferroptosis patterns in this study may be helpful for the further understanding of TME and guiding immunotherapy.