A Novel Ferroptosis-Related LncRNA Pair Prognostic Signature Predicts Immune Landscapes and Treatment Responses for Gastric Cancer Patients

Background: The construction of ferroptosis-related lncRNA prognostic models in malignancies has been an intense area of research recently. However, most of the studies focused on the exact expression of lncRNAs and had limited application values. Herein, we aim to establish a novel prognostic model for gastric cancer (GC) patients and discuss its correlation with immune landscapes and treatment responses.

Methods: The present study retrieved transcriptional data of GC patients from the Cancer Genome Atlas (TCGA) database. We identified differentially expressed ferroptosis-related lncRNAs between tumor and normal controls of GC samples. Based on a new method of cyclically single pairing, we constructed a 0 or 1 matrix of ferroptosis-related lncRNA pairs (FRLPs). A risk score signature consisting of 10 FRLPs was established using multi-step Cox regression analysis. Next, we performed a series of systematic analyses to investigate the association of the FRLP model and tumor microenvironment, biological function, and treatment responses. An alternative model to the FRLP risk score signature, the gene set score (GS) model was also constructed, which could represent the former when lncRNA expression was not available.

Results: We established a novel prognostic signature of 10 ferroptosis-related lncRNA pairs. High-risk patients in our risk score model were characterized by high infiltration of immune cells, upregulated carcinogenic and stromal activities, and heightened sensitivity to a wide range of anti-tumor drugs, whereas low-risk patients were associated with better responses to methotrexate treatment and elevated immunotherapeutic sensitivity. The practicability of the FRLP risk score model was also validated in two independent microarray datasets downloaded from Gene Expression Omnibus (GEO) using the GS model. Finally, two online dynamic nomograms were built to enhance the clinical utility of the study.

Conclusion: In this study, we developed a ferroptosis-related lncRNA pair-based risk score model that did not rely on the exact lncRNA expression level. This novel model might provide insights for the accurate prediction and comprehensive management for GC patients.

Introduction

Gastric cancer (GC) is one of the most common malignancies worldwide. In consideration of its frequently advanced stage at diagnosis, GC is also the fourth leading cause of cancer-related deaths, with 768,793 deaths globally in 2020 (Sung et al., 2021). In China, GC remains the second most deadly cancer. According to statistics, there were 679,100 new GC cases and approximately 498,000 deaths occurred in 2015 (Chen et al., 2016). Currently, surgical resection is still the first-line treatment for GC patients, and the administration of adjuvant or neoadjuvant chemotherapy and immunotherapy could further improve patients’ prognosis. Unfortunately, resistance to chemotherapy occurs frequently, which remains a major cause of GC metastasis or relapse. Besides, the effect of immunotherapy is quite uncertain for GC patients due to the highly variant tumor microenvironment (TME) (Coutzac et al., 2019). Recent studies have linked ferroptosis with immune factors and chemotherapy resistance. Thus, it is necessary to further explore such an association in detail in GC.

Ferroptosis was firstly proposed by Dixon et al., (2012) and was described as an iron-dependent non-apoptotic cell death triggered by the small molecule erastin. The team summarized in another essay that ferroptosis was associated with reactive oxygen species production and lipid peroxidation (Dixon, 2017). In the past few years, ferroptosis has been implicated in multiple diseases and functions as a tumor suppression mechanism (Stockwell et al., 2020). Exploring the mechanisms underlying susceptibility and resistance to ferroptosis in carcinogenesis has been an intense area of research in the past few years (Friedmann Angeli et al., 2019). It was revealed by some studies that some tumor suppressor genes exert their anti-tumor function by upregulating tumor cells’ sensitivity to ferroptosis. For example, p53 (Wang et al., 2016) and BAP1 (Zhang et al., 2018) were found to downregulate the expression of SLC7A11, a negative modulator of ferroptosis. It was thus hypothesized that these two tumor-suppressive genes executed their anti-tumor function partly by enhancing cancer cells’ sensitivity to ferroptosis. Conversely, carcinoma’s resistance to ferroptosis was also reported to be connected with some of the malignant signatures in tumorigenesis such as hypoxia (Luis et al., 2021), EMT (Behan et al., 2019), and stemness (Wang et al., 2020). Studies have already highlighted the possibility for ferroptosis to be a novel target for cancer treatment (Liang et al., 2019; Koppula et al., 2021). When it comes to gastric cancer, the induction of ferroptosis has also been achieved in several studies (Hao et al., 2017; Chen et al., 2020; Cai et al., 2021; Zhao et al., 2021). These findings may provide new insights into ferroptosis-mediated cancer treatment.

Long non-coding RNAs (lncRNAs) are defined as non-protein coding transcripts of over 200 nucleotides. It was revealed that lncRNAs were aberrantly expressed in tumor tissues and were involved in the carcinogenesis of various cancer types (Bhan et al., 2017). Up to now, there have been some studies establishing lncRNA-based, ferroptosis-related risk signatures in relation to GC prognosis (Pan et al., 2021; Wei et al., 2021; Xiao et al., 2021; Zhang et al., 2022). Table 1 lists recent works that constructed ferroptosis-related lncRNA risk score models in GC. It is worth noting that some recent studies cyclically singly paired the ferroptosis-related lncRNAs and constructed prognostic risk score models of lncRNA pairs (Li et al., 2021; Tang et al., 2021). Compared with the aforementioned risk score model consisting of single lncRNA, these models did not rely on the specific expression of lncRNA and had broader clinical practicability, for it is unnecessary to transform data format when the expression profiles were obtained from different sequencing platforms.

TABLE 1

TABLE 1. Recent Works Constructing Ferroptosis-related lncRNA Risk Score Models in GC patients.

In this study, we successfully established a novel ferroptosis-related lncRNA pair (FRLP) risk score model of clinical significance. An alternative gene set score (GS) model was also built to represent the FRLP model when the lncRNA expression profiles were not available. After confirming that the two models shared a high degree of compliance, we validated the FRLP risk score model in two external GEO cohorts using the GS model. Our present work could provide not only an accurate prediction for GC patients’ survival but also insightful strategies to optimize GC patient’ treatment.

Materials and Methods

Data Obtaining and Processing

