Public gene-expression data and complete clinical information were obtained through Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) and the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. Patients without survival information were removed from further evaluation. In total, five eligible GC cohorts [GSE15459, GSE34942, GSE62254/ACRG, GSE84437, TCGA-STAD (The Cancer Genome Atlas-Stomach Adenocarcinoma)] were gathered in this study for further analysis. Gene expression data and the clinical data of 90 Tianjin GC samples were obtained from previous studies (20). In short, six cohorts, including 1,404 samples, were integrated into the study as a training set. In addition, GSE26899 and GSE26901 datasets were included in the study as two independent validation cohorts. The baseline information of all eligible GC datasets was summarized in Supplementary Table S1. Ferroptosis-related genes (FRGs) were identified from published literature (21).

Unsupervised clustering analysis was conducted to identify distinct clusters according to the expression of 60 ferroptosis-related genes. We used the R package ConsensusClusterPlus (v 1.58.0) to execute the consensus clustering algorithm and repeat it 1,000 times to ensure the clustering stability (22).

Differentially expressed genes (DEGs) (Supplementary Table S2) between different clusters were determined using the empirical Bayesian approach of the limma (v 3.50.0) R package (23). The significance criterion for determining DEGs was |log fold change (FC)| > 1 and adjusted P-value < 0.05. The Benjamini and Hochberg method was applied for multiple-testing adjustments.

First, using univariate Cox regression analysis, selected 138 differentially expressed genes with prognostic value (Supplementary Table S3) for further analysis (P < 0.01). We named these 138 genes as “subtype-specific genes.”

Next, we performed unsupervised clustering of all samples using the expression values of 138 subtype-specific genes via the R package “ConsensusClusterPlus.” Two ferroptosis subtypes were obtained, named “Subtype-A” and “Subtype-B.”

Patients in the two subtypes were stratified according to gender, and the Kaplan-Meier survival analysis was used to assess the gender differences in overall survival (OS). The log-rank test was used to evaluate the statistical significance of the Kaplan-Meier curves at P < 0.05.

To further explore the factors of gender differences between the two ferroptosis subtypes, we analyzed the expressions of three sex hormone receptors (AR, ER, PR) and HER2 of the two subtypes using the Wilcoxon rank-sum test. Furthermore, we stratified the two subtypes by gender and compared their differences in hormone receptor levels, respectively.

