Identification and Validation of Ferroptosis-Related LncRNAs Signature as a Novel Prognostic Model for Chronic Lymphocytic Leukemia

Introduction

Chronic lymphocytic leukemia (CLL) is a B-cell hematological malignancy, which is characterized by the clonal expansion and the accumulation of CD5+/CD23+ mature B lymphocytes in the peripheral blood, bone marrow, and secondary lymphoid tissues.¹ Since the 1960s, more and more studies have revealed the complex pathogenesis and explored the treatment of CLL.² With great advances in the biological understanding of CLL and the emergence of the small molecular targeted agents, such as BTK inhibitors, PI3K inhibitors and BCL-2 inhibitors, the treatment of CLL has shifted from the era of traditional chemoimmunotherapy (CIT) to chemotherapy-free era.³ However, chemo-free therapy is not enough to achieve deep remission, high-risk patients still suffer from treatment resistance and encounter poor prognosis. Therefore, it is crucial to identify novel molecules to predict the prognosis and achieve individualized treatment in CLL.

Ferroptosis is an iron and lipid peroxidation-dependent programmed cell death, which is distinct from apoptosis, necrosis, autophagy, and proptosis.^4,5 Since the term was defined by Dr Stockwell in 2012, numerous studies have shown that ferroptosis is involved in many important pathological processes, including cancer development and drug resistance.^6,7 Currently, different degrees of ferroptosis sensitivity have been observed in cancer cells from different types of lymphoma and leukemia, suggesting the unique role of ferroptosis in hematological malignancies. For example, dimethyl fumarate induces ferroptosis by lipid peroxidation particularly in GCB DLBCL.⁸ Artesunate may regulate autophagy and ferroptosis via impairing the STAT3 signaling pathway in DLBCL cells as well as enhance ferroptosis in Burkitt’s lymphoma by activating ATF4-CHOP-CHAC1 pathway.^9,10 TXNRD1 is recognized as a key regulator involved in the ferroptosis of CML cells induced by cysteine depletion in vitro.¹¹ APR-246 induces cell death by ferroptosis in AML and triggers ferritinophagy and ferroptosis of diffuse large B-cell lymphoma cells with distinct TP53 mutations.^12,13 Additionally, the prognostic significance of ferroptosis-related genes is increasingly recognized in hematological malignancies.^14–16

Long noncoding RNAs (lncRNAs) are defined as RNAs longer than 200 nucleotides that are not translated into functional proteins. Studies have shown that lncRNAs are widely expressed and involved in gene regulation. Depending on their localization and specific interactions with DNA, RNA, and proteins, lncRNAs play an important role in cell activities, such as chromatin function regulation, assembly and function of nucleosome, stability and translation of cytoplasmic mRNA, and so on.¹⁷ Recently, the relationship between lncRNAs and ferroptosis sheds light on cancer research. Ferroptosis-related lncRNA signatures were established to predict the prognosis of cancer. All evidence suggests ferroptosis-related lncRNAs (FRLs) could be potential biomarkers and provide a rationale for clinical treatment in CLL.

Therefore, we obtained gene expression data of CLL from public databases to screen out differentially expressed prognostic FRLs. In this study, we first established a novel prognostic risk model with 6 FLRs (PRKCQ, TRG.AS1, LNC00467, LNC01096, PCAT6 and SBF2.AS1) and further validated its reliability and sensitivity in CLL. Moreover, our results revealed the difference in functional enrichment analysis and immune cell infiltration between the two clusters. Notably, we found similar results from the validation cohort from our center. Overall, the ferroptosis-related lncRNAs risk score prognostic model may accurately and effectively predict the prognosis of patients with CLL, thus providing new insights into clinical management.

Materials and Methods

Data Collection

