Ferroptosis-Related Gene Signature: A New Method for Personalized Risk Assessment in Patients with Diffuse Large B-Cell Lymphoma

Introduction

Diffuse large B-cell lymphoma (DLBCL) is a genetically heterogeneous disease that is divided into germinal center B cell and activated B cell subtypes in transcription.¹ The first-line treatment for DLBCL is chemoimmunotherapy using rituximab, cyclophosphamide, doxorubicin, vincristine, and prednisone (R-CHOP).² The 5-year survival rate of patients who undergo standard care can reach 60%–70%.³ Approximately 30–40% of patients respond poorly to the standard treatment and mortality is highest within the first 2 years after diagnosis.⁴ Despite a substantial increase in the investigation of molecular genetics and novel targeted agents in recent years, attempts to improve outcomes by combining standard therapy with novel targeted agents have thus far yielded disappointing results.⁵ One of the major barriers to the effective use of novel therapies targeting specific pathways is the biological heterogeneity of DLBCL.

In the era of rituximab, the most common international prognostic indicators used in DLBCL patients include the International Prognostic Index (IPI), the Revised IPI, and the National Comprehensive Cancer Network (NCCN-IPI).⁶ All of these systems use standard clinical and laboratory factors to predict outcomes without explicitly considering the molecular heterogeneity of DLBCL.^5,6 Hence, it is becoming increasingly important to define high-risk patients in the early stages of disease and to enter them into clinical trials investigating alternative therapies.

Rapidly developing ferroptosis research suggests the possibility for its use in the formulation of prognoses in cancer patients. Ferroptosis was first proposed by Dixon et al in carcinogen studies.⁷ It refers to the regulated form of iron-dependent cell death caused by the accumulation of lipid-based reactive oxygen species.⁸ Studies indicate that a network of long non-coding RNA and competitive endogenous RNA plays an important role in tumorigenesis and ferroptosis in lung cancer.⁹ The growth of tumor cells can evidently be inhibited by promoting ferroptosis and ferroptosis is also related to tumor resistance.¹⁰ Cancer-associated fibroblasts reportedly inhibit the ferroptosis of gastric cancer cells by targeting ALOX15 and blocking the accumulation of lipid-reactive oxygen species(ROS), and secrete exosomal miR-522 thereby promoting drug resistance.¹¹ In a sensitivity analysis conducted using 177 cancer cell lines, lymphoma cell lines were particularly susceptible to GPX4-regulated ferroptosis.¹² Erasing an inducer of ferroptosis induces hypertrophy and slows tumor growth in mouse lymphoma models.¹³

Some basic research suggests that the induction of ferroptosis is a potential therapeutic approach for lymphoma, but there is a lack of research on relationships between clinical lymphoma data and ferroptosis. In the current study the significance of ferroptosis-related mRNA expression was investigated with respect to prognostic evaluation in DLBCL patients. mRNA expression profiles and corresponding clinical data derived from DLBCL lymphoma patients were downloaded from public databases. An analysis of differentially expressed genes related to ferroptosis in the Gene Expression Omnibus cohort was then conducted, and it was subsequently verified in other databases. Lastly, functional enrichment analysis was conducted to investigate potential mechanisms involved. A prognostic model based on a large number of ferroptosis-related genes was thus developed. We believe that this robust prognostic model will improve risk stratification in DLBCL patients and facilitate more accurate assessment for their clinical management.

Materials and Methods

Data Collection

Clinical data and gene expression profiling of patients with DLBCL were obtained from the National Center for Biotechnology Information Gene Expression Omnibus (NCBI Gene Expression Omnibus). Data series were downloaded in a normalized expression matrix file format and used directly in analyses. A training cohort containing 470 patients from the GSE31312 dataset was generated,¹⁴ as was a validation cohort containing 414 patients from the GSE10846 dataset.¹⁵

Identification of Ferroptosis-Related Genes Associated with Overall Survival