The transcriptional data of 375 stomach adenocarcinoma (STAD) tissues and 32 normal tissues were downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov) database. We also obtained the corresponding clinical data from the University of California Santa Cruz (UCSC; https://xena.ucsc.edu) database. Patients without clinical information were excluded. Overall, a total of 334 cases were included. The clinical information of patients, including age, gender, TNM status, stage, and tumor location, is listed in Table 2. Next, we transformed the Ensembl ID to gene symbols and picked out the lncRNAs using GTF annotation document downloaded from the ensemble human genome browser GRCh38. p13 (http://asia.ensembl.org/index.html). Ferroptosis-related genes were retrieved from the FerrDb (http://www.zhounan.org/) database.

TABLE 2

TABLE 2. Characteristics of patients with STAD from TCGA database, GSE62254 and GSE84437.

Screening of Ferroptosis-Related lncRNAs

The frlncRNAs were identified by the Pearson correlation analysis (|R|>0.4 and p < 0.001). Next, we filtered out differentially expressed lncRNAs between tumor and normal controls using the “limma” package of R software under the criteria of |log fold change (FC)|> 2 and adjusted p value < 0.05. The volcano plot and heatmap of DEfrlncRNAs were depicted by the “ggplot2” and “pheatmap” packages of R software.

LncRNA Pairing

All the DEfrlncRNAs were cyclically single-paired and a 0 or 1 matrix was established according to the following procedure (Luis et al., 2021): if the expression of lncRNA a was higher than its pairing lncRNA, lncRNA b, the value of this lncRNA pair was defined as 1; conversely, we regarded the pair’s value as 0 if lncRNA A’s expression was lower than lncRNA B. FrlncRNA pairs (FRLPs) were excluded if there were too many (more than 80%) 0 or 1 values. In other words, only the lncRNA pairs with an optimal 0 or 1 range (20%–80%) were sorted out for further analysis. Univariate Cox regression analyses were performed to evaluate the prognostic value of FRLPs (p < 0.01).

Risk Score Model Construction

Using R package “glmnet”, the least absolute shrinkage and selection operator (LASSO) regression analysis was conducted in TCGA cohort, and the lncRNA pairs filtered out were further subjected to multivariate Cox regression analysis. The multicollinearity of lncRNA pairs was estimated through the variance inflation factor (VIFs), and we defined that VIF ≥2 was considered to indicate multicollinearity in the study (Yan et al., 2021). LncRNA pairs that did not violate the multicollinearity assumption were filtered out for model construction. Risk scores were calculated based on the formula below:

$Riskscore={\displaystyle \sum _{i=1}^n\mathbf{C}\mathbf{o}\mathbf{e}\mathbf{f}\left(\ln \mathbf{c}\mathbf{R}\mathbf{N}\mathbf{A}\mathbf{p}\mathbf{a}\mathbf{i}\mathbf{r}\mathbf{i}\right)}\ast \mathbf{V}\mathbf{a}\mathbf{l}\left(\ln \mathbf{c}\mathbf{R}\mathbf{N}\mathbf{A}\mathbf{p}\mathbf{a}\mathbf{i}\mathbf{r}\mathbf{i}\right)$

Coef (lncRNA pair i) and Exp (lncRNA pair i) indicated the coefficient and value (0 or 1) of an individual lncRNA pair I, respectively. Time-dependent receiver operating characteristic (ROC) curve (Kamarudin et al., 2017) and decision curve analysis (DCA) (Vickers and Elkin, 2006) were employed to access the predictive efficacy of the FRLP risk score model. Then, patients were divided into high- and low-risk groups based on the maximum reflection point of the most suitable ROC curve. Comparisons of differences in overall survival (OS) were conducted using the Kaplan–Meier method.

Development and Evaluation of Nomogram Based on the Risk Score

To quantitatively estimate our FRLP risk score model’s prognostic potential in clinical practice, we created a nomogram using “Survival” and “RMS” packages in R to predict 1, 3-, or 5-year OS. The nomogram integrated risk score and other clinicopathological variables, which were also identified as independent risk factors by the Cox regression analysis. In addition, using the “coxph” function of R package “survival,” we performed Schoenfield’s residual test to decide whether these risk factors met the equally proportional risk hypothesis. Finally, the robustness of the nomogram was also evaluated by the calibration curves, time-dependent ROC curves, and DCA analysis (Vickers and Elkin, 2006).

Clinical Correlation Analysis

To explore the clinical significance of our risk score model, we analyzed the difference in the clinicopathological characteristics between the two risk groups by the chi-square test. The result was depicted in a heatmap form using R packages “pheatmap” and “ggpubr.” The clinicopathological characteristics include age, gender, tumor grade, tumor location, and TNM status. Moreover, we conducted a stratification analysis comparing the risk score difference in different clinical groups: age (>65 and ≤65), gender (male and female), tumor stage (I–II and III–IV), grade (G1–G2 and G3), tumor location (proximal, body/fundus and distal), T status (T1–T2 and T3–T4), N status (N0 and N1–N3), and M status (M0 and M1).

Somatic Variation Analysis

Using VarScan2 annotation files downloaded from TCGA database, the tumor mutation burden (TMB) of each sample was calculated through the VarScan2 pipeline somatic mutation calling workflow (Binatti et al., 2021). A comparison of TMB between the two groups was carried out. Survival analysis combining the risk score and TMB was performed. Spearman correlation analysis was performed to test the relation between TMB and risk scores. The top 20 genes with the highest mutation frequencies were visualized using the “maftools” package in R.

Tumor Microenvironment Analysis

The stromal and immune cell content in the microenvironment was quantified for every STAD patient using the “ESTIMATE” algorithm, and a comparison was performed between the two risk groups. Next, several currently acknowledged algorithms include TIMER, xCell, quanTIseq, MCP-counter, EPIC, CIBERSORT-ABS, and CIBERSORT. to calculate the contents of tumor infiltrating immune cells (TIICs). The whole analysis process was performed via the online platform Tumor Immune Estimation Resource 2.0 (TIMER2.0, http://timer.cistrome.org/). Spearman correlation analysis was further performed to analyze the correlations between the risk score and immune cells. We set “R > 0.1, p < 0.05” as the criterion and visualized our results in the form of a lollipop diagram using the R package “ggplot2.”

Immunotherapeutic Sensitivity Analysis

To explore the correlation between risk score and response toward the immune checkpoint blockage (ICB) treatment, we obtained 24 immune checkpoint genes (ICGs) from previous literatures (Garris et al., 2018; Han et al., 2019; Nishino et al., 2017; Patel et al., 2017; Yang et al., 2017). The tumor immune dysfunction and exclusion (TIDE) score was defined by Jiang and his colleagues (Jiang et al., 2018) and has been proved to be a reliable indicator for predicting the ICB treatment response. We obtained the TIDE score, dysfunction score, and exclusion score of each patient from the TIDE website (http://tide.dfci.harvard.edu/). In addition, the immunophenoscore (IPS), calculated by unbiased machine learning methods, is another parameter reflecting immunogenicity. Higher IPS represents better accuracy for the more corresponding result (Jiang et al., 2021). The IPS results for anti-CTLA4 and anti-PD1 treatments of TCGA STAD patients were downloaded from The Cancer Immunome Atlas (TCIA, https://tcia.at/home). Two groups’ differences in IPS and TIDE scores were compared.

Biological Function Analysis

To shed light on the difference in biological functions between the two risk groups, “KEGG” gene sets “c2. cp.kegg.v7.0. symbol.gmt” and “GO” gene sets “c5. go.v7.4. symbols.gmt” were downloaded from the MsigDB (http://www.gsea-msigdb.org) database for gene set enrichment analysis (GSEA) (Subramanian et al., 2005). Besides, “Hallmark” gene set “h.all.v7.4. symbols.gmt,” which contains 50 well-defined biological signatures, were also available from MsigDB (Subramanian et al., 2005). We estimated the enrichment score (ES) of each signature using single sample gene set enrichment analysis (ssGSEA).

Chemotherapy Response Analysis

Half-inhibitory concentration (IC50) of different types of chemotherapeutic drugs was estimated using R package “pRRophetic.” The “pRRophetic” package selected 138 kinds of drugs from more than 700 cell lines in the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/) and developed a ridge regression algorithm to predict treatment responses. Here, we only chose gastrointestinal cell lines to ensure a more accurate prediction by setting the parameter “tissueType” of the “pRRopheticPredict” function as “digestive_system.” We compared the difference of IC50, and the results were displayed in the form of boxplots using R package “ggpubr.”

External Validation of the FRLP Risk Score Model

Two microarray datasets for GC patients, GSE84437 and GSE62254, were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) database and used for external validation of the FRLP risk score model. The clinical information of patients from these two cohorts is also listed in Table 2. An alternative model to the FRLP risk score, the GS model was established and used to conduct the external validation (Zhang et al., 2021). Differentially expressed genes (DEGs) between high- and low-risk score groups were screened out based on |logFC|>0.585 (|FC|>1.5), p < 0.05 criteria. Gene set A comprised DEGs upregulated in the high-risk group, while low-risk group’s highly expressed DEGs were classified into gene set B. We performed ssGSEA to calculate the enrichment score (ES) of gene sets A and B in each sample. Consequently, the gene set score (GS), namely the subtraction of gene set B ES from gene set A ES, was estimated and used to represent the risk score when lncRNA expression profiles were not available. We performed the Spearman correlation analysis to show the relevance of the two scores (Zhang et al., 2021). Using the cut-off of the ROC curves for GS model, we classified patients into high- and low-risk groups. Then, we compared the patients’ stratification between the FRLP risk score model and the alternative GS model. A Sankey diagram was drawn to demonstrate the intersection of the two models using R package “Director” (Icay et al., 2018; Xiang et al., 2021; Tang et al., 2021). Finally, we calculated the GS of each sample in two GEO cohorts and stratified the patients using the same cut-off value mentioned above. We analyzed the differences in survival, TIDE score, and drug sensitivity between high- and low-risk groups. Additionally, a nomogram based on the GS model was constructed to enhance the clinical utility.

Development of Dynamic Nomograms Based on FRLP Risk Score Model and GS Model

To further enhance the clinical utility of the two aforementioned nomograms, we sought to visualize them in an interactive form. Using R package “DynNom” (Jalali et al., 2019; Yin et al., 2022), we generated two dynamic nomograms with interactive interfaces for clinical application based on the FRLP risk score model and GS model, respectively.

Statistical Analysis

All statistical analyses were performed in R version 4.1.2. Wilcoxon test was performed to conduct numerical difference comparisons of two groups. Log-rank test was performed to evaluate the differences in survival, and Kaplan–Meier plots were drawn to visualize the comparison. Univariate and multivariate Cox proportional hazard regression analyses were used to assess the predictive efficacy of the risk score. Statistical significance was set at p < 0.05.

Results

Identification of FRLPs in TCGA Cohort

The entire analytical process of the study is presented in Supplementary Figure S1. Our study included the transcriptome data of 375 STAD samples and 32 normal samples from the TCGA database and identified 14,086 lncRNAs. A total of 382 ferroptosis-related genes were downloaded from the FerrDb, and their expression was analyzed in the TCGA STAD expression matrix. We conducted the Pearson correlation analysis and screened out 1,343 ferroptosis-related lncRNAs. Then, 112 differentially expressed lncRNAs (DELs) were filtered out according to the |logFC|>2 criteria, of which 92 were upregulated and 20 were downregulated in tumor samples, respectively (Figures 1A,B). After pairing analysis of the 112 DELs, a 0 or 1 matrix of 4,722 FRLPs was established. Finally, the univariate Cox analysis identified 47 prognostic FRLPs (Figure 1C).

FIGURE 1

FIGURE 1. Identification of differentially expressed ferroptosis-related lncRNAs (DElncRNAs) and differentially expressed ferroptosis-lncRNA pairs (FRLPs). (A) Heat map of 112 DElncRNAs. Red dots and blue dots represent upregulated and downregulated lncRNAs in tumor samples, respectively. (B) Volcano map of 112 DElncRNAs. (C) Forest plot showing 47 prognostic FRLPs identified by univariate Cox analysis.

Construction of FRLP Risk Score Model

To further explore the prognostic value of FRLPs in GC, a stepwise Cox regression procedure was performed. Firstly, to reduce the overfitting among 46 prognostic FRLPs, we conducted a LASSO Cox regression analysis and 27 FRLPs were filtered out according to the minimum partial likelihood deviance (Figures 2A,B). Then, we performed multivariate Cox analysis on these 27 FRLPs and finally obtained 10 prognostic FRLPs (Figure 2C). The list of 10 FRLPs and their corresponding calculation coefficients are shown in Table 3. As VIFs were all <2, there wasn’t any multicollinearity relationship between these FRLPs (Table 4). Compared with other clinical factors, our risk score model was more accurate to predict the overall survival according to multi-factor ROC analysis (Figure 2D). Then, we conducted DCA to explore the clinical significance of the risk score. As shown in Figure 2E, the risk score displayed better net benefit than other factors, which indicated that the risk score was competent to help clinicians make more accurate assessment of patient prognosis compared with the age-, gender-, grade-, location-, or stage-only model. Moreover, the time-dependent ROC analysis revealed that the area under curves (AUCs) for 1-, 3-, and 5-year overall survival (OS) were 0.738, 0.795, and 0.787, respectively (Figure 2F).

FIGURE 2

FIGURE 2. Construction and evaluation of the FRLP risk score model in TCGA cohort. (A) Least absolute shrinkage and selection operator (LASSO) coefficients of 27 prognosis-related FRLPs. (B) Tenfold cross-validation for tuning parameter selection in the LASSO model. (C) Forest plot showing 10 FRLPs identified by multivariate Cox regression analysis. (D) Receiver operating characteristic (ROC) comparing the risk score and other clinical factors in predicting total OS. (E) Decision curve analysis (DCA) curves estimating the predictive efficacy of the risk score from the perspective of clinical benefit. The y-axis refers to the net benefit. The x-axis refers to the predicted OS. The black line represents the hypothesis that all patients survive in 5 years. The gray line represents the hypothesis that no patients stay alive for more than 1 year. (F) ROC curve for predicting 1-, 3-, and 5-year overall survival (OS) of the FRLP risk score model. (G) Cut-off point of the risk score model. (H) Kaplan–Meier plot of high- and low-risk patients. (I) The risk score distribution. Green dots represent risk scores for low-risk patients; red dots represent risk scores for high-risk patients. (J) The relationship between survival status and risk score. The horizontal ordinate represents the number of patients; the vertical ordinate represents risk score (AUC: area under curve; CI: confidence interval).

TABLE 3

TABLE 3. lncRNA pairs used for the construction of the FRLP risk score model.

TABLE 4

TABLE 4. The variance inflation factors (VIFs) of 10 lncRNA pairs.

Predictive Assessment of the FRLP Risk Score Model

We identified the maximum inflection point of 1.925 as the optimal cut-off point on the 3-year ROC curve (Figure 2G). Subsequently, 148 patients were classified into the low-risk group and 186 into high-risk group. The OS of the low-risk group was significantly longer than that of the high-risk group (Figure 2H). The distribution of patients’ risk score is depicted in Figure 2I. Based on the scatter plot, the number of deaths increased as the risk score increased (Figure 2J). To figure out whether the risk score model is independent of other clinicopathological parameters, univariate and multivariate Cox analyses incorporating age, gender, tumor location, and tumor stage were adopted. It was revealed that the FRLP model is an independent prognostic factor for STAD patients (univariate Cox p < 0.001, HR = 1.265, 95% CI 1.195–1.339; multivariate Cox p < 0.001, HR = 1.241, 95% CI: 1.171–1.315; Figures 3A,B). Subsequently, we built a nomogram that integrated tumor stage and risk score to predict patients’ 1-, 3-, and 5-year OS (Figure 3C). The AUCs of the nomogram for predicting 1-, 3-, and 5-year OS rates were 0.768, 0.775, and 0.766, respectively (Figure 3D). According to Schoenfield’s residual test, the individual p-value for age, stage, and risk score was 0.97, 0.80, and 0.80, respectively, while the global p-value was 0.95 (Figure 3E). These results indicated that the nomogram model met the equally proportional risk hypothesis. The total point of the nomogram (nomoscore) had higher efficacy in predicting GC patients’ 1, 3, and 5-year OS than age or stage (Figure 3F). Moreover, DCA analysis also revealed that the nomoscore displayed better net benefit than the age- or stage-alone model (Figure 3G). The calibration curve also confirmed the predictive value of the nomogram in predicting 1-, 3-, and 5-year OS (Figure 3H).

FIGURE 3

FIGURE 3. Evaluation of FRLP risk score model’s predictive efficacy. (A) Univariate Cox analysis of risk score and other clinicopathological factors. (B) Multivariate Cox analysis of risk score and other clinicopathological factors. (C) Nomogram integrating risk score, age, and tumor stage for predicting 1-, 3-, and 5-year OS. (D) Time-dependent ROC curve of the nomogram for predicting 1-, 3-, and 5-year OS. (E) Schoenfield’s residual test of age, stage, and risk score. (F) Calibration curves of the nomogram for predicting 1-, 3-, and 5-year OS. The gray lines represent the ideal predictive model, and the red lines represent the observed model. (G) Time-dependent ROC curves evaluating the efficacy of the nomogram to predict 1-, 3-, and 5-year OS. (H) DCA curves estimating the predictive efficacy of the nomogram from the perspective of clinical benefit.

The Clinical Significance of the FRLP Risk Score Model

To evaluate the relationship risk score and clinicopathological characteristics, chi-square test was performed and the results are demonstrated in Figure 4A, which indicated that tumor grade and tumor stage were closely linked to the FRLP risk level. In addition, we analyzed the risk score differences between groups stratified by different clinicopathological factors. As shown in Figures 4B–H, a high-risk score with statistical significance was more common to see in patients with higher tumor grades as well as more advanced N stages and TNM stages. Nevertheless, patients in different gender, age, M status, and T status groups exhibited no differences in risk scores.

FIGURE 4

FIGURE 4. Clinical correlation of the FRLP risk score model. (A) Heatmap showing the clinical relevance of the risk score model (*p < 0.05). (B–H) Boxplots showing risk score differences in different age, gender, tumor grade, tumor stage, T status, N status, and M status groups.

The Correlation of the FRLP Risk Score Model and Somatic Variance

To investigate the gene mutation landscape of the FRLP risk score model, we performed TMB and gene mutation frequency analysis. As shown in Figure 5A, low-risk STAD patients had higher TMB levels than high-risk patients. And a negative correlation between TMB levels and risk score was found according to the Spearman correlation analysis (R = −0.17, p = 0.0016) (Figure 5B). The Kaplan–Meier analysis suggested that STAD patients with high TMB had longer OS than low TMB patients (Figure 5C). Stratified survival analysis further confirmed that the survival benefits for high TMB patients still exist in both high- and low-risk groups (Figure 5D). In general, gene mutation frequency was higher in the low-risk group than that in the high-risk group (Figures 5E,F).

FIGURE 5

FIGURE 5. Tumor mutation burden (TMB) analysis of the FRLP risk score model. (A) TMB difference between high- and low-risk groups. (B) Correlation between the risk score and TMB. (C) Kaplan–Meier plots of patients with high and low TMB. (D) Kaplan–Meier curves of patients stratified by both TMB and the risk score. (E,F) Gene mutation analysis of patients in low- and high-risk groups.

The Correlation of the FRLP Risk Score Model and Tumor Microenvironment

To shed light on the model’s association with tumor microenvironment, we estimated the stromal scores, immune scores, and ESTIMATE scores of two risk groups by the R package “ESTIMATE”. These three scores were reflections of stromal contents, immune cell contents, and the aggregation of the two contents. Then, we calculated the content of tumor infiltrating immune cells (TIICs) using diversified algorithms online and discussed the correlation between the risk score and tumor-infiltrating immune cells via Spearman correlation analysis. The results demonstrated that the stromal score and ESTIMATE score were higher in the high-risk group, while there was no significant difference in the immune score between the two groups (Figure 6C). High-risk scores were more closely linked with high TIICs (Figure 6A). We specified the MCPcounter algorithm result, and it was indicated that immunosuppressive cells (e.g., endothelial cells, fibroblasts, monocytes, and neutrophils) were more densely infiltrated in the high-risk group, whereas there was no statistical significance of anti-tumor immune cells (e.g., B lineage, CD8 T cells, cytotoxic lymphocytes, NK cells, and T cells) infiltration between the two groups (Supplementary Figure S1).

FIGURE 6

FIGURE 6. Tumor infiltrating immune cells (TIICs) and immunotherapeutic sensitivity analysis of the FRLP risk score model. (A) The correlation between risk score and TIICs analyzed by seven different quantification methods of immune infiltration estimations including TIMER, xCell, quanTIseq, MCP-counter, EPIC, CIBERSORT-ABS, and CIBERSORT. (B) Expression of 24 immune checkpoint genes in high- and low-risk groups. (C) Boxplots showing immune score, stromal score, and ESTIMATE score in high- and low-risk groups. (D) Boxplots showing dysfunction score, exclusion score, and tumor immune dysfunction and exclusion (TIDE) score differences between high- and low-risk score groups. (E) Immunophenoscore (IPS) differences for ICB treatment between high- and low-risk groups; ips_ctla4_neg_pd1_pos refers to CTLA4-negative response and PD1-positive response; ips_ctla4_pos_pd1_neg refers to CTLA4-positive response and PD1-negative response; ips_ctla4_pos_pd1_pos refers to CTLA4-positive response and PD1-positive response (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001).

The Correlation of the FRLP Risk Score Model and Immunotherapeutic Sensitivity

Of all the 24 selected ICGs, CD28, CD86, FGL2, HAVCR2, PDCD1LG2, TNFSF4, and TNFSF18 were highly expressed in the high-risk group (Figure 6B). These results further confirmed the immunosuppressive phenotype of the high-risk patients. Intriguingly, according to the TIDE analysis, which suggested that high-risk patients had a heightened level of dysfunction score, exclusion score, and TIDE score (Figure 6D), high-risk patients might not actually benefit from ICB treatment though highly expressed in ICGs. On the contrary, the IPS scores of anti-CTLA4+ anti-PD1+, anti-CTLA4- anti-PD1+, and anti-CTLA4+ anti-PD1- were all higher in low-risk subgroups (Figure 6E), implying low-risk patients’ better responses toward anti-CTLA4 and/or anti-PD1 immunotherapy.

The Correlation of the FRLP Risk Score Model and Biological Function

To investigate different risk groups’ enriched biological function, GSEA using “KEGG” and “GO” gene sets was conducted to compare the enrichment differences between two risk groups. As shown in Figures 7A,B, KEGG pathways in relation to stromal activity and diseases including “ECM–receptor interaction,” “focal adhesion,” and “dilated cardiomyopathy” were enriched in the high-risk group. Meanwhile, GO items in relation to cell migration such as “ameboidal-type cell migration,” “cell junction assembly,” and “cell matrix adhesion” were enriched in the high-risk group. For the low risk group, we found an enrichment tendency toward KEGG metabolism-related pathways (e.g., “nitrogen metabolism,” “ribosome,” and “peroxisome”) and GO metabolism-related items (e.g., “rRNA metabolic process,” “mitochondrial-protein containing complex,” and “structural constituent of ribosome”). In addition, we performed ssGSEA using the “Hallmark” gene sets. As revealed by Figure 7C, the high-risk group was markedly correlated with carcinogenic activities (e.g., “Wnt beta catenin,” “myc targets,” “kras signaling,” etc.) and stromal pathways (e.g., “hypoxia,” “angiogenesis,” and “epithelial–mesenchymal transition”), whereas the low-risk group was characterized by its enrichment in cell cycle-related events (e.g., “DNA repair,” “G2M checkpoint,” and “E2F targets”).

FIGURE 7

FIGURE 7. Biological function analysis of the FRLP risk score model. (A) Gene set enrichment analysis (GSEA) analysis using KEGG gene sets for high- and low-risk groups. (B) GSEA analysis using GO gene sets for high- and low-risk groups. (C) ssGSEA analysis using Hallmark gene sets for high- and low-risk groups. Gene sets markedly enriched in high- or low-risk groups were marked red and blue, respectively (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ns: not significantly).

External Validation of the FRLP Risk Score Model

The lncRNAs in our FRLP risk score model could not be detected in two GEO microarray datasets due to limited sequencing depth. Therefore, we constructed gene set A and gene set B, which were composed of genes significantly upregulated in high- and low-risk groups, respectively (Supplementary Figure S3A, Supplementary Table S1). KEGG analysis suggested that the two gene sets targeted at different pathways, with gene set A associated with stromal and carcinogenic pathways and gene set B targeted at immune-related pathways. (Supplementary Figures S3B,C). Then, we developed a novel gene set score (GS) model that was defined as the subtraction of gene set B enrichment score (ES) from gene set A ES. The distribution of GS and the ES of two gene sets was significantly between high- and low-FRLP risk score groups, and the high-FRLP risk score group also exhibited a higher GS score (Supplementary Figure S3D). Moreover, there was a positive correlation between risk score and GS score (Spearman correlation R = 0.57; Supplementary Figure S3E). Subsequently, patients were divided into high- and low-GS groups according to the cut-off value of 3-year ROC curves, which was identical to the method we used to construct the FRLP risk score model (Supplementary Figure S3F). There was a significant overlap between the established FRLP risk score model and the GS model: over 70% of the high-FRLP risk score patients could be classified into a high-GS group, and the proportion of low-FRLP risk score patients falling to the category of low-GS group is over 75% (Supplementary Figures S3G,H). Since the GS model had a high degree of compliance with the FRLP risk score model, it is reasonable to regard the GS score as an alternative to FRLP risk score in distinguishing GC patients’ survival, biological function, tumor microenvironment, etc.

We calculated the GS of each sample in two GEO cohorts, and the classification of distinct risk groups was also conducted using the same cut-off value in the TCGA cohort (Supplementary Figure S3D). High-risk patients were validated to have worse OS in GSE84437 (n = 433, log rank p < 0.001), GSE62254 (n = 300, log rank p < 0.001), and shortened disease-free survival (DFS) in GSE62254 (log rank p < 0.001) (Figures 8A–C). Time-dependent ROC analysis also identified the efficacy of the alternative GS model in predicting patients’ survival. In GSE62254, the AUCs for 1-, 3-, and 5-year’s OS were 0.681, 0.660, and 0.661, respectively; the AUCs for 1-, 3-, and 5-year’s DFS were 0.632, 0.660, and 0.701, respectively. In GSE84437, the AUCs for 1-, 3-, and 5-year’s OS were 0.578, 0.612, and 0.630, respectively (Figures 8D–F). Specifically, in GSE62254, stage III/IV patients tended to exhibit higher GS (Figure 8G). Moreover, GSE62254 contains information of two classification systems for GC patients: the Asian Cancer Research Group (ACRG) subtype (Cristescu et al., 2015) and Lauren subtype. By inspecting the intersection of our GS model and these two subtypes, we found that there were a larger proportion of high-risk patients classified into “EMT” ACRG subtype and “diffuse” Lauren subtype (Figures 8H,I), which were both subtypes indicating worse clinical outcomes. In addition, TIDE analysis also showed elevated dysfunction, exclusion, and TIDE scores in high-risk groups in both of the two cohorts (Supplementary Figure S4).

FIGURE 8

FIGURE 8. External validation of the FRLP risk score model in two GEO cohorts using the alternative GS model. Kaplan–Meier plots comparing OS of high- and low-risk patients in GSE84437 (A) and GSE62254 (B). (C) Kaplan–Meier plots comparing disease-free survival (DFS) of high- and low-risk patients in GSE62254. ROC analysis of the GS model in predicting 1, 3, and 5 year OS in (D) GSE84437 and (E) GSE62254. (F) ROC analysis of the GS model in predicting 1, 3, and 5 year DFS in GSE62254. (G) Boxplots comparing GS score differences among different stages of patients in GSE62254. (H,I) Barplots showing the relative proportion of four Asian Cancer Research Group (ACRG) subtypes and three Lauren subtypes in high- and low-risk patients in GSE62254.

The Correlation of the FRLP Risk Score Model or GS Model and Chemotherapy Sensitivity

The ssGSEA analysis revealed a higher enrichment level of carcinogenic pathways in the high-risk group. This indicated that high-risk patients may exhibit better responses toward chemotherapy. We estimated the IC50 of six anti-tumor drugs in samples from the TCGA cohort and two GEO cohorts. For TCGA patients, we compared the sensitivity between high- and low-risk groups of the FRLP risk score model. For GEO patients, we compared the sensitivity between high- and low-risk groups of the alternative GS model. As demonstrated in Figure 9, the sensitivity to imatinib, bosutinib, vinblastin, doxorubicin, and cisplatin was significantly higher in the high-risk patients than that of the low-risk patients. On the other hand, low-risk patients only displayed higher sensitivity to methotrexate.

FIGURE 9

FIGURE 9. Chemotherapeutic response analysis. Boxplots comparing differences in half-inhibitory concentration (IC50) values of six anti-tumor drugs between high- and low-risk score groups in (A) TCGA, (B) GSE84437, and (C) GSE62254 (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001).

Construction of a Nomogram Based on GS Model

Considering the inconvenience of using GS score directly in predicting patients' prognosis, we constructed a nomogram integrating GS score, age, and tumor stage based on TCGA cohort and GSE62254 (GSE84437 was excluded for lack of tumor stage profile) (Figure 10A). AUC analysis (Figures 10B,C) and calibration curves (Figures 10D,E) confirmed the high accuracy of the nomogram for predicting OS at 1, 3, and 5 years in both of the two cohorts. Furthermore, time-dependent ROC (Figures 10F,G) and DCA analysis (Figures 10H,I) was also performed in each of the cohort, which all confirmed the superior predicting ability of the nomogram compared with the age- or stage-only model.

FIGURE 10

FIGURE 10. Construction of a nomogram using GS model. (A) Nomogram integrating GS score, age, and tumor stage for predicting 1, 3, and 5 year OS in TCGA and GSE62254. 1, 3, and 5 year ROC analyses of the nomogram in (B) TCGA and (C) GSE62254. Calibration curves of the nomogram for predicting 1, 3, and 5 year OS in (D) TCGA and (E) GSE62254. Time-dependent ROC curves evaluating the efficacy of the nomogram to predict 1, 3, and 5 year OS in (F) TCGA and (G) GSE62254. DCA curves estimating the predictive efficacy of the nomogram in (H) TCGA and (I) GSE62254.

Development of Dynamic Nomograms Based on FRLP Risk Score Model and GS Model

At the end of the study, we generated two online dynamic nomograms based on FRLP risk and GS score model, respectively (https://ljzwhdx.shinyapps.io/FRLPdynanomo/; https://ljzwhdx.shinyapps.io/GSdynanomo/). Both of the two nomograms were in interactive forms, which could facilitate clinician’s prediction for prognosis. Supplementary Figure S5 displays the interface of the nomogram based on the FRLP risk score model. Predicted values of time-dependent survival probability could be easily obtained after selecting risk score or GS score, stage, and age of a specific GC patient.

Discussion

Ferroptosis is a form of programmed cell death intimately associated with iron metabolism and peroxidation of polyunsaturated fatty acids (Stockwell et al., 2017). LncRNAs were able to regulate gene expression at both transcriptional and post-transcriptional levels (Statello et al., 2021). In recent years, the construction of prognostic prediction models based on ferroptosis-related lncRNAs has attracted particular attention of researchers (Table 1). However, these risk models concentrated on the specific expression of lncRNAs, which limited their clinical practicability. Two recently established risk score models independent of the exact lncRNA expressions caught our attention (Li et al., 2021; Tang et al., 2021). In the present study, we demonstrated that the similar risk score model of ferroptosis-related lncRNA pairs (FRLP) could also be applied to GC patients. The AUCs for the risk score to predict 1, 3, and 5 years’ OS were all over 0.70 according to the ROC curves. Aside from the high accuracy in predicting GC patients’ clinical outcomes, the FRLP risk score model was also closely connected with tumor microenvironment, biological function, and responses to chemo- and immunotherapies. High-risk patients of our FRLP risk score model were characterized by higher infiltration of immune cells (especially immunosuppressive and pro-tumorigenic cells), higher carcinogenic and stromal activities, and better sensitivity to some types of anti-tumor drugs. On the other hand, low-risk patients displayed better treatment responses to methotrexate and higher immunotherapy sensitivity. Considering that lncRNA expression profiles in the FRLP model were not available in microarray datasets, we also introduced an alternative gene set score (GS) model. We performed a series of analyses to confirm the high degree of compliance between FRLP and GS models. Subsequently, we validated the robustness of our FRLP risk score model in two external GEO cohorts using the alternative GS score, which confirmed the potential of applying the risk score model to a wider range of patients. Moreover, to enhance the clinical utility of the FRLP and the alternative GS risk score model, we also built two nomograms based on the two models, respectively. Finally, two dynamic nomograms with interactive interfaces based on the FRLP model and GS model were also constructed to further facilitate clinician’s prediction for prognosis.

The interaction between tumor cells and tumor infiltrating immunes (TIICs), namely the tumor microenvironment (TME), has been proved to play a crucial role in tumorigenesis. Some types of infiltrating immune cells have anti-tumor activities, such as CD8⁺ cytotoxic T lymphocytes (CTL), CD4⁺ T helper cells, natural killer cells, and dendritic cells. Conversely, regulatory T cells (Treg) and myeloid-derived suppressor cells (MDSCs) are regarded as “bad guys” in the TME (Pitt et al., 2016). Besides, some immune cells undergo polarization in the development of cancer and can exhibit either anti-tumor or tumor-promoting function depending on different cancer stages, such as neutrophils and macrophages (Ngambenjawong et al., 2017; Giese et al., 2019). Some works still regard these two kinds of cells tumor-promoting as their high infiltration level was frequently correlated with poor prognosis of many human tumor types (Mantovani et al., 2017; Shaul and Fridlender, 2019). Cancer-associated fibroblasts (CAFs) are another important kind of immunosuppressive cells that promote tumor progression by modulating angiogenesis and epithelial–mesenchymal transmission of cancer cells (Sahai et al., 2020).

Currently, the relationship between ferroptosis and immunity is quite ambiguous. Experimental evidence has already shown that interferon gamma (IFN-γ) released by CD8⁺ T cells promotes tumor cell ferroptosis (Wang et al., 2019). The negative association between CD8⁺ T cell infiltration, IFN-γ expression, and the expression of SLC3A2 and SLC7A11 was also shown in human melanoma tissues (Wang et al., 2019). A recent study revealed that CAFs exerted the tumor promoting role in gastric cancer by inhibiting ferroptosis in GC cells (Zhang et al., 2020). However, it is still immature to classify ferroptosis into the category of “immune-related cell death.” In the current study, we identified that the TME of high-risk group was characterized by higher stromal contents, as demonstrated by ESTIMATE and biological function analysis: high-risk patients exhibited higher stromal scores, and KEGG stromal-related pathways including “ECM–receptor interaction,” “focal adhesion,” “dilated cardiomyopathy,” “epithelial–mesenchymal transition,” “hypoxia,” and “angiogenesis” had a higher enrichment level in the high-risk group. As for the immune cells in TME, despite a positive correlation between the risk score and TIIC infiltration, there were no significant differences in the anti-tumor immune cell content between the two risk groups. On the other hand, pro-tumorigenic cells were more densely infiltrated in the TME of high-risk patients. It is also noteworthy that “GO” items related to cell migration and “Hallmark” items related to oncogenic activities were also enriched in the high-risk group. Thus, we speculated that the TME of the high-risk group was stroma-related and immunosuppressive, which may closely be linked with tumorigenesis.

Immune checkpoint is capable of inhibiting the over-activation of T cells and preventing autoimmune diseases. But under tumor circumstances, it will prevent T cells from approaching the tumor, weakening the ability of the immune system to recognize and destroy tumor cells (Tan et al., 2020). In recent years, immunotherapy targeting immune checkpoint modulation, namely the immune checkpoint blockage (ICB), has shown promising efficacy in cancer treatment (de Miguel and Calvo, 2020). However, the benefit, to date, has been limited to a minority of patients with certain cancer types (Sharma et al., 2017). TMB and TIDE scores are two effective methods to evaluate responses to ICB treatment. TMB is a quantitative measure of the total number of somatic non-synonymous mutations per coding area of a tumor genome (Melendez et al., 2018). Generally, high TMB suggests better OS in cancer patients after receiving ICB (Cao et al., 2019; Samstein et al., 2019; Valero et al., 2021). TIDE score was a computational method that combines two primary mechanisms in tumor evasion, T cell dysfunction and T cell exclusion. A higher tumor TIDE score is associated not only with worse ICB response but also with worse patient survival under anti-PD1 and anti-CTLA4 therapies (Jiang et al., 2018). In the present study, we observed a heightened expression of seven ICGs in the high-risk group. To figure out which group of GC patients displayed better ICB treatment response, we performed TMB and TIDE analyses in two groups of patients. Intriguingly, we found that high-risk STAD patients exhibited lower TMB and higher TIDE, which indicated that high-risk STAD patients may not actually benefit from ICB treatment though highly expressed seven ICGs. On the contrary, patients in low-risk groups may exhibit better immunotherapeutic response due to the relatively high TMB and low TIDE scores. IPS data collected from TCIA database also demonstrated that TCGA low-risk score or GEO low-GS patients were more likely to benefit from anti-CTLA4 and (or) anti-PD1 immunotherapy.

Finally, high-risk patients were demonstrated to exhibit higher sensitivity toward five kinds of chemotherapeutic or targeted drugs including vinblastin, cisplatin, doxorubicin, immatinib, and bosutinib. In contrast, low-risk patients solely responded to methotrexate. This perhaps could be attributed to two group’s discrepancy in ssGSEA analysis. The overactivity of carcinogenic pathways in high-risk patients, on the one hand, contributed to the shortened survival. On the other hand, however, this may also confer high-risk patients exhibiting better responses to a wider range of anti-tumor drugs. For instance, protein tyrosine kinase-associated pathways, such as IL6/JAK/STAT and PI3K/AKT, were highly enriched in high-risk patients, which may account for high-risk patients’ elevated sensitivity to tyrosine kinase inhibitors (TKIs) such as immatinib and bosutinib. Item “E2F target” was markedly enriched in the low-risk group according to ssGSEA analysis. Transcriptional factor E2F family plays a critical role in determining the time point of cell division. The expression of E2F targets gradually increases during G1 and must reach a threshold level in order for cells to pass the restriction point and progress to S phase. It was revealed that the activity of E2F peaks at G1-S transition gradually decreases during the phase and is totally repressed in both G2 and M phases (Kent and Leone, 2019). Therefore, the above analysis probably could explain why low-risk patients exhibited higher sensitivity to G1/S specific drug, methotrexate. Besides, the proportion of tumor cells arrested in the G2 phase was higher in low-risk patients due to higher enrichment level of “G2M checkpoint.” Hence, low-risk patients were less sensitive to the M phase-specific drug vinblastin. Our findings may help to optimize current chemotherapeutic strategy for GC patients.

There are several limitations of our present study. First, limited by the sequencing depth of GEO dataset, our external validation of the FRLP risk score signature used an alternative GS score, instead of the risk score, which may affect the robustness of the validation results. Second, all of our conclusions were drawn based on public databases; future large-scale and real-world studies are thus warranted.

Conclusion

In this study, we developed an FRLP risk model composed of 10 differentially expressed ferroptosis lncRNA pairs that does not rely on the exact lncRNA expression level in the TCGA cohort. An alternative GS model, which shared a high degree of compliance with the FRLP model, was also constructed and used it to validate the application value of the FRLP model in two GEO microarray datasets. Two online dynamic nomograms with interactive interfaces were finally generated to facilitate clinician’s prediction for prognosis. The novel FRLP risk score signature established in the current study might provide insights for the accurate prediction and comprehensive management for GC patients.

Data Availability Statement

The datasets supporting the conclusions of this article are available in The Cancer Genome Atlas (https://portal.gdc.cancer.gov), Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo) and University of California Santa Cruz database (https://xena.ucsc.edu). The two online nomograms generated in the study could be accessed via https://ljzwhdx.shinyapps.io/FRLPdynanomo/ and https://ljzwhdx.shinyapps.io/GSdynanomo/.

Author Contributions

TF designed and supervised the project. JL and CK curated and analyzed the data. JL, RX, and JW interpreted and visualized the results. JL, RX, and WS drafted and reviewed the article. All authors read and approved the final article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.899419/full#supplementary-material

Supplementary Figure S1 | Analytic process of the study.

Supplementary Figure S2 | Violin plots exhibiting the infiltration differences of 10 immune cell types between high- and low-risk groups.

Supplementary Figure S3 | Construction of the gene set score (GS) model, an alternative to the FRLP risk score model in TCGA cohort. (A) Volcano plot displaying genes upregulated and downregulated in high-risk patients, which were defined as gene sets A and B, respectively. (B,C) KEGG analysis of gene set A and gene set B. (D) Boxplots showing the differences in FRLP risk score, enrichment score of gene sets A and B, GS score in TCGA high- and low-FRLP risk groups. (E) Spearman correlation analysis of FRLP risk score and GS score. (F) Cut-off of the GS model. (G) Barplot and (H) Sankey diagram showing GS model’s representing ability of FRLP risk score model (*** p<0.001, **** p<0.0001).

Supplementary Figure S4 | Boxplots showing differences in dysfunction score, exclusion score, and TIDE score of high- and low-GS groups in (A) GSE84437 and (B) GSE62254 (**** p<0.0001).

Supplementary Figure S5 | The interface of online dynamic nomogram (https://ljzwhdx.shinyapps.io/FRLPdynanomo/) integrating FRLP risk score, tumor stage, and age for predicting time-independent survival probabilities in TCGA. (A) Input area for users to select stage (stages I–IV) or age (>65 or ≤65) as well as input the risk score and the follow-up time (futime). (B) Survival plots, showing patients' survival probabilities at different time points. (C) Predicted survival probabilities with a 95% confidence interval (CI), which could be obtained after inputting patient information in the input area. For example, when selecting “>65” for age, “stage III” for stage, and entering “3” for risk score, then a patient's 1-, 3-, and 5-year survival probabilities with a 95% CI are displayed using a black line, a blue line, and a red line, respectively. (D) Numerical summary showing the exact values of survival probabilities with a 95% CI.

Supplementary Table S1 | Lists of genes in gene set A and gene set B.

Abbreviations

AUC, area under the curve; DEL, differentially expressed lncRNA; DFS, disease-free survival; FRLP, ferroptosis-related lncRNA pair; GC, gastric cancer; GSEA, gene set enrichment analysis; GEO, Gene Expression Omnibus; GO, gene ontology; IC50, half-inhibitory concentration; IPS, immunophenoscore; ICG, immune checkpoint gene; ICB, immune checkpoint blockage; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; ssGSEA, single sample gene set enrichment analysis; STAD, stomach adenocarcinoma; TCGA, the cancer genome database; TME, tumor microenvironment; TIIC, tumor infiltrating immune cell; TIDE, tumor immune dysfunction and exclusion; TCIA, The Cancer Immunome Atlas; UCSC, University of California Santa Cruz.

References

Behan, F. M., Iorio, F., Picco, G., Gonçalves, E., Beaver, C. M., Migliardi, G., et al. (2019). Prioritization of Cancer Therapeutic Targets Using CRISPR-Cas9 Screens. Nature 568 (7753), 511–516. doi:10.1038/s41586-019-1103-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhan, A., Soleimani, M., and Mandal, S. S. (2017). Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res. 77 (15), 3965–3981. doi:10.1158/0008-5472.can-16-2634

PubMed Abstract | CrossRef Full Text | Google Scholar

Binatti, A., Bresolin, S., Bortoluzzi, S., and Coppe, A. (2021). iWhale: a Computational Pipeline Based on Docker and SCons for Detection and Annotation of Somatic Variants in Cancer WES Data. Brief. Bioinform 22 (3), bbaa065. doi:10.1093/bib/bbaa065

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, S., Fu, S., Zhang, W., Yuan, X., Cheng, Y., and Fang, J. (2021). SIRT6 Silencing Overcomes Resistance to Sorafenib by Promoting Ferroptosis in Gastric Cancer. Biochem. Biophysical Res. Commun. 577, 158–164. doi:10.1016/j.bbrc.2021.08.080

CrossRef Full Text | Google Scholar

Cao, D., Xu, H., Xu, X., Guo, T., and Ge, W. (2019). High Tumor Mutation Burden Predicts Better Efficacy of Immunotherapy: a Pooled Analysis of 103078 Cancer Patients. Oncoimmunology 8 (9), e1629258. doi:10.1080/2162402x.2019.1629258

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Qiao, L., Bian, Y., and Sun, X. (2020). GDF15 Knockdown Promotes Erastin-Induced Ferroptosis by Decreasing SLC7A11 Expression. Biochem. Biophysical Res. Commun. 526 (2), 293–299. doi:10.1016/j.bbrc.2020.03.079

CrossRef Full Text | Google Scholar

Chen, W., Zheng, R., Baade, P. D., Zhang, S., Zeng, H., Bray, F., et al. (2016). Cancer Statistics in China, 2015. CA A Cancer J. Clin. 66 (2), 115–132. doi:10.3322/caac.21338

PubMed Abstract | CrossRef Full Text | Google Scholar

Coutzac, C., Pernot, S., Chaput, N., and Zaanan, A. (2019). Immunotherapy in Advanced Gastric Cancer, Is it the Future? Crit. Rev. Oncology/Hematology 133, 25–32. doi:10.1016/j.critrevonc.2018.10.007

CrossRef Full Text | Google Scholar

Cristescu, R., Lee, J., Nebozhyn, M., Kim, K.-M., Ting, J. C., Wong, S. S., et al. (2015). Molecular Analysis of Gastric Cancer Identifies Subtypes Associated with Distinct Clinical Outcomes. Nat. Med. 21 (5), 449–456. doi:10.1038/nm.3850

PubMed Abstract | CrossRef Full Text | Google Scholar

de Miguel, M., and Calvo, E. (2020). Clinical Challenges of Immune Checkpoint Inhibitors. Cancer Cell 38 (3), 326–333. doi:10.1016/j.ccell.2020.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, S. J. (2017). Ferroptosis: Bug or Feature? Immunol. Rev. 277 (1), 150–157. doi:10.1111/imr.12533

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, S. J., Lemberg, K. M., Lamprecht, M. R., Skouta, R., Zaitsev, E. M., Gleason, C. E., et al. (2012). Ferroptosis: an Iron-dependent Form of Nonapoptotic Cell Death. Cell 149 (5), 1060–1072. doi:10.1016/j.cell.2012.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedmann Angeli, J. P., Krysko, D. V., and Conrad, M. (2019). Ferroptosis at the Crossroads of Cancer-Acquired Drug Resistance and Immune Evasion. Nat. Rev. Cancer 19 (7), 405–414. doi:10.1038/s41568-019-0149-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Garris, C. S., Arlauckas, S. P., Kohler, R. H., Trefny, M. P., Garren, S., Piot, C., et al. (2018). Successful Anti-PD-1 Cancer Immunotherapy Requires T Cell-Dendritic Cell Crosstalk Involving the Cytokines IFN-γ and IL-12. Immunity 49 (6), 1148–1161. doi:10.1016/j.immuni.2018.09.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Giese, M. A., Hind, L. E., and Huttenlocher, A. (2019). Neutrophil Plasticity in the Tumor Microenvironment. Blood 133 (20), 2159–2167. doi:10.1182/blood-2018-11-844548

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., et al. (2019). Anti-tumour Immunity Controlled through mRNA m6A Methylation and YTHDF1 in Dendritic Cells. Nature 566 (7743), 270–274. doi:10.1038/s41586-019-0916-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, S., Yu, J., He, W., Huang, Q., Zhao, Y., Liang, B., et al. (2017). Cysteine Dioxygenase 1 Mediates Erastin-Induced Ferroptosis in Human Gastric Cancer Cells. Neoplasia 19 (12), 1022–1032. doi:10.1016/j.neo.2017.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Icay, K., Liu, C., and Hautaniemi, S. (2018). Dynamic Visualization of Multi-Level Molecular Data: The Director Package in R. Comput. Methods Programs Biomed. 153, 129–136. doi:10.1016/j.cmpb.2017.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jalali, A., Alvarez-Iglesias, A., Roshan, D., and Newell, J. (2019). Visualising Statistical Models Using Dynamic Nomograms. PLoS One 14 (11), e0225253. doi:10.1371/journal.pone.0225253

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24 (10), 1550–1558. e0225253. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Q., Chen, H., Tang, Z., Sun, J., Ruan, Y., Liu, F., et al. (2021). Stemness-related LncRNA Pair Signature for Predicting Therapy Response in Gastric Cancer. BMC Cancer 21 (1), 1067. doi:10.1186/s12885-021-08798-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamarudin, A. N., Cox, T., and Kolamunnage-Dona, R. (2017). Time-dependent ROC Curve Analysis in Medical Research: Current Methods and Applications. BMC Med. Res. Methodol. 17 (1), 53. doi:10.1186/s12874-017-0332-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kent, L. N., and Leone, G. (2019). The Broken Cycle: E2F Dysfunction in Cancer. Nat. Rev. Cancer 19 (6), 326–338. doi:10.1038/s41568-019-0143-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Koppula, P., Zhuang, L., and Gan, B. (2021). Cystine Transporter SLC7A11/xCT in Cancer: Ferroptosis, Nutrient Dependency, and Cancer Therapy. Protein Cell 12 (8), 599–620. doi:10.1007/s13238-020-00789-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Liu, L., Huang, T., Jin, M., Zheng, Z., Zhang, H., et al. (2021). Establishment of a Novel Ferroptosis-Related lncRNA Pair Prognostic Model in Colon Adenocarcinoma. Aging 13 (19), 23072–23095. doi:10.18632/aging.203599

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, C., Zhang, X., Yang, M., and Dong, X. (2019). Recent Progress in Ferroptosis Inducers for Cancer Therapy. Adv. Mater 31 (51), e1904197. doi:10.1002/adma.201904197

PubMed Abstract | CrossRef Full Text | Google Scholar

Luis, G., Godfroid, A., Nishiumi, S., Cimino, J., Blacher, S., Maquoi, E., et al. (2021). Tumor Resistance to Ferroptosis Driven by Stearoyl-CoA Desaturase-1 (SCD1) in Cancer Cells and Fatty Acid Biding Protein-4 (FABP4) in Tumor Microenvironment Promote Tumor Recurrence. Redox. Biol. 43 (51), 102006. doi:10.1016/j.redox.2021.102006

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantovani, A., Marchesi, F., Malesci, A., Laghi, L., and Allavena, P. (2017). Tumour-associated Macrophages as Treatment Targets in Oncology. Nat. Rev. Clin. Oncol. 14 (7), 399–416. doi:10.1038/nrclinonc.2016.217

PubMed Abstract | CrossRef Full Text | Google Scholar

Meléndez, B., Van Campenhout, C., Rorive, S., Remmelink, M., Salmon, I., and D'Haene, N. (2018). Methods of Measurement for Tumor Mutational Burden in Tumor Tissue. Transl. Lung Cancer Res. 7 (6), 661–667. doi:10.21037/tlcr.2018.08.02

PubMed Abstract | CrossRef Full Text | Google Scholar

Ngambenjawong, C., Gustafson, H. H., and Pun, S. H. (2017). Progress in Tumor-Associated Macrophage (TAM)-targeted Therapeutics. Adv. Drug Deliv. Rev. 114, 206–221. doi:10.1016/j.addr.2017.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishino, M., Ramaiya, N. H., Hatabu, H., and Hodi, F. S. (2017). Monitoring Immune-Checkpoint Blockade: Response Evaluation and Biomarker Development. Nat. Rev. Clin. Oncol. 14 (11), 655–668. doi:10.1038/nrclinonc.2017.88

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, J., Zhang, X., Fang, X., and Xin, Z. (2021). Construction on of a Ferroptosis-Related lncRNA-Based Model to Improve the Prognostic Evaluation of Gastric Cancer Patients Based on Bioinformatics. Front. Genet. 12, 739470. doi:10.3389/fgene.2021.739470

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, S. J., Sanjana, N. E., Kishton, R. J., Eidizadeh, A., Vodnala, S. K., Cam, M., et al. (2017). Identification of Essential Genes for Cancer Immunotherapy. Nature 548 (7669), 537–542. doi:10.1038/nature23477

PubMed Abstract | CrossRef Full Text | Google Scholar

Pitt, J. M., Marabelle, A., Eggermont, A., Soria, J.-C., Kroemer, G., and Zitvogel, L. (2016). Targeting the Tumor Microenvironment: Removing Obstruction to Anticancer Immune Responses and Immunotherapy. Ann. Oncol. 27 (8), 1482–1492. doi:10.1093/annonc/mdw168

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahai, E., Astsaturov, I., Cukierman, E., DeNardo, D. G., Egeblad, M., Evans, R. M., et al. (2020). A Framework for Advancing Our Understanding of Cancer-Associated Fibroblasts. Nat. Rev. Cancer 20 (3), 174–186. doi:10.1038/s41568-019-0238-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Samstein, R. M., Lee, C.-H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor Mutational Load Predicts Survival after Immunotherapy across Multiple Cancer Types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, P., Hu-Lieskovan, S., Wargo, J. A., and Ribas, A. (2017). Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell 168 (4), 707–723. doi:10.1016/j.cell.2017.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaul, M. E., and Fridlender, Z. G. (2019). Tumour-associated Neutrophils in Patients with Cancer. Nat. Rev. Clin. Oncol. 16 (10), 601–620. doi:10.1038/s41571-019-0222-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Statello, L., Guo, C.-J., Chen, L.-L., and Huarte, M. (2021). Gene Regulation by Long Non-coding RNAs and its Biological Functions. Nat. Rev. Mol. Cell Biol. 22 (2), 96–118. doi:10.1038/s41580-020-00315-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Stockwell, B. R., Friedmann Angeli, J. P., Bayir, H., Bush, A. I., Conrad, M., Dixon, S. J., et al. (2017). Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism, Redox Biology, and Disease. Cell 171 (2), 273–285. doi:10.1016/j.cell.2017.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Stockwell, B. R., Jiang, X., and Gu, W. (2020). Emerging Mechanisms and Disease Relevance of Ferroptosis. Trends Cell Biol. 30 (6), 478–490. doi:10.1016/j.tcb.2020.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. U.S.A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660

CrossRef Full Text | Google Scholar

Tan, S., Li, D., and Zhu, X. (2020). Cancer Immunotherapy: Pros, Cons and beyond. Biomed. Pharmacother. 124, 109821. doi:10.1016/j.biopha.2020.109821

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, R., Wu, Z., Rong, Z., Xu, J., Wang, W., Zhang, B., et al. (2021). Ferroptosis-related lncRNA Pairs to Predict the Clinical Outcome and Molecular Characteristics of Pancreatic Ductal Adenocarcinoma. Brief. Bioinform 23 (1), bbab388. doi:10.1093/bib/bbab388

CrossRef Full Text | Google Scholar

Valero, C., Lee, M., Hoen, D., Wang, J., Nadeem, Z., Patel, N., et al. (2021). The Association between Tumor Mutational Burden and Prognosis Is Dependent on Treatment Context. Nat. Genet. 53 (1), 11–15. doi:10.1038/s41588-020-00752-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Vickers, A. J., and Elkin, E. B. (2006). Decision Curve Analysis: a Novel Method for Evaluating Prediction Models. Med. Decis. Mak. 26 (6), 565–574. doi:10.1177/0272989x06295361

CrossRef Full Text | Google Scholar

Wang, C., Shi, M., Ji, J., Cai, Q., Zhao, Q., Jiang, J., et al. (2020). Stearoyl-CoA Desaturase 1 (SCD1) Facilitates the Growth and Anti-ferroptosis of Gastric Cancer Cells and Predicts Poor Prognosis of Gastric Cancer. Aging 12 (15), 15374–15391. doi:10.18632/aging.103598

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S.-J., Li, D., Ou, Y., Jiang, L., Chen, Y., Zhao, Y., et al. (2016). Acetylation Is Crucial for P53-Mediated Ferroptosis and Tumor Suppression. Cell Rep. 17 (2), 366–373. doi:10.1016/j.celrep.2016.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Green, M., Choi, J. E., Gijón, M., Kennedy, P. D., Johnson, J. K., et al. (2019). CD8+ T Cells Regulate Tumour Ferroptosis during Cancer Immunotherapy. Nature 569 (7755), 270–274. doi:10.1038/s41586-019-1170-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, J., Zeng, Y., Gao, X., and Liu, T. (2021). A Novel Ferroptosis-Related lncRNA Signature for Prognosis Prediction in Gastric Cancer. BMC Cancer 21 (1), 1221. doi:10.1186/s12885-021-08975-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, R., Ge, Y., Song, W., Ren, J., Kong, C., and Fu, T. (2021). Pyroptosis Patterns Characterized by Distinct Tumor Microenvironment Infiltration Landscapes in Gastric Cancer. Genes (Basel) 12 (10), 1535. doi:10.3390/genes12101535

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, S., Liu, X., Yuan, L., and Wang, F. (2021). A Ferroptosis-Related lncRNAs Signature Predicts Prognosis and Therapeutic Response of Gastric Cancer. Front. Cell Dev. Biol. 9, 736682. doi:10.3389/fcell.2021.736682

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Z., He, M., He, L., Wei, L., and Zhang, Y. (2021). Identification and Validation of a Novel Six-Gene Expression Signature for Predicting Hepatocellular Carcinoma Prognosis. Front. Immunol. 12, 723271. doi:10.3389/fimmu.2021.723271

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Pan, W., Chen, S., Trendel, N., Jiang, S., Xiao, F., et al. (2017). Dynamic Regulation of CD28 Conformation and Signaling by Charged Lipids and Ions. Nat. Struct. Mol. Biol. 24 (12), 1081–1092. doi:10.1038/nsmb.3489

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, T., Shi, S., Zhu, X., Cheang, I., Lu, X., Gao, R., et al. (2022). A Survival Prediction for Acute Heart Failure Patients via Web-Based Dynamic Nomogram with Internal Validation: A Prospective Cohort Study. Jir Vol. 15, 1953–1967. doi:10.2147/jir.s348139

CrossRef Full Text | Google Scholar

Zhang, H., Deng, T., Liu, R., Ning, T., Yang, H., Liu, D., et al. (2020). CAF Secreted miR-522 Suppresses Ferroptosis and Promotes Acquired Chemo-Resistance in Gastric Cancer. Mol. Cancer 19 (1), 43. doi:10.1186/s12943-020-01168-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Ping, L., Du, T., Liang, G., Huang, Y., Li, Z., et al. (2021). A Ferroptosis-Related lncRNAs Signature Predicts Prognosis and Immune Microenvironment for Breast Cancer. Front. Mol. Biosci. 8, 678877. doi:10.3389/fmolb.2021.678877

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Zheng, N., Chen, X., Du, K., Yang, J., and Shen, L. (2022). Establishment and Validation of a Ferroptosis-Related Long Non-coding RNA Signature for Predicting the Prognosis of Stomach Adenocarcinoma. Front. Genet. 13, 818306. doi:10.3389/fgene.2022.818306

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Shi, J., Liu, X., Feng, L., Gong, Z., Koppula, P., et al. (2018). BAP1 Links Metabolic Regulation of Ferroptosis to Tumour Suppression. Nat. Cell Biol. 20 (10), 1181–1192. doi:10.1038/s41556-018-0178-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Peng, Y., He, S., Li, R., Wang, Z., Huang, J., et al. (2021). Apatinib Induced Ferroptosis by Lipid Peroxidation in Gastric Cancer. Gastric Cancer 24 (3), 642–654. doi:10.1007/s10120-021-01159-8

PubMed Abstract | CrossRef Full Text | Google Scholar