Firstly, we obtained gene sequencing profiles and related clinical information of 151 CLL patients in the GSE22762 dataset from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) as a training cohort. According to Table S1, two hundred fifty-nine ferroptosis genes including markers, drivers and suppressors were comprehensively acquired from FerrDb (http://www.zhounan.org/ferrdb/). Then, the lima package in R was used to analyze the correlation between the expression of ferroptosis-related genes and lncRNAs in CLL patients. According to the limitation of R (Relevant Coefficient) >0.4 and P<0.001, 250 FRLs were found. Through univariate Cox regression, 69 lncRNAs with prognostic value (P<0.05) were screened out. Subsequently, another dataset GSE50006 containing 188 CLL patients and 32 normal samples was used for the identification of 291 differentially expressed genes (DEGs). Finally, 20 prognostic ferroptosis-related differently expressed lncRNAs (PFDELs) were considered for further study. All data from GEO and FerrDb are available.

Generation and Evaluation of the Ferroptosis-Related Prognostic Score Model

To avoid the potential risk of overfitting, the least absolute contraction and selection operator (LASSO) regression analysis was performed with the “glmnet” R package.¹⁸ Based on the expression of 20 FRLs in CLL, the following equation was used to calculate the risk score of individuals in the cohort: FPS=βlncRNA1 × Expression lncRNA1 + βlncRNA2 × Expression lncRNA2 + · · · + βlncRNAn × Expression lncRNAn (the value of β represents the corresponding coefficients). The samples in the training cohort were classified into high- and low-risk groups based on the suitable cutoff value. The “pheatmap” R package was applied to generate the heatmap of the 6 FRLs expression, while the “ggbiplot” R package was used to generate the principal component analysis (PCA) and t-distributed stochastic neighbor embedding (tSNE) for dimensionality analysis. Kaplan–Meier curve analysis and long-rank tests were used to compare the difference of the overall survival (OS) time, treatment-free survival (TFS) time and progression-free survival (PFS) time between the high- and low-risk groups with the “survminer” R package. Finally, the R package “ROC survival” was performed to describe ROCs for the prediction of the OS, TFS and PFS at different times.

Functional Enrichment Analysis and Immune Cell Infiltration

Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) enrichment analyses of DEGs were conducted by “cluster Pofolier” R package to explore the biological process and potential signaling pathways of two risk groups in our model. Subsequently, the CIBERSORT (https://cibersort.stanford.edu/) algorithm was applied to assess the population of 22 distinct immune cells in the tumor microenvironment (TME). We continued to assess the immune infiltration score of 28 immune cell subtypes and 13 immune pathways with single-sample gene set enrichment analysis (ssGSEA) to compare the immune features between the high- and low-risk groups.¹⁹

Ethics Statement

This study was approved by the Ethics Committee of the Jiangsu Province Hospital and written informed consent was obtained from all patients according to the Declaration of Helsinki. All methods were carried out in accordance with relevant guidelines.

Sample Collection in Validation Cohort and RNA-Seq Analysis

We collected 36 peripheral blood samples of CLL in the First Affiliated Hospital of Nanjing Medical University (Jiangsu Province Hospital). Peripheral blood mononuclear cells (PBMC) from 36 CLL samples were first collected, and a CD19⁺ B cell selection kit (Miltenyi Biotech, Gladbach, Germany) was used to extract B-CLL cells. Total RNA was extracted and assessed with RNeasy Micro Kit (QIAGEN, Hilden, Germany) to ensure quality. The prepared sequencing libraries were sequenced by a HiSeq X Ten high-throughput sequencing system. The sequences were mapped to hg38 (human genome 38) and aligned with Bowtie and BLAT (the BLAST-like alignment tool).

Statistical Analysis

All statistical data were analyzed by R software (version 4.11) and IBM SPSS software (version 21.0). Pearson’s correlation analysis was used to analyze the correlation between the lncRNAs and ferroptosis-related genes. The t-test or Kruskal–Wallis test was conducted for continuous variables to determine significant factors. Survival analysis was estimated with the Kaplan–Meier methods and the differences between the high- and low-risk groups were evaluated by long-rank tests. Univariate and multivariate regression analyses were performed to identify independent prognostic factors. P<0.05 was defined as statistically significant: * P < 0.05; ** P < 0.01; *** P < 0.001.

Results

Flowchart of the Study

The flowchart provided a clear summary of the construction of the prognostic model and further analysis (Figure 1).

Figure 1 Flow chart of the design in the study.

Construction of the Novel Ferroptosis-Related lncRNA Risk Score Prognostic Model for CLL Patients

Firstly, to construct the prognostic model, twenty differentially expressed FRLswere obtained after intersecting 69 lncRNAs with prognostic value, 291 DEGs and 250 FRLs, which were significantly associated with the OS of CLL patients through univariate COX regression analysis (P<0.05) (Figure 2A). The forest plot showing the prognostic value of the 20 lncRNAs is presented the Figure 2B. We revealed the gene expression correlation by correlation analysis of 20 FRLs (Figure 2C). Then, LASSO Cox regression analysis was conducted to establish a new ferroptosis-related lncRNA prognostic score (FPS) model (Figure 2D‒E), which was calculated as follows: FPS=0.4720* LNC00467 expression+0.3359* LNC01096 expression+0.1326* PCAT6 expression+0.0831* SBF2.AS1 expression‒0.6319* PRKCQ.AS1 expression‒0.1653* TRG.AS1 expression. The risk score of PRKCQ.AS1 and TRG.AS1 had negative coefficients, indicating the protective effect on prognosis, whereas LNC00467, LNC01096, PCAT6 and SBF2.AS1 were observed as risk factors with positive coefficients (Figure 2F). We also described the relationship between the lncRNAs, ferroptosis-related genes and the risk types (Figure 2G). The Kaplan–Meier analysis showed that patients with the high level of protective factors (PRKCQ.AS1 and TRG.AS1) had favorable OS (P < 0.0001 and P < 0.0001, respectively) and TTFT (P < 0.0001 and P = 0.0097, respectively), while patients with the high level of risk factors (LNC00467, LNC01096, PCAT6 and SBF2.AS1) had worse outcomes (Figure S1A‒B). Furthermore, the Sankey diagram shows the relationship between the lncRNAs, ferroptosis-related genes and the risk types.

Figure 2 Construction of the Novel Ferroptosis-Related lncRNA Risk Score Prognostic Model for CLL patients (A)Venn diagram of the generation of 20 PFDELs. (B)Forest plots showed the HR and p value of the 20 PFDELs. (C) Correlation analysis of the 20 FRGs. (D and E) The LASSO coefficient profiles of the 6 FRLs signature. (F) Coefficients (λ) of the 6 FRLs in the FPS model. (G) Sankey diagram of the relationship between the lncRNAs, ferroptosis-related genes and the risk types.

Evaluation of the Novel Ferroptosis-Related LncRNA Risk Score Prognostic Model in the GEO Dataset

The samples were divided into high and low-risk groups based on the optimal risk score of 1.08 as the cutoff (Figure 3A). As shown in Figure 3B, the patients in the high-risk group tended to have worse prognosis compared to those in low-risk group. In addition, the heatmap showed lower expression levels of PRKCQ and TRG.AS1, and the higher expression levels of LNC00467, LNC01096, PCAT6 and SBF2.AS1 in the high-risk group, whereas the expression levels were opposite in the low-risk group (Figure 3C). Then, PCA and t-SNE indicated that CLL patients were well distributed into two groups (Figure 3D‒E). Moreover, Kaplan–Meier curve analysis indicated that patients in the high-risk group had shorter OS and TTFT than those in the low-risk group (P < 0.0001, respectively) (Figure 3F‒G). The ROC curves and the corresponding AUC values showed superior performance in predicting prognosis. The AUC values of the ROC curves for predicting the 1-, 3- and 5-year OS were 0.829, 0.824 and 0.821, respectively, and the 1-, 3- and 5-year TFS were 0.788, 0.857 and 0.839, respectively (Figure 3H and I). The data in the training cohort suggested that the FPS model could provide an accurate prediction of the outcomes of CLL patients.

Figure 3 Evaluation of the Novel Ferroptosis-Related LncRNA Risk Score Prognostic Model in the GEO Dataset (A) The distribution and optimal cut-off value of risk scores in the training cohort. (B) The distributions of OS status, OS and risk score. (C) Heatmap of the expression levels of the 6 selected FRLs. (D and E) t-SNE and PCA plot of the CLL cohort visualizing the distribution of the low and high-risk groups. (F and G) Kaplan–Meier survival curves for OS and TFS of CLL patients stratified by FPS risk score. (H and I) Time-dependent ROC curves of the risk model for predicting the 1-, 3- and 5-year OS and TFS.

Functional Enrichment Analyses and Immune Cell Infiltration Between the High and Low-Risk Groups in CLL

As shown in Figure 4A, the heatmap displayed the top 40 of 226 DEGs between the FPS-compliant low- and high-risk groups. To explore the role of the FPS model on the occurrence and progression of CLL patients, Go and KEGG enrichment analyses were conducted. Go analysis indicated that the DEGs were associated with biological processes (BP) such as T cell activation, lymphocyte differentiation and T cell differentiation (Figure 4B). Meanwhile, KEGG analysis suggested that the chemokine signaling pathway, cytokine-cytokine receptor interaction, hematopoietic cell lineage, T helper (Th) cell differentiation, TCR pathway and NF-κB signaling pathway were significantly enriched (Figure 4C). Considering the enrichment of immune cells in the DEGs, we continued to investigate the relationship between the risk scores and immune infiltration. CIBERSORT algorithm was implemented to assess the infiltration level of 22 distinct immune cells. Different infiltrations were found between high- and low-risk groups, including naive CD4⁺ T cells, resting natural killer (NK) cells, monocytes, naive B cells, plasm cells, activated mast cells, follicular helper T cells, regulatory T cells and activated NK cells (Figure 4D, all P<0.05). Based on this, we proceeded to use ssGSEA algorithm to estimate differences in the enrichment scores of immune cells and immune-related functional pathways in high- and low-risk groups. The scores of activated B cells, activated CD4⁺ T cells, activated CD8⁺ T cells, CD56 bright/dim natural killer cells, central memory CD8⁺ T cells, effector memory CD8⁺ T cells, gamma delta T cells, mast cells, memory B cells, monocytes, myeloid-derived suppressor cell, natural killer cells, natural killer T cells, plasmacytoid dendritic cells, regulatory T cells, follicular helper T cells, type 1 T helper cells, type 2 T helper cells and type 17 T helper cells were significantly different in the two groups (Figure 4E, all P<0.05). Furthermore, APC co-suppression/co-stimulation, chemokine receptor (CCR), checkpoints, cytokine activity, inflammation promotion, T-cell co-suppression or co-stimulation and type II interferon (IFN) response were significantly decreased in the high-risk group, while HLA pathway was increased (Figure 4F, all P<0.05). In summary, the evidence of functional analyses and immune cell infiltration could indicate that the FPS model may be correlated with the immune microenvironment in CLL.

Figure 4 Functional Enrichment Analyses and Immune Cell Infiltration Between the High and Low-risk groups in CLL(A) The heatmap displaying the differentially expressed genes between the low-risk group and high-risk group. (B and C) The representative results of GO and KEGG analysis of DEGs between two groups. (D) Boxplots of the relative fraction of 22 immune cell types between two groups. (E and F) The score of immune cell types and immune-related functions using ssGSEA analysis, respectively.

Validation of the Novel Ferroptosis-Related LncRNA Risk Score Prognostic Model

To further demonstrate the feasibility of the FPS model, thirty-six patients with CLL diagnosed between 2012 and 2017 in our center were included in this study as a validation cohort. We calculated the FPS based on the expression of 6 FLRs, and then the patients were also separated into high- and low-risk groups based on the optimal cutoff. It was consistent with the training cohort from the GEO database (Figure 5A‒C). Besides, Kaplan–Meier curves were also applied to evaluate the prognostic role of the 6 FRLs by OS, FPS and TTFT (Figure S2). Our model showed that the higher scores were associated with poorer survival outcomes. We observed significant differences in OS, PFS and TTFT between the high- and low-risk groups (OS: P=0.0066; FPS: P=0.016; TTFT: P=0.038, Figure 5D‒F). Moreover, the AUC values of the ROC curves showed high accuracy to predict the 3- and 5-year OS (OS: AUC = 0.805 and 0.851) and the 1-, 3- and 5-year PFS (PFS: AUC = 0.744, 0.618 and 0.720, respectively), and TTFT as well (TTFT: AUC = 0.660, 0.593 and 0.754, respectively) (Figure 5G‒I).

Figure 5 Validation of the Novel Ferroptosis-Related LncRNA Risk Score Prognostic Model (A) The distribution and the value of risk scores in the validation cohort. (B) The distributions of OS status, OS and risk score in the validation cohort. (C) Heatmap of the expression levels of the 6 FRLs in the validation cohort. (D-F) Kaplan–Meier survival curves for OS, PFS and TFS of CLL patients stratified by FPS risk score. (G-I) Time-dependent ROC curves of the risk model for predicting the 1-, 3- and 5-year OS, PFS and TFS.

Independent Prognostic and Clinical Characteristics of the 6-FRLs Prognostic Model

Next, univariate analysis Cox regression analyses were applied to explore the relationship between the FPS model and clinical characteristics, such as stage, B symptoms, gender, age, lymphocytes, hemoglobin, platelets, LDH, β2-MG, TP53 disruption, IGHV status and CLL International Prognosis Index (CLL-IPI). The results indicated that FPS high risk (P=0.013), Binet stage B or C (P=0.023), B symptoms (P=0.030), age >65years (P=0.034), hemoglobin <100g/L (P=0.008), platelets (P=0.002), β2-MG >3.5 mg/L (P=0.006), TP53 disruption (P=0.016), IGHV unmutated (P=0.001), and CLL International Prognosis Index (CLL-IPI) (P<0.001) were significantly related to poor OS (Figure 6A). Finally, FPS (P=0.034) and CLL-IPI (P=0.030) were selected as independent prognostic factors for OS prediction by multivariate Cox regression analysis (Figure 6B). Notably, a prognostic nomogram containing the FPS and CLL-IPI score was plotted to quantitatively predict the 2-, 3- and 5-year OS (Figure 6C). Furthermore, calibration curves of 2-, 3- and 5-year OS were described to confirm the consistency between the predicted survival rate and the actual survival rate (Figure 6D-F). Overall, the FPS model displayed benefited performance in predicting the survival of CLL.

Figure 6 Independent Prognostic and Clinical Characteristics of the 6-FRLs Prognostic Model (A and B) The univariate and multivariate Cox regression model analyses for assessment of the prognostic value of clinical characteristics and FPS. (C) The nomogram model for predicting 2-, 3- and 5-year OS rate of CLL patients. (D-F) The calibration plot analysis to assess the nomogram accuracy for OS prediction at 2-, 3- and 5-year.

Discussion

CLL is a highly heterogeneous disease. In the past 50 years, relevant advances have been made in the understanding of genetics and molecular biology, clinical characteristics, prognostic stratification and treatment of CLL.²⁰ To date, different prognostic models have been developed for CLL prognostic stratification,^21,22 but CLL-IPI is widely used with the most convincing and the highest acceptance.²³ However, with the emergence of the novel oral inhibitors, the CLL-IPI might be partly insufficient.²⁴ Therefore, it is essential to develop a novel prognostic model to better predict outcomes.

Since the 21st century, various novel types of unprogrammed cell death including apoptosis, necrosis, autophagy, proptosis, ferroptosis and cuproptosis have been gradually recognized to reveal the complexity of biology.^4,25–27 Among them, a unique modality of cell death, ferroptosis is induced by iron-dependent lipid peroxidation and regulated by multiple cellular metabolic events, including redox homeostasis, mitochondrial activity, and metabolism.⁶ Although lncRNAs have no ability to encode proteins, it can participate in competing endogenous RNAs, interacting with proteins or inhibiting translation to regulate ferroptosis.^16,28–31 However, the potential value of ferroptosis-related lncRNAs in CLL is still unclear. Hence, based on the data from GEO database, we developed a novel ferroptosis-related lncRNAs prognostic model.

In the prognostic model, PRKCQ.AS1 and TRG.AS1 are considered as protective factors, whereas LINC00467, LINC01096, PCAT6 and SBF2.AS1 were observed as risk factors. Actually, lncRNA PRKCQ.AS1 plays an important role in the process of lymphoma including DLBCL and HL.^32,33 Consistent with our study, it was also shown to be a protective factor in DLBCL.³² Another study showed that up-regulation of TRG.AS1 could protect against glucocorticoid-induced osteoporosis in a rat model by regulating miR-802-mediated CAB39/AMPK/SIRT-1/NF-κB axis.³⁴ To date, no studies of TRG.AS1 have been reported in hematological malignancies except in this study. In previous studies, LINC00467, LINC01096, PCAT6 and SBF2.AS1 have been reported to be overexpressed and play oncogenic roles in both solid cancers and hematological malignancies. Through competing endogenous RNA, they promoted tumorigenesis, tumor progression and tumor metastasis by regulating various signaling pathways, such as Wnt/β-catenin, NK-κB and STAT3 pathway. LINC00467 was reported by Bo to play a role in signaling pathways of peroxisomal lipid metabolism and immunity.³⁵ PCAT6 was included in a novel ferroptosis-related lncRNA model to accurately predict the prognosis and immune microenvironment of hepatocellular carcinoma.³⁶ SBF2.AS1 was identified and validated to be a ferroptosis-related gene signature for predicting survival in glioma.³⁷ The biological role of LINC01096 remains poorly understood. Although all lncRNAs in our model were first reported in CLL, they unexpectedly assisted us in establishing an accurate and convincing prognostic model.

Notably, the biological function of the DEGs between the FPS-compliant low- and high-risk groups was analyzed in our study. By KEGG analysis, we found several pathways were significantly enriched, including the chemokine signaling pathway, TCR pathway and NF-κB signaling pathway. Besides, T helper (Th) cell differentiation also plays a unique role. Previously, some studies have reported that ferroptosis and tumor-infiltrating immune cells could influence each other.^38,39 Chen indicated that follicular helper T cells (Tfh) cells are susceptible to undergoing ferroptosis with the persistent stimulation of TCR, which results in higher ROS production in Tfh cells,⁴⁰ which may explain the enrichment of T cell activation and differentiation between the DEGs. Besides, the NF-κB pathway has been reported to be aberrantly hyperactivated in chronic lymphocytic CLL development and evolution.⁴¹ Considering the TME plays an essential role in the development, growth, and survival of the malignant B-cell clone in CLL, the infiltrating abundance of tumor immune cells also differed between high- and low-risk groups by measuring the proportions of immune-infiltrating cells with CIBERSORT. Compared to the high-risk group, resting NK cells, monocytes, resting mast cells and eosinophils, which were recognized to play anti-tumor roles, were infiltrated at higher levels in the low-risk group. However, Tregs and Tfhs, considered as immunosuppressive roles, were infiltrated at higher levels in high-risk populations.^42–44 To some degree, the immunosuppressive microenvironment could explain the fact that patients with high-risk scores undergo worse prognosis than those with low-risk scores in our study. All the genetic complexity and diverse and complex tissue TME interactions participate in the heterogeneous course of CLL. Although the introduction of small molecule inhibitors of BTK, PI3K and BCL-2 has greatly improved the prognosis of CLL patients, prolonging the survival time, limitations still exist on achieving a deep and permanent remission. Combination strategies of small molecule inhibitors and immune checkpoint inhibitors or Immunomodulatory drugs, targeting the TME, may provide light on patients with high-risk.

In conclusion, we identify ferroptosis-related lncRNAs and construct a novel and accurate prognostic model. The reliability and sensitivity of the model were verified both in the external database and in our data. We sincerely hope that our findings would shed insights for further study and clinical practice.

Indubitably, our study still has some limitations. First, the 6-FRLs signature prognostic model may be not easy to be widely used in clinical practice because of the price of RNA-sequencing, but it could be accurate and useful in clinical trials to some degree. Second, the model was only validated with 36 patients in our center, and more external validations are needed to assess the reliability and sensitivity. Finally, the molecular mechanism of the lncRNAs included in the prognostic model has not been fully recognized. In future, functional analysis and mechanistic studies of the 6 ferroptosis-related lncRNAs are expected to be carried out to explore the specific value in CLL progression.