Univariate Cox regression analysis was used to investigate relationships between gene expression and prognostic values, and the results were strengthened by retaining only genes with p values < 0.05. A total of 259 ferroptosis-related genes were retrieved from FerrDb (http://www.zhounan.org/ferrdb). Lastly, ferroptosis-related genes that were significantly associated with overall survival (OS) in the datasets were extracted.

Ferroptosis-Related Risk Model

The GSE31312 dataset was used to construct a ferroptosis-related model. Using least absolute shrinkage and selection operator (LASSO) Cox regression analysis, a risk formula was generated to predict prognoses. The R package “glmnet” statistics software was used to perform the LASSO regression. An individualized risk score was then calculated for each patient, and the median risk score was used as the cut-off to divide all patients into a high-risk group or a low risk-group. Kaplan-Meier survival analysis and Log rank testing were used to evaluate the difference in OS between the high-risk group and the low-risk group. Time-dependent receiver operating characteristic (ROC) curves were constructed to assess prognostic value,¹⁶ and relationships between risk scores and clinical characteristics were assessed. Differences in clinical factors between the high-risk and low-risk groups were evaluated via stratified analysis. The analysis was performed using the Sangerbox tool (http://sangerbox.com/Tool).

Univariate and Multivariate COX Regression

SPSS 26 was used for univariate and multivariate COX regression analyses of clinical information and the risk prediction model.

Nomogram Construction and Validation

Based on the GSE31312 cohort, a nomogram integrating ferroptosis-related gene features and IPI components was established to assess the probability of individualized OS at 1, 3, and 5 years using the R package (https://cran.r-project.org/web/packages/rms/). The discriminant capacity of the nomogram was graphically evaluated via calibration mapping.

Functional Enrichment Analyses and Protein-Protein Interaction Network Generation

Gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the “clusterProfiler” R package, then the results were displayed in R language via the ggplot2 package. Adjusted p values < 0.05 were considered statistically significant. A protein–protein interaction (PPI) network was obtained from the STRING online database (https://string-db.org/) and visualized via Cytoscape (version 3.7.2).

Results

Identification of Ferroptosis-Related Genes Associated with OS

The clinical characteristics of all patients in the GSE31312 and GSE10846 datasets are shown in Table 1. Univariate Cox regression was used to assess associations between ferroptosis-related genes and OS. A total of 104 ferroptosis-related genes that were significantly associated with OS in the training dataset were extracted (Supplementary Table S1). A LASSO Cox method was used to select the most useful predictors among the 104 genes in the GSE31312 dataset, and 8 prognostic ferroptosis-related genes with the highest predictive values based on OS were identified (Figure 1A and B).

Table 1 The Clinical Characteristics of the DLBCL Patients from GEO

Figure 1 LASSO Cox regression (A) performed using 104 ferruginous disease-related genes to obtain the prognostic characteristics of 8 genes (B). A total of 104 ferroptosis-related genes was used to construct the PPI network (C).

Enrichment Analysis and PPI Networks

To gain insight into the mechanisms ferroptosis-related genes are involved in DLBCL, 104 ferroptosis-related genes were used to perform gene ontology enrichment and KEGG pathway analyses. As expected, in gene ontology enrichment analysis 104 genes were related to oxidative stress-related gene ontology terms, such as “cellular response to oxidative stress” and “response to oxidative stress” (Supplementary Figure S1A). KEGG pathway analysis also indicated that these genes are primarily involved in the ferroptosis pathway (Supplementary Figure S1B). A total of 104 ferroptosis-related genes were used to establish a PPI network (Figure 1C).

Construction of a Prognostic Gene Signature Based on Ferroptosis-Related Gene

LASSO Cox regression was used to determine the optimal weight coefficients for the eight prognostic ferroptosis-related genes, then a prognostic risk score model was established with the following formula:

Risk score = expression level of HIF1A*(−0.05455544) + expression level of NFE2L2*(−1.87952976) + expression level of FH*(−0.09480619) + expression level of LAMP2*(−0.34575386) + expression level of ZEB1*(−0.43828031) + expression level of NGB*(0.16447965) + expression level of PSAT1*(0.47117607) + expression level of CXCL2*(−0.20840707)

Risk scores were calculated for all 470 patients in the training set, who were then assigned to a high-risk group or a low-risk group based on whether their risk score was higher or lower than the median risk score. The distributions of risk scores, death status, survival time, and expression levels of the eight genes used to construct the prognostic signature were visualized to evaluate prognostic differences between the two groups (Figure 2A–C). The expression of 8 genes was different between high and low risk groups (Supplementary Figure S2). The respective areas under time-dependent ROC curves of the eight-gene model at 1, 3, and 5 years were 0.71, 0.66, and 0.66 (Figure 2D), indicating that the model had high sensitivity and high specificity with regard to predicting the prognosis of DLBCL patients. In Kaplan-Meier curve analysis, patients with high-risk scores had significantly poorer OS than patients with low-risk scores (hazard ratio 8.62; 95% confidence interval 5.39–13.77; p < 0.0001) (Figure 2E). In order to explore the correlation among the eight core genes, we constructed a PPI network, and the results showed that six genes out of eight were closely related (Supplementary Figure S1C).

Figure 2 Distributions of risk scores (A), death status, survival time (B), and expression levels of the eight genes used to derive the prognostic signature tool (C) were visualized to evaluate the prognostic difference between a high-risk group and a low-risk group after LASSO Cox regression. Time-dependent ROC curve and area under the curve of the signature (D). Kaplan-Meier plots of overall survival in the high-risk and low-risk groups in the training set determined via Log rank testing (E).

Relationships Between Risk Score and Clinical Features

Correlations between risk score and age, Eastern Cooperative Oncology Group (ECOG), IPI score, lactate dehydrogenase (LDH), stage, and extranodal sites were calculated. Risk score was positively correlated with age (Figure 3A), IPI score (Figure 3C), LDH (Figure 3D), DLBCL stage (Figure 3E), and extranodal sites (Figure 3F), but not with ECOG (Figure 3B).

Figure 3 Correlations between risk score and age (A), ECOG (B), IPI score (C), LDH (D), DLBCL stage (E), extranodal sites (F).*p < 0.05, ***p < 0.001, ****p < 0.0001, ns, p > 0.05.

Subgroup Analysis of the Eight-Gene Signature in the Training Dataset

The prognostic capacity of the eight-gene signature was further validated in DLBCL subgroups. The signature could predict survival outcomes in clinical subgroups defined by parameters such as activated B cell or germinal center B cell subtype, clinical stage, and IPI score. Patients in the high‐risk group exhibited poorer prognoses irrespective of activated B cell or germinal center B cell subtype (Figure 4A and D), clinical stage (Figure 4B and E), and IPI score (Figure 4C and F). Thus, the signature exhibited independent predictive capacity in clinical application.

Figure 4 Kaplan-Meier plots of overall survival in high-risk and low-risk activated B cell subtype (A), germinal center B cell subtype (D), stage I–II (B), stage III–IV (E), IPI < 2 (C), and IPI ≥ 2 (F) derived via Log rank testing. Red is the high-risk group and blue is the low-risk group.

Validation of the Eight-Gene Signature in Patients in the GSE10846 Dataset

In the GSE10846 dataset patients classified as low-risk based on the median risk score had longer OS than patients classified as high-risk (hazard ratio 2.1; 95% confidence interval 1.50–2.93; p < 0.0001) (Figure 5A). Respective areas under the ROC curve for predicting OS at 1, 3, and 5 years were 0.610, 0.654, and 0.647, respectively (Figure 5B).

Figure 5 Kaplan-Meier plots of overall survival in high-risk and low-risk subgroups in the GSE10846 dataset derived via Log rank testing (A). The time-dependent ROC curve and AUC in the validation set (B).

Univariate and Multivariate COX Regression

In analyses performed solely on datasets with IPI components (GSE31312 and GSE10846), the eight-gene signature was an independent prognostic indicator (Table 2).

Table 2 Univariable and Multivariable Cox Regression Analysis in DLBCL

Construction and Validation of the Nomogram

IPI components and the eight-gene model were used to create nomograms to predict 1, 3, and 5-year survival (Figure 6A). A calibration chart indicated acceptable agreement between the predicted survival rate and the actual survival rate, suggesting that the nomogram may reliably predict DLBCL patient prognoses in clinical practice (Figure 6B–D).

Figure 6 Construction of a nomogram to forecast 12, 36, and 60-month survival using the IPI components and the eight-gene model (A). The calibration chart shows that the predicted survival rate was consistent with the actual survival rate for 12 months (B), 36 months (C), and 60 months (D). IPI, International Prognostic Index.

Discussion

DLBCL is a distinct disease, but exhibits extensive molecular heterogeneity and different clinical manifestations, morphological characteristics, immunophenotypes, molecular subsets, and clinical prognoses.¹⁷ DLBCL patient risk stratification and prognosis remain a challenge for clinicians. In the rituximab era it was found that IPI was somewhat deficient with respect to distinguishing risk groups, especially high-risk groups, and its prognostic value was found to be lacking.⁶ The study formulated the Revised-IPI, which divides patients treated via the R-CHOP regimen into three groups. This overcame the problem that the survival curves of patients in low-risk, low-medium-risk, medium-high-risk, and high-risk groups merged.⁶ Notably however, there was still an inability to adequately identify high-risk patients. Studies indicate that NCCN-IPI classification can quantify the degree of patient risk to an extent. However, all three scoring systems fail to identify high risk group with long-term OS clearly below 50%.⁶ Moreover, the progression of tumors is affected by numerous biological parameters in addition to those ascertained via clinical and laboratory investigations.

Ferroptosis may have a profound impact on carcinogenesis and prognosis. It is an iron-dependent non-apoptotic form of cell death that can inhibit tumor growth and increase the sensitivity of chemotherapy.¹⁸ Multiple small molecular compounds have been identified that induce ferroptosis by targeting iron metabolism and lipid peroxidation,¹⁹ but the role and significance of mRNA expression of ferroptosis-related genes in the prognostic assessment of DLBCL are currently unclear. In the current study a prognostic signature based on eight markers that strongly and reliably predicted individualized prognoses in DLBCL patients on the basis of ferroptosis-related genes was established. In conjunction with the incorporation of an IPI component a comprehensive prognosis nomogram was developed that exhibited a good capacity to accurately predict prognoses and positive benefit of DLBCL.

The eight gene products included in the signature developed in the current study were ZEB1, PSAT1, NGB, NFE2L2, LAMP2, HIF1A, FH, and CXCL2. ZEB1 is the main component of the transcription factor network controlling epithelial-to-mesenchymal transition, which is associated with tumor growth and metastasis of malignant tumors via effects on critical cellular functions and states including differentiation, proliferation, responses to DNA damage, and cell survival.²⁰ c-Myc is a common transcription factor of epithelial–mesenchymal transition (EMT) including Snail, ZEB1, and ZEB2 via direct associations with their promoter regions containing E-box motifs, suggesting close correspondence between c-Myc and ZEB1 in the regulation of EMT.²¹ In DLBCL, ZEB1 activates SNHG14 and PD-L1, which in turn promotes immunoevasion by certain tumors.²² Recent research suggests that ZEB1 overexpression increases susceptibility to ferroptosis, whereas ZEB1 silencing decreases susceptibility to it.²³ Many studies indicate that NRF2 plays an important role in the therapeutic response to ferroptosis-targeted therapies in hepatocellular carcinoma.^24,25 In signal pathway, NRF2/CBS signal transduction confers elastin-induced ferroptosis resistance.²⁶ The NRF2-Keap1 pathway is essential for cancer cell growth.²⁷ Interaction between ARF and NRF2 is critical for p53-independent ferroptosis, and inhibition of the NRF2-ARE pathway reversed resistance to ferroptosis in head and neck cancer.²⁸ NRF2 is reportedly an effective therapeutic target for chemotherapeutic drugs to make cancer sensitive.²⁵ The hypoxia-inducible factor (HIF) pathway was the driving factor for the vulnerability of ferroptosis.²⁹ That study demonstrated that HIF1A can act as a tumor promoter by degrading p53 protein via binding to five response elements in the p53 promoter, and increases their invasive and metastatic activity.³⁰ HIF1A also activates autophagy, a lysosomal degradation pathway that may promote DLBCL cell survival.³¹ PSAT1 is upregulated in Burkitt lymphoma via a MYC-dependent ATF4 transcription factor mechanism, and overexpression of PSAT1 is reportedly associated with poor prognoses in patients with Burkitt lymphoma,³² nasopharyngeal carcinoma,³³ colorectal carcinoma,^34,35 and epithelial ovarian cancer.³⁶ Studies suggest that NGB may be a tumor suppressor that inhibits cell growth, cell aggregation, and tumorigenicity.³⁷ LAMP2 deficiency increases the risk of ferroptosis in retinal pigment epithelial cells.³⁸ ZEB1 and HIF1A can reportedly influence ferroptosis in DLBCL, but the molecular functions of other ferroptosis-related genes in DLBCL patients remain to be elucidated.

The current study had some limitations. Due to the unavailability of control group samples in the Gene Expression Omnibus databases and the retrospective study design, there was a degree of bias. Thus, a well-designed, international, prospective, multicenter clinical trial is needed to validate our findings in the future. In addition, further functional studies are warranted to explore the molecular functions of the identified ferroptosis-related genes during DLBCL progression.

Conclusion

This is the first time a ferroptosis-related prognostic model has been established in DLBCL patients that is capable of serving as a promising prognostic indicator in combination with clinical IPI components. The model may also facilitate the development of new clinical ferroptosis-targeted therapies in patients with DLBCL.