We recognized the significantly mutated genes (SMG) with the GenVisR (v 1.26.0) tool in the R package. Mutation signature analysis of two ferroptosis subtypes was conducted using the R package MutationalPatterns (v 3.4.0) and maftools (v 2.10.0). We extracted the mutational signature of GC data and compared them with the mutation database (COSMIC V2) by using the cosine similarity method (https://cancer.sanger.ac.uk/cosmic).

Gene set variation analysis (GSVA) (24) was performed with the R package “GSVA”(v 1.42.0) to evaluate pathway enrichment for two ferroptosis subtypes. All hallmark gene sets were downloaded from the Molecular Signature Database (MSigDB). Adjusted P with a value < 0.05 was considered statistically significant. The Benjamini and Hochberg method was applied for multiple-testing adjustments.

The single-sample gene set enrichment analysis (ssGSEA) (25) was introduced to quantify the relative infiltration of 28 immune cell types in the tumor microenvironment. Unique feature gene panels for each immune cell subset were obtained from the latest literature (26). An enrichment score in ssGSEA analysis represented the relative abundance of each immune cell type. The ssGSEA score was normalized to unity distribution, for which zero is the minimal and one is the maximal score for each immune cell type. The bio-similarity of the immune cell filtration was estimated by multidimensional scaling (MDS) and a Gaussian fitting model.

In addition, we assessed the abundances of tumor immune infiltration in subpopulations of M0 and M1 macrophages using the CIBERSORT algorithm. Then, the infiltration abundance of M2a, M2b, M2c, and M2d macrophage subpopulations was analyzed using the ssGSEA algorithm based on four subtype-related markers (27).

Principal component analysis (PCA) was conducted to construct the ferroptosis relevant risk score, which we called was ferroptosis score. Principal components 1 and 2 were selected as signature scores. The advantage of this method is that the score is concentrated on the set with the largest block of well-correlated (or anticorrelated) genes in the set, while down-weighting contributions from genes that do not track with other set members (28, 29):

ferroptosis score = ∑(PC1_i + PC2_i)

where i is the expression of subtype-specific genes.

The GC samples were categorized into high- and low-ferroptosis score groups according to the optimal cut-off value determined by the survminer R package. Then, we performed Kaplan-Meier analysis of ferroptosis score in GC cohort (n = 1,404) and calculated the Hazard Ratio (HR) by univariate and multivariate cox regression analysis to further assess the predictive accuracy of the model. Furthermore, we did external validation in two independent cohorts: GSE26899 (n = 93) and GSE226901 (n = 109). We fitted the ferroptosis score with the meaningful clinical indicators obtained by multivariate cox regression, and the accuracy of the ferroptosis score was evaluated using the receiver operating characteristic (ROC) curves. The R package “glmnet” was used for Logistic regression analysis for multivariate fitting, and the R package “timeROC” was used to draw ROC curves.

To investigated the association between ferroptosis score and biological factors, we performed correlation analysis to reveal the association between ferroptosis score and AR, ER, PR, HER2, and 28 immune cell types.

We compared three remarkable molecular biomarkers, including TMB, MSI, and TIDE, between high- and low-ferroptosis score groups. Microsatellite instability (MSI) status was classified as Microsatellite stable (MSS), MSI-low (MSI-L, one marker unstable), and MSI-high (MSI-H, over two markers unstable). The data of MSI status were obtained from the UCSC Xena database.

In order to validate the efficiency of the ferroptosis score in males and females, an immunotherapeutic cohort was included in our study: advanced urothelial cancer with the intervention of atezolizumab, an anti-PD-L1 antibody (IMvigor210 cohort) (30). Clinical information and gene expression data were extracted from the IMvigor210 data set (http://research-pub.gene.com/IMvigor210CoreBiologies).

The R package “pRRophetic” (v 0.5) (31) was used to predict chemotherapeutic response in GC patients. The half-maximal inhibitory concentration (IC50) of the samples was demonstrated using ridge regression, and the prediction accuracy was assessed using 10-fold cross-validation based on the GDSC training set (32).

Clustering of ferroptosis was performed using the consensus clustering method to identify structure over many clustering realizations robustly (33, 34). The survival curves for the prognostic analysis were generated via the Kaplan-Meier method, and log-rank tests were used to identify the significance of differences. Univariate and multivariate Cox regression models were employed for calculating hazard ratio (HR) and forest plots visualized the coefficients of those regression models. Continuous variables were compared between two groups through was performed by Wilcoxon rank-sum test. Classified variables were compared between two groups through was performed by chi-square test. Kruskal-Wallis test was used to conduct difference comparisons of three or more groups. Correlations between variables were assessed via Pearson and Spearman coefficient. Differences were considered statistically significant when a P-value was < 0.05, and all P-values reported were two-sided. All statistical analyses were conducted using R 4.0.0 (https://www.r-project.org/).

To systematically describe our study, we developed a flow chart (Figure 1). A total of 1,404 samples from six independent GC cohorts accompanied with complete survival information were retained for the present study. All GC cases were divided into two clusters by unsupervised clustering of the expression of 60 ferroptosis-related genes (FRGs). Survival analysis indicated a difference in prognosis between the two clusters (log-rank test, P = 0.039; Supplementary Figure S1).

To construct a more precise ferroptosis subtype with prognostic significance, we recognized 278 differentially expressed genes (DEGs) among clusters and certified 138 DEGS with prognostic value by univariate Cox regression analysis. Subsequently, we carried out secondary clustering that the consensus cluster analysis indicated that the optimal number of clusters was two, defined by CDF curves (Figures 2A,B). Unsupervised clustering analysis based on the expression of these 138 genes also divided GC patients into two clusters, which we called Subtype-A and Subtype-B.

Therefore, the 138 genes are called “subtype-specific genes” (Figure 2C).

Additionally, survival analysis indicated that Subtype-A had a better prognosis (log-rank test, P = 0.0012; Figure 2D). Considering that there may be gender differences in individual ferroptosis levels (19), we stratified the two subtypes by sex and conducted a survival analysis to observe the prognosis. As shown in Figure 2E, male patients with gastric cancer had the worst prognosis in Subtype-B (log-rank test, P = 0.0037).

To further dissect the factors involved in gender differences, we examined three sex hormone receptors (AR, ER, and PR) and HER2 between Subtype-A and Subtype-B. The results presented in Figure 3A showed that all hormone receptors differ significantly between the two subtypes (Wilcoxon rank-sum test, all P < 2.2e-16).

Furthermore, we compared three sex hormone receptors (AR, ER, and PR) and HER2 between the two subtypes in males and females, respectively. The results were consistent with the results of the above overall analysis, and there were significant differences in the four hormone receptors among the two ferroptosis subtypes in both males and females (Figures 3B,C).

We performed SMG analysis for GC samples in Subtype-A and Subtype-B. The SMG mutational landscapes of these two subgroups (Figures 4A,B) exhibited a distinct mutation ratio, and the top three were the same between subtypes: TTN (59.4 vs. 36.2%), TP53 (49.0 vs. 35.2%), and MUC16 (37.7 vs. 28.6%). The fourth place in Subtype-A was LRP18 (31.0%), while in Subtype-B was SYNE1 (22.9%).

To gain further insights into the operative mutational processes in the two subtypes, we extracted the mutational signatures from the COSMIC database using GC's genomic somatic mutation data. We extracted four mutational signatures (i.e., signatures 17, 10, 6, and 21; Figure 4C) from mutation data of Subtype-A. We extracted four mutational signatures (i.e., signatures 3, 1, 6, and 17; Figure 4D) from mutation data of Subtype-B. Subtype-A had the independent characteristics of signature 10 and signature 21, while Subtype-B had the independent characteristics of signature 3 and signature 1. These results suggested that the mutation pattern of Subtype-B was associated with DNA damage and repair pathways such as homologous recombination, leading to DNA double-strand break repair failure.

To explore the differences in biological behaviors between the two subtypes, we performed a GSVA enrichment analysis. As shown in Figure 5A, the two subtypes showed some differences, E2F_TARGETS as the top enriched signature in Subtype-A and UV_RESPONSE_DN as the top enriched signature in Subtype-B. The results from GSVA analyses have revealed Subtype-B was significantly associated with genes down-regulated in response to ultraviolet (UV) radiation, Subtype-A was significantly associated with genes encoding cell cycle-related targets of E2F transcription factors. Therefore, we speculated that Subtype-B has more obvious DNA damage.

We next composed a heatmap to visualize the relative abundance of immune infiltrating cell subpopulations from the discovery dataset (Figure 5B). We found a significant difference in immune cell infiltration between the two subtypes (Figure 5C). In addition, it can be observed that there are differences in immune infiltration among macrophage subpopulations of ferroptosis subtypes. There was no difference in M1 macrophage abundance between the two subtypes. M0 and M2b subpopulations were significantly enriched in the Subtype-A (both P < 0.001), while M2a (P < 0.01), M2c (P < 0.001), and M2d (P < 0.001) subpopulations were significantly enriched in the Subtype-B.

To better predict the prognosis and guide individualized treatment, we construct a ferroptosis-related risk model. The ferroptosis score was used to quantify the ferroptosis levels of individual patients with gastric cancer, considering the individual heterogeneity and complexity of gene expression levels. The ferroptosis score was modeled based on the expression of subtype-specific genes by PCA algorithms. Supplementary Figure S2 visualizes the contribution of two components of the PCA algorithm.

Next, we sought to identify the value of the ferroptosis score in predicting the outcome of GC patients. Patients were divided into high- and low-ferroptosis score groups with the optimal cut-off value determined by the survminer R package. We performed Kaplan-Meier survival analysis and univariate cox regression analysis, and the result showed that patients with low-ferroptosis scores demonstrated a prominent survival benefit [log-rank test, P < 0.0001; HR 1.54 (1.32-1.79); Figure 6A]. So “low-ferroptosis score” means “risk-low.” We then tested whether the ferroptosis score could serve as an independent prognostic biomarker for gastric cancer. Multivariate Cox regression analysis, which included the factors of patients' age, gender, tumor stage, confirmed ferroptosis score is an independent prognostic indicator of evaluating patient outcomes [HR 1.33 (1.08-1.65); Figure 6B]. To further test the stability of the risk model, the predictive value was externally validated in the GSE26899 and GSE26901 datasets. Consistent with the training set, the survival analysis results and univariate and multivariate cox regression analysis of GSE26899 (Figures 6C,D) and GSE26901 (Figures 6E,F) datasets were all statistically significant. Hence, the ferroptosis score could be an independent prognostic factor for gastric cancer. ROC curve showed that the ferroptosis score in our present study exhibited a good or acceptable predictive value in the aspect of 1 year AUC, 3 years AUC, and 5 years AUC, which was 0.68, 0.67, and 0.67, respectively (Supplementary Figure S3A). Moreover, in the two verification sets, the ROC curve also validated a good predictive value in the survival time of 1 year, 3 years, or 5 years (Supplementary Figures S3B,C).

We examined the specific correlation between ferroptosis score and four hormone receptors using Pearson correlation analyses. The correlation analysis revealed that the ferroptosis score was significantly associated with three sex hormone receptors (AR, ER, PR) and HER2. In detail, ferroptosis score was positively associated with AR (r = 0.611, P = 2.03e-144; Figure 7A), ER (r = 0.281, P = 7.64e-27; Figure 7B), PR (r = 0.747, P = 1.02e-250; Figure 7C), and negatively associated with HER2 (r = −0.247, P = 5.46e-21; Figure 7D).

Then, we performed Spearman correlation analysis to explore the relationship between ferroptosis score and immune cell infiltration. As shown in Figure 7E, there was a widespread correlation between the ferroptosis score and 27 out of 28 kinds of immune cell types. And we further analyzed the correlation between ferroptosis score and macrophage subpopulations. Ferroptosis score was positively correlated with M2c and M2d macrophages and negatively correlated with M2b and M0 macrophages. Moreover, we found that the strongest positive correlation was the plasmacytoid dendritic cell (r = 0.51, P < 0.0001; Figure 7F), and the strongest negative correlation was the activated CD4 T cell (r = −0.48, P < 0.0001; Figure 7G).

We analyzed intrinsic immune escape mechanisms between high-and low-ferroptosis score groups in GC patients. The results showed that TMB was relatively higher in the low-ferroptosis score group (Wilcoxon rank-sum test, P = 1.239e-07; Figure 7H), indicating the relatively high immunogenicity. Meanwhile, the low-ferroptosis score group tended to have more MSI-H status (chi-square test, P = 0.00207; Figure 7I). The above results suggested that patients with low ferroptosis scores were more suitable for immunotherapy. Immunotherapy has manifested improved survival in treating multiple tumor types, and it is urgent to identify patients who will benefit most. Thus, we used the TIDE algorithm to further explored immunotherapy response prediction (35). The results showed that the low-ferroptosis score group had a better response to immunotherapy (Figure 8A). The result was consistent with those of TMB and MSI.

Afterward, we included an immunotherapy cohort to investigate whether ferroptosis scores predicting the patients' responsiveness to immunotherapy differ in males and females. The results showed that males are more sensitive to ferroptosis scores and more likely to benefit from immunotherapy (Figures 8B,C).

Additionally, we performed a prediction analysis of response to chemotherapy in the two groups by applying the pRRophetic R package. Patients in the low-ferroptosis score group had a higher sensitivity for the following chemotherapy drugs: paclitaxel (Wilcoxon rank-sum test, P < 2.2e-16; Figure 8D), gemcitabine (Wilcoxon rank-sum test, P = 1.6e-09; Figure 8E), methotrexate (Wilcoxon rank-sum test, P < 2.2e-16; Figure 8F). Patients in the high-ferroptosis score group had a higher sensitivity for bicalutamide (Wilcoxon rank-sum test, P = 0.00016; Figure 8G).

The Asian Cancer Research Group (ACRG) divides GC into four subtypes using gene expression data: microsatellite instability (MSI), epithelial-to-mesenchymal transition (EMT), MSS/TP53+ and MSS/TP53– (36). The subsets described by ACRG were associated with distinct genomic alterations, survival outcomes, and postoperative recurrence patterns (36). Nevertheless, so far, the possibility of GC molecular classification in daily practice and therapeutic implementation is still challenging (37). Our study found no ferroptosis score specificity among ACRG's GC molecular subtypes (Supplementary Figure S4), indicating that our score system differed from ACRG's subtype classification system. It may also mean that better treatment and prognosis can be suggested, but more exploration is needed in the future. Overall, our study strongly confirmed that the ferroptosis score could be used to assess patients' prognosis and immunotherapy response and may be different for males and females. Our work may be helpful to deepen the understanding of ferroptosis and provide a reference for clinical management and targeted treatment of GC.