Background. Ferroptosis is a special form of cell death with extensive biological associations with various cancers; however, the role of aberrantly expressed ferroptosis-related long non-coding RNAs (lncRNAs) in hepatocellular carcinoma (HCC) remains unclear.
Objectives. To explore the role and prognostic value of ferroptosis-related lncRNAs in HCC and to screen potential therapeutic targets.
Materials and methods. The RNA-seq data for 424 HCC patients and clinical data for 377 HCC patients were obtained from The Cancer Genome Atlas (TCGA) and evaluated using the Pearson’s test to identify differentially expressed lncRNAs. The univariate analysis, least absolute shrinkage and selection operator Cox regression analysis were performed to construct and validate a prognostic risk-score model. The prognostic capacity was evaluated using the Kaplan–Meier method, univariate and multivariate Cox regression, and receiver operating characteristic (ROC) curve analyses. The enrichment analysis was performed to explore the functions of ferroptosis-related lncRNAs from the perspective of tumor immunology.
Results. Seventeen differentially expressed lncRNAs were identified (AL139384.1, AL928654.1, MKLN1-AS, AC145207.8, LINC00205, ZFPM2-AS1, LINC00942, POLH-AS1, AC090772.3, AL603839.3, AC012073.1, AC099850.1, AC026401.3, AP001469.3, AL031985.3, SNHG10, SNHG21) to establish the risk-score model. Survival analyses demonstrated poorer survival in high-risk patients, and the risk score was shown to be a better independent prognostic factor than conventional clinical characteristics. Additionally, we found close correlations between the risk score and immune cell subpopulation functions, as well as between the expression of immune checkpoint and carcinogenic N6-methyladenosine (m6A)-related mRNAs.
Conclusions. The novel ferroptosis-related lncRNA signature may serve as an individualized prognostic prediction tool for patients with HCC.
Key words: immune response, hepatocellular carcinoma, long non-coding RNA, ferroptosis, genes
Hepatocellular carcinoma (HCC) is the main subtype of primary liver cancer and the 4th most common cause of cancer-related deaths worldwide.1 The incidence of HCC in USA increased from 4.4/100,000 in 2000 to 6.7/100,000 in 2012.2 Multidisciplinary treatment strategies have greatly improved the efficacy of HCC treatment, while novel immunotherapy and molecular targeted therapies, such as programmed cell death protein 1 inhibitors, have shown promising results in selected patients.3, 4 However, the survival outcomes of patients with advanced HCC remain poor, with a 5-year overall survival (OS) rate <30%.5, 6 Thus, there is a need to identify novel potential therapeutic targets and accurate prognostic biomarkers at the gene level.
Mounting evidence has indicated the involvement of ferroptosis, a special mode of programmed cell death, in several pathophysiological processes of HCC.7 Ferroptosis is characterized by the over-accumulation of lethal reactive oxygen species (ROS) and production of lipid peroxides,8 and has shown suppressive effects in head and neck, colorectal and ovarian cancer.9, 10, 11 In the context of HCC, ferroptosis induced by the multikinase inhibitor sorafenib and its potential targets has attracted the attention of researchers. Louandre et al. demonstrated that the deactivation of retinoblastoma protein has increased the efficacy of sorafenib in HCC by exacerbating ferroptosis, both in vitro and in vivo.12 Additionally, Sun et al. found that the expression of metallothionein-1G after sorafenib treatment contributed to drug resistance and negatively regulated ferroptosis.13 The increasing appreciation of the importance of ferroptosis in terms of the pathogenesis and treatment of HCC highlights the need to understand the mechanisms underlying ferroptosis-related genetic alterations in HCC.
Long non-coding RNAs (lncRNAs) perform a wide range of regulatory functions at the gene expression level. They have been shown to play a vital role in promoting the progression of HCC through the involvement in proliferative signaling, invasion and metastasis, angiogenesis, as well as immune escape.14 Recent studies have indicated the potential role of lncRNAs in regulating ferroptotic cell death in relation to cancer biology. For example, the overexpression of lncRNA H19 was reported to function as a competing endogenous RNA to enhance the activity of ferroptosis.15 These findings suggest that the ferroptosis–lncRNA relationship may have relevant and important functions in terms of prognostic prediction and treatment in patients with HCC.
The role of ferroptosis-related lncRNAs is currently not fully understood and sequence-based systematic evaluations has been limited. We therefore aimed to explore the role and prognostic value of ferroptosis-related lncRNAs in HCC by constructing a lncRNA signature based on the Cancer Genome Atlas (TCGA), thereby screening potential targets for HCC treatment.
Materials and methods
Transcriptome data for 424 HCC patients (374 tumor and 50 non-cancerous (normally differentiated)) and clinicopathological data for 377 HCC patients (age, sex, tumor grade, TNM stage, follow-up time, and survival outcome) were obtained from the TCGA-HCC database up to October 1, 2021. The clinical characteristics of patients are shown in Table 1. The corresponding ferroptosis-related genes were downloaded from FerrDb database (http://www.zhounan.org/ferrdb/; the first database dedicated to ferroptosis regulators and ferroptosis-disease associations).16 Free data on gene regulators, including drivers, suppressors and markers, could be downloaded for further management and investigation. We identified a total of 382 genes (drivers: 150; suppressors: 109; markers: 123; detailed data given in Supplementary Table 1). The correlations between lncRNAs and HCC were analyzed using the Pearson’s test. We first evaluated the functions of the ferroptosis-related differentially expressed genes (DEGs), and then assessed the related biological pathways with Gene Ontology analysis, and analyzed biological processes (BP), molecular functions (MF), and cellular components (CC) regulated by the DEGs using Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, using R software (ggplot2 package; R Foundation for Statistical Computing, Vienna, Austria).
Construction of ferroptosis-related
lncRNA prognostic signature
We identified the ferroptosis-related lncRNA signature using univariate Cox regression analyses, least absolute shrinkage and selection-penalized Cox regression. The risk scores were calculated as the sum of coefficient lncRNA × expression of lncRNA. The HCC patients were divided into high-risk and low-risk group according to the median risk score, and we compared the OS between the 2 groups using the Kaplan–Meier analysis, and evaluated the prognostic capacity of the risk-score model using receiver operating characteristic (ROC) curves.
Gene set enrichment analysis
and predictive nomogram
We performed gene set enrichment analysis (GSEA) to assess the enriched pathways and to define the lncRNA signature with the KEGG gene set. A nomogram combining the prognostic signatures was created to predict the 1-, 3- and 5-year OS of HCC patients.
Gene expression and immunity analysis
We compared the TMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC, and single-sample GSEA (ssGSEA) algorithms (http://timer.comp-genomics.org/) to estimate cellular immune responses and cellular components between the low-risk and high-risk group according to the risk score. The tumor-infiltrating immune cell subgroups were quantified using ssGSEA and the immunological functions were assessed. A box diagram was used to assess the differential expression of N6-methyladenosine (m6A)-related genes.
Data were analyzed using Bioconductor packages (R software v. 4.1.0; R Foundation for Statistical Computing). Gene expression was compared between normal and HCC tissues using the Wilcoxon test, and correlations between lncRNAs and HCC were analyzed using the Pearson’s test (correlation coefficient R2 > 0.4 at p < 0.01). We identified the differential expression of lncRNAs using the Benjamini–Hochberg method, with the criteria of a false discovery rate <0.05 and log2 fold change ≥ 1. We compared the ssGSEA-normalized HCC DEGs with a genome using the Gene Set Variation Analysis (GSVA) R-package. The sensitivity and specificity of the prognostic signature were assessed using ROC and decision curve analyses,17 and the associations between lncRNAs and clinical characteristics were evaluated with logistic regression analyses and a heatmap graph. The survival of HCC patients based on the ferroptosis-related signature was evaluated using Kaplan–Meier survival analysis with a 2-sided log-rank test. Additionally, immune cell infiltration, immune checkpoint and m6A-related gene expression were compared between the groups using Wilcoxon test. Statistical significance was set at p < 0.05 unless otherwise noted.
of ferroptosis-related genes
We identified 84 ferroptosis-related DEGs, including 71 upregulated and 13 downregulated genes, using KEGG-based analysis. Among these, BP-related genes were involved with organic acid transport and oxidative stress reaction, as well as changes in nutrient levels, toxic substances and extracellular stimuli. The CC-related genes were primarily related to the nicotinamide adenine dinucleotide phosphate hydrate (NADPH) oxidase complex, apical plasma membrane synthesis pathway, lysosomes, and melanin granules, while the MF were connected to the regulation of NADPH production, along with the binding of molecular oxygen and ferric iron (Supplementary Figure 1A). The KEGG-based analysis indicated that the overexpressed genes were commonly involved in ferroptosis, hypoxia-inducible factor-1 (HIF-1) signaling pathway, multiple cancers, synthesis of microRNAs in cancer cells, diabetes, atherosclerosis, viral infection, endocrine system, cell senescence, leukemia, and the mechanistic target of rapamycin signaling pathway (Supplementary Figure 1B).
of ferroptosis-based lncRNAs
A total of 781 ferroptosis-related lncRNAs, including 17 downregulated and 764 upregulated lncRNAs, were found. Sixty-six significant ferroptosis-related lncRNAs that were incorporated into the multivariate Cox analysis were also analyzed using the univariate Cox analysis. Seventeen differentially expressed lncRNAs were finally confirmed as independent prognosis predictors of HCC: AL139384.1, AL928654.1, MKLN1-AS, AC145207.8, LINC00205, ZFPM2-AS1, LINC00942, POLH-AS1, AC090772.3, AL603839.3, AC012073.1, AC099850.1, AC026401.3, AP001469.3, AL031985.3, SNHG10, and SNHG21 (Supplementary Table 2). The network diagram shows the relationship between mRNAs and lncRNAs (Supplementary Figure 1C). On the basis of these findings, we established risk scores and prognostic signatures for the lncRNAs.
Survival results and multivariate analysis
The Kaplan–Meier analysis indicated that the expression of a high-risk lncRNA signature was related to lower survival outcomes (p < 0.001; Figure 1A). Similar results were obtained in a subgroup analysis according to tumor stage (Supplementary Figure 2A). The area under the ROC curve for the lncRNA signature was higher than those for the clinical indices (Figure 1B) in all patients, and in subgroups according to tumor stage (Supplementary Figure 2B). The risk survival status plot (Figure 1C) verified that patients with high-risk scores were prone to low survival outcomes. Additionally, the heatmap demonstrated that most candidate lncRNAs were strongly expressed in the high-risk group (Figure 1C). Taken together, the current risk model demonstrated that specific ferroptosis-associated lncRNAs provided a useful prognostic signature for HCC.
The areas under the curve for 1-, 3- and 5-year OS rates predicted by the novel lncRNA signature were 0.751, 0.779 and 0.801, respectively (Figure 1D). The signature still exhibited good specificity and sensitivity when analyzed in relation to HCC stage (Supplementary Figure 2C). The decision curve analysis confirmed that the lncRNA signature performed better predictive ability than other conventional clinical and pathological characteristics (Figure 1E).
The univariate analysis revealed that lncRNA-based signature (hazard ratio (HR): 1.618, 95% confidence interval (95% CI): 1.589–1.647) together with tumor stage (HR: 1.480, 95% CI: 1.169–1.862) were independent prognostic factors for HCC (Figure 2A, Supplementary Table 3). The multivariate Cox analysis confirmed that the lncRNA signature was an independent prognostic risk factor for HCC (Figure 2B, Supplementary Table 3).
Heatmap analyses illustrated the association between the prognostic signature of ferroptosis-related lncRNA and clinical manifestations (Figure 2C). Predictors in the hybrid nomogram incorporated the risk-score model and clinicopathological characteristics (Figure 2D).
The GSEA revealed that most of the prognostic signatures played roles in metabolism and tumor-related pathways, including the metabolism of amino acids, salts and drugs (such as β-alanine), butanoate, fatty acids, propanoate, pyruvate, and tryptophan, as well as being involved in the citrate cycle, tricarboxylic acid cycle, complement and coagulation cascades, and the peroxisome proliferator-activated receptor signaling pathway (Supplementary Figure 3).
Immunity and gene expression
We evaluated the immune response using the MCPCOUNTER, XCELL, EPIC, TIMER, CIBERSORT, CIBERSORT-ABS, and QUANTISEQ algorithms to determine the association between the ferroptosis-related lncRNA signature and tumor immunity (Figure 3A). The heatmap demonstrated that macrophage and T cell functions were relatively active in the high-risk group. We further investigated the association between immune cell subpopulations and their related functions, and found that the immune cell functions of antigen-presenting cell (APC) co-inhibition, chemokine receptor (CCR), checkpoint, cytolytic activity, inflammation-promoting, T cell co-inhibition/stimulation, and type II interferon response differed significantly between the high-risk and low-risk group (Figure 3B, Supplementary Table 4).
Because immune checkpoint inhibition is an important clinical therapeutic strategy in HCC, we also compared the expression of immune checkpoints between the 2 groups.3, 4 All selected immune checkpoints showed a significantly differential expression between the 2 groups (Figure 4A, Supplementary Table 5). Notably, the expression levels of the tumor necrosis factor superfamily (TNFSF)/TNF receptor superfamily (TNFRSF) signaling pathways were all higher in the high-risk group, as compared with the low-risk group.
The m6A as the cofactor of lncRNAs in regulating post-transcriptional modifications has recently come to the spotlight in lncRNA studies.18 We compared the m6A-related mRNA status between the 2 groups and found that the expression levels of RBM15, YTHDF1/2, ALKBH5, WTAP, METTL3, HNRNPC, and YTHDC1/2 were all significantly higher in the high-risk group (Figure 4B, Supplementary Table 6).
The poor survival and high recurrence rate of HCC have highlighted the need for effective prognostic biomarkers and therapeutic targets. In this study, we constructed a novel ferroptosis-related lncRNA prognostic risk model and explored the roles of the relevant lncRNAs in the immune response and checkpoint regulation in HCC. These findings allowed us to screen possible ferroptosis-related biomarkers and treatment targets, which might inform the development of new therapeutic approaches.
The prognostic model in our study integrated 17 ferroptosis-related lncRNAs. A subset of lncRNAs (MKLN1-AS, LINC00205, ZFPM2-AS1) was reported to be involved in several HCC cellular activities, including cell proliferation, migration and invasion, cell cycle progression, and apoptosis, through microRNAs, signaling pathways and other biological components or proteins (Table 2). Notably, the biological functions of some of the identified lncRNAs (LINC00205, ZFPM2-AS1, LINC00942) have also been associated with other malignancies and might thus be nonspecific to HCC.22, 23, 24, 27, 28, 29, 30, 33, 34 However, numerous other lncRNAs (AL139384.1, AL928654.1, AC145207.8, POLH-AS1, AC090772.3, AL603839.3, AC012073.1, AC099850.1, AC026401.3, AP001469.3, AL031985.3, SNHG10, SNHG21) lack corresponding basic studies to validate their prognostic roles in the development of HCC. Further studies are therefore needed to establish the associations of these identified lncRNAs with ferroptosis and their contribution to the regulation of HCC pathogenesis through ferroptosis.
Based on the risk-score method, the ferroptosis-related lncRNA signature classified HCC patients into high-risk and low-risk groups, with significantly different OS rates. Both univariate and multivariate Cox regression analyses identified the risk score as an independent prognostic factor. Additionally, ROC curve and nomogram analyses showed that the risk score had higher sensitivity and better clinical applicability than traditional standards for predicting HCC survival. Furthermore, stratification analyses by stage confirmed the prognostic predictive value of the risk score signature for each tumor stage.
Numerous studies have established the role of the lncRNA prognostic signature for HCC using bioinformatics analysis. The biological processes related to the lncRNAs in the signatures in these studies included pyroptosis,35 hypoxia,36 metabolism,37 and epithelial–mesenchymal transition.38 Chen et al. recently reported a ferroptosis-related signature model that integrated both mRNAs and lncRNAs,39 with outstanding HCC predictive efficiency. Notably, lncRNA ZFPM2-AS1 was also included in this signature, indicating its strong connection with the development and progression of HCC. In the current study, we focused solely on ferroptosis-related lncRNAs and attempted to explore their role from the perspective of the immune response, immune checkpoints and m6A status.
Mounting evidence has revealed that crosstalk between lncRNAs and immune cells modulates the tumor immune response. However, both immune response and enrichment analyses revealed no significant differences in activities of any immune cell species between the high-risk and low-risk groups. Moreover, levels of all examined immune cell subpopulations were relatively lower in the high-risk group, suggesting a possible regulatory role for lncRNAs in escape from tumor immunity.
Notably, the expression levels of all the immune checkpoints examined in this study were significantly higher in the high-risk group. Among these, TNFSF/TNFRSF signaling pathways showed potential as immunotherapy targets in combination with ferroptosis induction. Interactions among TNFSF/TNFRSF mediate signaling that controls immune cell survival, proliferation and differentiation.40 As the most important superfamily member, the function of tumor necrosis factor alpha (TNF-α) in promoting tumor growth and its association with a poor prognosis of HCC have been widely demonstrated, and TNF-α inhibition has accordingly demonstrated appreciable anti-tumor effects in various studies.41, 42, 43, 44 Therefore, the relationship between ferroptosis and TNF-α signaling, as well as the combined efficacy of TNF-α inhibition and ferroptosis in HCC, warrant further investigation.
The importance of reversible chemical modifications of RNA, particularly methylation, has recently gained a renewed interest. The m6A is the most prevalent form of internal mRNA modification and has been implicated in HCC carcinogenesis.45 Accordingly, the current study showed that several m6A-related genes were highly expressed in the high-risk group (RBM15, YTHDF1/2, ALKBH5, WTAP, METTL3, HNRNPC, and YTHDC1/2). The METTL3 was reported to process oncogenic functions in HCC by promoting m6A modification of the tumor suppressor gene SOCS2 in a YTHDF2-dependent pathway.46 At the same time, METTL3 was shown to work collaboratively with YTHDF1 to activate the protein translation of Snail, leading to increased epithelial–mesenchymal transition and metastasis of HCC.47 The investigation of the biological links between m6A and lncRNAs suggests that the modification of m6A might affect the structure, stability, expression, and subcellular distribution of lncRNAs, thereby promoting tumor growth.48, 49, 50 From this perspective, we hope that these findings might identify feasible targets for HCC intervention.
This study had some limitations. It was based mainly on integrative bioinformatics, and histological examination of corresponding tissue samples, follow-up clinical data and experimental validation of the findings are currently lacking. Additionally, the biological functions of the signature components might be nonspecific to HCC, given that they have also been associated with other malignancies. Therefore, further validation of the clinical applicability and clarification of the exact role of ferroptosis in HCC based on our risk score are needed.
In summary, we constructed a ferroptosis-related lncRNA signature to predict the prognosis of HCC. These results might provide new indications for understanding the mechanisms of ferroptosis-related lncRNAs in regulating the immune response, and for the development of individualized treatments.
The supplementary tables and figures are available at https://doi.org/10.5281/zenodo.6492096. The package contains the following files:
Supplementary Figure 1. Gene Ontology (GO), KEGG analysis of DEGs, and the relationship between identified lncRNAs and mRNAs expression.
Supplementary Figure 2. Kaplan–Meier analysis and the area under the curve (AUC) predicting value results of HCC patients in TNM stage I, II and III.
Supplementary Figure 3. GSEA for ferroptosis-related lncRNAs based on TCGA.
Supplementary Table 1. Information of driver, marker and suppressor.
Supplementary Table 2. Statistical information on univariate and multivariate Cox analysis of identified lncRNAs.
Supplementary Table 3. Statistical information on univariate and multivariate Cox analysis of clinical characteristics and the risk score.
Supplementary Table 4. Statistical information on immune function between the high-risk and low-risk group.
Supplementary Table 5. Statistical information on immune checkpoint between the high-risk and low-risk group.
Supplementary Table 6. Statistical information on m6A-related gene expression between the high-risk and low-risk group.