Construction and Validation of a Ferroptosis-Related lncRNA Signature as a Novel Biomarker for Prognosis, Immunotherapy and Targeted Therapy in Hepatocellular Carcinoma

Recently, immunotherapy combined with targeted therapy has significantly prolonged the survival time and improved the quality of life of patients with hepatocellular carcinoma (HCC). However, HCC treatment remains challenging due to the high heterogeneity of this malignancy. Sorafenib, the first-line drug for the treatment of HCC, can inhibit the progression of HCC by inducing ferroptosis. Ferroptosis is associated with the formation of an immunosuppressive microenvironment in tumours. Moreover, long non-coding RNAs (lncRNAs) are strongly associated with ferroptosis and the progression of HCC. Discovery of ferroptosis-related lncRNAs (FR-lncRNAs) is critical for predicting prognosis and the effectiveness of immunotherapy and targeted therapies to improve the quality and duration of survival of HCC patients. Herein, all cases from The Cancer Genome Atlas (TCGA) database were divided into training and testing groups at a 6:4 ratio to construct and validate the lncRNA signatures. Least Absolute Shrinkage and Selection Operator (LASSO) regression and Cox regression analyses were used to screen the six FR-lncRNAs (including MKLN1-AS, LINC01224, LNCSRLR, LINC01063, PRRT3-AS1, and POLH-AS1). Kaplan–Meier (K–M) and receiver operating characteristic (ROC) curve analyses demonstrated the optimal predictive prognostic ability of the signature. Furthermore, a nomogram indicated favourable discrimination and consistency. For further validation, we used real-time quantitative polymerase chain reaction (qRT-PCR) to analyse the expression of LNCSRLR, LINC01063, PRRT3-AS1, and POLH-AS1 in HCC tissues. Moreover, we determined the ability of the signature to predict the effects of immunotherapy and targeted therapy in patients with HCC. Gene set enrichment analysis (GSEA) and somatic mutation analysis showed that ferroptosis-related pathways, immune-related pathways, and TP53 mutations may be strongly associated with the overall survival (OS) outcomes of HCC patients. Overall, our study suggests that a new risk model of six FR-lncRNAs has a significant prognostic value for HCC and that it could contribute to precise and individualised HCC treatment.


INTRODUCTION
Liver cancer ranks sixth in the world in terms of incidence and is the fourth most common cause of cancer-related deaths (Villanueva, 2019). HCC is the most common malignancy of the liver, accounting for approximately 90% of liver cancer cases (Llovet et al., 2021). In the early stages of HCC, curative resection, transplantation, and local ablation are the preferred treatment options, but most patients with HCC are diagnosed at an advanced stage and cannot be cured using the above approaches (Llovet et al., 2021). Despite the many advances in the treatment of HCC, such as immunotherapy and targeted therapies, the median survival time for patients with advanced HCC is only 11-13 months (Villanueva, 2019). Given the high mortality rate of HCC, it is imperative to further explore biomarkers that can predict the effectiveness of HCC immunotherapy or targeted therapy, discover new prognosisrelated biomarkers, and identify new therapeutic targets to prolong the survival time of patients with HCC.
Ferroptosis is a type of regulated cell death which, unlike apoptosis, pyroptosis, and necroptosis, is caused by the accumulation of lipid peroxide (Stockwell et al., 2017). Recent studies have shown that ferroptosis is associated with the progression and treatment response of various types of tumours (Chen et al., 2021a). For instance, ferroptosis can induce inflammation-related immunosuppression in the tumour immune microenvironment (TIME), thereby promoting tumour progression, whereas mutations in TP53 and RAS genes are closely related to ferroptosis (Chen et al., 2021a). However, the regulatory mechanisms of ferroptosis in HCC remain to be investigated and are far from being applied for the treatment of HCC (Nie et al., 2018;Capelletti et al., 2020). Therefore, identifying the main regulators associated with ferroptosis is a crucial way to broaden the therapeutic approach to HCC.
Long non-coding RNAs (LncRNAs) are most commonly defined as RNAs that are more than 200 nucleotides in length and do not encode proteins (Quinn and Chang, 2016). In HCC, lncRNAs are specifically involved in protein synthesis, degradation, and epigenetic regulation (Mai et al., 2019). HCC-related lncRNAs can regulate the protein modifications of transcription factors and influence protein function by regulating the corresponding cellular signalling pathways (Li et al., 2017;Yan et al., 2017;Zhang et al., 2018). Furthermore, recent studies have demonstrated that GABPB1 and its antisense lncRNA GABPB1-AS1 are closely linked to erastininduced ferroptosis and that GABPB1 and GABPB1-AS1 can be used as treatment targets for HCC patients (Qi et al., 2019). However, the mechanism underlying the role of FR-lncRNAs in HCC is unclear, and its importance in the therapeutic and prognostic value of HCC needs to be further elucidated.
In this study, we identified and validated a prognostic signature based on ferroptosis-related (FR)-lncRNAs. We also explored the association between the risk model and immune cell infiltration, immune function, sensitivity of immune checkpoint inhibitor (ICI) treatment, and sorafenib treatment. Furthermore, we used gene set enrichment analysis (GSEA) to reveal the pathways enriched in the high-and low-risk groups of the signature. Finally, we validated the significantly differentially expressed lncRNAs in this signature between the tumour and para-tumour tissues in vitro. In conclusion, this new prognostic signature is more accurate and convenient than previous signatures in predicting the prognosis and treatment outcomes of HCC patients.

Patient Datasets and Processing
Gene expression profiles and clinical data from patients with HCC were downloaded from The Cancer Genome Atlas (TCGA) database. Considering the likelihood of non-cancer deaths, patients with a survival time of ≤30 days and missing expression data were excluded (n = 35), leaving 342 HCC patients in the final cohort. The flow chart of the analysis process is presented in Figure 1. A total of 342 HCC cases were randomly assigned to the training and testing sets at a ratio of 6:4 for systematic analysis employing the R project "caret" package. Both the training and testing sets needed to comply with the following requirements: 1) cases were randomly classified into training and test groups, and 2) the clinical characteristics of the subjects in both groups were similar. Gene transfer format (GTF) files were downloaded from Ensembl for annotation to distinguish messenger RNAs (mRNAs) from lncRNAs for subsequent analysis. The data of ferroptosis genes (FRGs) were downloaded from FerrDb (Zhou and Bao, 2020) and led to the identification of 259 FRGs. (Supplementary Table  S1). Correlation analysis was performed between the FRGs and all the lncRNAs. Those with ferroptosis gene correlation coefficients higher than 0.4 and p values lesser than 0.001 were considered FR-lncRNAs. To identify the differentially expressed FR-lncRNAs, we used the R package "limma" for differential expression analysis of FR-lncRNAs. Significantly differentially expressed FR-lncRNAs fulfilled the following criteria: p < 0.05 and |log2FC| ≥ 1.

Construction of a Prognostic Risk Signature
We screened lncRNAs associated with OS in patients with HCC using univariate Cox regression analysis. p < 0.05 was considered statistically significant. For the training group, LASSO regression analysis was performed using the R project "glmnet" package to further select the screened lncRNAs. To prevent over-fitting, 1,000 rounds of cross-validation were used to adjust the parameter selection, and the partial likelihood deviation met the minimum criteria. We then performed a multivariate Cox regression analysis and identified six FR-lncRNAs. Subsequently, we calculated the six FR-lncRNAs corresponding coefficients to construct a prognostic risk score profile for HCC. The formula was established as follow: Risk score coef lncRNA1 × lncRNA1 expression + coef lncRNA2 × lncRNA2 expression + · · · · + coef lncRNA N × lncRNA N expression. dependent ROC curves were constructed, and the area under the time-dependent ROC curve (AUC) values for 1, 3, and 5 years were calculated. Moreover, multivariate Cox regression analysis was conducted to assess whether the risk score model could be used as an independent predictor of OS in HCC patients. The predictive efficacy of this six-FR-lncRNA signature was further confirmed in the testing and whole cohorts.

Nomogram Construction and Assessment
The R package "rms" and "regplot" were used to construct the nomogram based on the six-FR-lncRNA signature, which integrated the signature, age, and stage. We have estimated the prognosis of HCC patients at 1, 3, and 5 years in the nomogram. Finally, calibration curves were plotted to assess the accuracy and reliability of the nomogram. frozen in liquid nitrogen (−196°C) immediately after tissue excision. Total cellular RNA was extracted using the TRIzol reagent (Invitrogen, California, United States). Complementary DNA (cDNA) was synthesised using the PrimeScript Reverse Transcriptase Reagent Kit (Takara Bio, Inc., Japan). Amplification and detection were performed using TB Green ® Premix Ex Taq ™ II (Takara, Tokyo, Japan) in the ABI Step One Plus Real-Time PCR system (Applied Biosystems). β-Actin was used as an endogenous control. The expression level of PRRT3-AS1, LNCSRLR, LINC01063, and POLH-AS1 was normalised to that of βactin using the 2 −ΔΔCt method. The primer sequences are listed in Supplementary Table S2. This study was approved by the ethics committee of the Chinese PLA General Hospital (approval no. S2018-111-01). Written informed consent was obtained from all patients.

Estimation of the Immune Cell Types and Analysis of Tumour Immune Infiltration Cell Types
Quantitative translation of the tumour tissue transcriptome data into the absolute abundance of immune and stromal cells by Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis was conducted to assess the proportion of 22 human immune cell subpopulations (Newman et al., 2015;Becht et al., 2016). Differences in each type of immune cell were compared between the high-risk and low-risk groups to assess differences in the tumour immune microenvironment between the two groups. Then, we implemented a single sample GSEA (ssGSEA) approach using the gene set variation analysis (GSVA) and "GSEABase" packages to analyse the immune functions and inflammatory infiltration profiles (Newman et al., 2019). Immune infiltrating cells and functions were compared between the high-and low-risk groups using the Wilcoxon test.

Analysis of the Expression of Immune Checkpoint Genes and Sensitivity to Clinical Treatment in the High-and Low-Risk Groups
We analysed the expression of 47 immune checkpoint genes between the high-and low-risk groups. Then, we used the tumour immune dysfunction and exclusion (TIDE) algorithm to predict the effectiveness of immunotherapy (Jiang et al., 2018). We also calculated the half-inhibitory concentration (IC50) of sorafenib in the entire dataset. The difference in the IC50 of sorafenib between the two groups was compared using Wilcoxon signedrank test and the results were obtained by using the R packages "pRRophetic" and "ggplot2."

Somatic Variant Analysis
Gene somatic mutation data based on the whole-exome sequencing platform of TCGA-Liver Hepatocellular Carcinoma (LIHC) datasets were downloaded from the Genomic Data Commons (GDC) database on August 22, 2021. The downloaded Mutation Annotation Format (MAF) files of simple nucleotide variation (workflow type: varScan2 variant aggregation and masking) were analysed using the R package "maftools."

Construction of the lncRNA-FRG Co-Expression Network, Copy Number Variations of FRGs Associated With the Signature, and Functional Analysis
According to the results of the previous co-expression analysis, Cytoscape (version 3.8.2) was used to visualise the co-expression network of prognostic FR-lncRNAs and FRGs. CNV data were obtained through the UCSC Xena program, which is publicly available under specific guidelines. The landscape of genomic CNV in chromosomes was plotted using the "RCircos" R package. GSEA software, version 4.1.0 was used to reveal the signal transduction pathways (Subramanian et al., 2005). The gene sets used in this work were c2.cp.kegg.v7.4.symbols.gmt (Kanehisa et al., 2017) and h.all.v7.4.symbols.gmt (Subramanian et al., 2005). A nominal p-value of <0.05 was used as the screening criterion.

Statistical Analysis
We used R version 4.1.0 (Institute for Statistics and Mathematics, Vienna, Austria), GraphPad Prism 8 (GraphPad Software Inc., La Jolla, CA, United States), and SPSS 25.0 (SPSS, IL, United States) to analyse our data. Survival analysis was performed using K-M and log-rank tests. The t-test or Mann-Whitney U test was used to compare two independent groups. Categorical data were analysed using the χ 2 test. TIDE scores between the high-risk group and low-risk group were compared using the Wilcoxon test. p < 0.05 was considered statistically significant (*p < 0.05, **p < 0.01, ***p < 0.001).

Screening of FR-lncRNAs and Differential Expression Analysis
The data of 374 HCC samples and 50 normal samples were downloaded from the LIHC project of TCGA. Then, data annotation was performed according to the GTF file from Ensembl and divided into lncRNA and mRNA data, followed by FRG and lncRNA co-expression analysis and FR-lncRNA coexpression network visualization using the Sankey diagram in Figure 2A. Finally, 548 FR-lncRNAs were identified totally (|r | > 0.4, p < 0.001), of which 336 were identified as differentially expressed FR-lncRNAs according to the standard p < 0.05, and | log2FC| > 1 ( Figure 2B). Among these differentially expressed FR-lncRNAs, 332 were upregulated and four were downregulated.

Construction and Confirmation of the FR-lncRNAs Signature
A total of 342 HCC patients were included according to the inclusion criteria. We performed univariate Cox regression analysis and screened 115 differentially expressed FR-lncRNAs associated with OS (p < 0.05) (Supplementary Table S3). The 342 HCC patients (entire dataset) were assigned to the training and testing groups at a ratio of 6:4 randomly, and the clinical features were similar between the three groups (Table 1). Then, LASSO regression and Cox proportional hazard model analyses were applied to identify the best model in the training dataset ( Figures 2C,D). Furthermore, multivariate analysis identified six FR-lncRNAs risk scores for HCC (Supplementary Table S4), according to the following formula: Figure 2E shows the relationship between each lncRNA and OS. The expression of six FR-lncRNAs was significantly upregulated in HCC by comparing HCC and normal tissues in the TCGA transcriptome profiles ( Figure 2B).
We then calculated risk values for each case and categorised all cases as low and high risk based on the median threshold of the training dataset. As shown in the heat map, the expression of the six lncRNAs was higher in the high-risk group of HCC patients than in the low-risk group, gradually increasing patients' mortality as risk score increased ( Figure 3A). In the training dataset, K-M analysis showed that the survival rate of the low-risk group was significantly higher than that of the high-risk group (p = 9.603e-06) ( Figure 3B). Then, we constructed an ROC curve and the results showed that the AUC values of 1, 3, and 5 years OS were 0.812, 0.758, and 0.709, respectively ( Figure 3C). Furthermore, univariate and multivariate Cox regression analyses were performed on age, gender, tumour stage, histological grade, and risk score to verify whether the FR-related risk score could be used independently as an indicator of survival in OS. The results of univariate Cox regression analysis showed that tumour stage (p < 0.001) and risk score (p < 0.001) were associated with prognosis ( Figure 3D). Multivariate Cox regression analysis illustrated that tumour stage (p < 0.001) or risk score (p < 0.001) could be used as an independent predictor of prognosis for HCC patients.

Validation of the FR-lncRNA Signature
We used the test dataset and entire dataset to further validate the accuracy and reliability of the prognostic signature. In the test dataset, the distributions of expression profiles of the six-lncRNA, risk score, and OS status were consistent with those of the training dataset ( Figure 4A). Similarly, K-M analysis of the testing set ( Figure 4B) and entire data set ( Figure 4D) showed that patients in the low-risk group had a higher survival rate than those in the high-risk group. The AUC values of the 1, 3, and 5 years OS in the test dataset were 0.845, 0.787, and 0.700, respectively ( Figure 4C) and the AUC values of the 1, 3, and 5 years OS in the entire dataset were 0.826, 0.768, and 0.714, respectively ( Figure 4E). Subsequently, we compared the accuracy of our signature with that of published articles according to the values of the AUC and concordance index (C-index) (Supplementary Figure S1). The efficacy of our signature is better than that of other studies using an identical protocol. Consistent with the results of the training group, therefore, the FR-lncRNA risk score could be used as an independent predictor of OS based on univariate and multivariate Cox regression analyses of the test dataset and the entire dataset ( Figures 4F-I).
Furthermore, to determine whether the risk score model could be used as the best predictor of survival, we included age, sex, tumour stage, and pathological grade as candidate predictive indicators. We analysed the AUC curves for 1-year prognosis in the three datasets and found that our signature possessed the highest AUC values among these factors ( Figures 5A-C).
Significantly higher expression of MKLN1-AS and LINC01224 in HCC tissues compared to normal liver tissues has been demonstrated in previous studies (Dan Gao et al., 2020). We used qRT-PCR to detect the expression of the other four FR-lncRNAs in HCC tissues and adjacent normal tissues. The results showed that the expression levels of LNCSRLR ( Figure 5D), LINC01063 ( Figure 5E), PRRT3-AS1 ( Figure 5F), and POLH-AS1 ( Figure 5G) were higher in HCC tissues than in normal tissues.
A Nomogram Combining the Risk Score, Age, and Tumour Stage to Predict Survival Time in Patients With HCC Our prognostic nomogram integrated the six-lncRNA signature, age, and tumour stage. Then, we used the prognostic nomogram to predict the prognosis of HCC patients at one, three, and 5 years after diagnosis.
( Figure 6A). Moreover, calibration plots of 1, 3 and 5 years survival probabilities also showed good agreement between nomogram predictions and actual observations ( Figure 6B). Similar results have been observed for the testing set ( Figure 6C) and the entire dataset ( Figure 6D).

Differences in Immune Cell Infiltration, Expression of Immune Checkpoint Genes and Sensitivity to Clinical Treatments Between High-and Low-Risk Groups
The entire dataset was used to explore the relationship between the signature and immune cell infiltration, the expression of immune checkpoint genes, and clinical treatments. The differences in the infiltration of 22 immune cell types in the tumour tissue of all HCC patients are shown in Figure 7A. Patients with HCC in the high-risk group had higher ratios of M0 macrophages, follicular helper T cells, memory B cells, memory activated CD4 T cells, and regulatory T cells (Tregs) than those in the low-risk group (p < 0.05); meanwhile, patients with HCC in the low-risk group had higher ratios of activated NK cells, M1 macrophages, monocytes, naïve B cells, resting mast cells, and memory resting CD4 T cells than those in the high-risk group (p < 0.05) ( Figure 7B). To explore the relationship between the expression of the six lnc-RNAs and immune infiltration level in HCC, lollipop plots were constructed showing Spearman's correlation coefficient and statistical significance in Supplementary Figure S2. Notably, the six lnc-RNAs were positively correlated with the populations of M0 macrophages and Tregs and negatively correlated with the populations of monocytes and M2 macrophages (p < 0.05). We also applied the ssGSEA method to the RNA sequencing data of HCC samples to assess immune cell infiltration and related functions. The populations of immune cells, including activated dendritic cells (aDCs), B cells, immature dendritic cells (iDCs), macrophages, mast cells, neutrophils, Natural Killer (NK) cells, plasmacytoid dendritic cells (pDCs), and Tregs, were found to be markedly different between the two groups ( Figure 7C). Moreover, immune signature comparisons revealed that low-risk patients had higher cytolytic activity, type-I IFN response, and type-II IFN response than high-risk patients, whereas the opposite was observed for MHC class I ( Figure 7D). Given the importance of ICIs in the treatment of HCC, we further analysed the differential expression of immune checkpoint genes between the high-and low-risk groups. We found that patients in the high-risk group had higher expression of immune checkpoint genes, such as PDCD-1 (PD-1), CTLA4, LAG3, HAVCR2 (TIM3), and TIGIT ( Figure 8A), compared to that in the low-risk group. We then used TIDE to assess the potential clinical efficacy of immunotherapy in the two groups. A higher TIDE prediction score represents a higher potential for immune evasion, suggesting that patients are less likely to benefit from ICI treatment. Our results indicated that the high-risk group had a lower TIDE score than that of the low-risk group, implying that patients in the high-risk group could benefit more from ICI therapy than patients in the low-risk group ( Figure 8B). In addition to ICI treatment, we attempted to identify the association between the signature and the efficacy of sorafenib for the treatment of HCC. We found that the low-risk group was associated with lower IC50 for sorafenib (p < 0.01), which suggested that the signature could be used as a potential predictor of sorafenib treatment sensitivity ( Figure 8C).

Gene Mutation Analysis
We processed simple nucleotide variation data using the "maftools" package in R. A waterfall plot displayed the top 20 mutated genes in patients with HCC ( Figures 9A,B). A higher proportion of somatic mutations (TP53) were found in the high-risk group (FDR <0.01, p < 0.001) after a Fisher's exact test was used to analyse the mutation differences between two groups. Figures 9C,D summarise the mutation information for the high-and low-risk gene groups, respectively.

Construction of Co-Expression Network, CNV Analysis of FR Genes Associated With the Signature, and GSEA
To explore the interaction between the six FR-lncRNAs and gene expression in HCC, Cytoscape was used to visualise the lncRNA and mRNA co-expression network ( Figure 10A). The Sankey diagram showed the relationship between FR-lncRNAs, FRGs, and OS in patients with HCC ( Figure 10B). It was found that the FRGs associated with signature had different levels of CNV events. RB1 harboured the most CNV events. Meanwhile, YY1AP1, RPL8, HSF1, MAFG, PGD, and STMN1 also had relatively high CNV events ( Figure 10C). Figure 10D shows the location of CNV events in the FRGs related to the signature on the chromosome. Then, we performed GSEA to explore the biological effects of the six FR-lncRNAs signature, and the entire group dataset was used for GSEA analysis. The results revealed that the high-risk group showed significant enrichment in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to cancer processes, such as base excision repair, cell cycle, endocytosis, mismatch repair, nucleotide excision repair, and the WNT signal transduction pathway. Correspondingly, ferroptosis and metabolism-related pathways, such as betaalanine metabolism, drug metabolism cytochrome P450, fatty acid metabolism, and retinol metabolism were significantly enriched in the low-risk group ( Figure 10C). Furthermore, GSEA of the hallmark gene sets indicated that PI3K-AKT-MTOR, TGF-β, NOTCH, P53, and WNT-β-Catenin pathways were enriched in the high-risk group significantly, whereas fatty acid metabolism, bile acid metabolism, and xenobiotic metabolism pathways were highly enriched in the low-risk group ( Figure 10D).

DISCUSSION
The prognosis for patients with HCC is poor, mainly because a large proportion of patients are diagnosed with HCC at an advanced stage (Llovet et al., 2021). Recently, despite the great success achieved using a combination of anti-PD-L1 with anti-VEGF therapies in advanced HCC, there are still a majority of patients who either do not respond to these treatments or do not have a lasting clinical benefit because of high tumour heterogeneity or treatment resistance (Finn et al., 2020;Pinter et al., 2020). Therefore, to maximize the benefits to patients and improve the effectiveness of systematic therapy, it is necessary to explore reliable molecular biomarkers to predict the effectiveness of immunotherapy, targeted therapies, and the prognosis of HCC patients.
LncRNAs play a key role in chromatin structure, cell growth, gene expression, differentiation, and development, and mutations or dysregulation of their expression are associated with a variety of diseases, particularly malignancies (Esteller, 2011;Bhan et al., 2017). Recent studies have revealed that lncRNAs are associated with the following six cancer hallmarks: proliferation, growth suppression, motility, immortality, angiogenesis, and viability (Schmitt and Chang, 2016). For the treatment of HCC, lncRNAs can be used as biomarkers to predict the efficacy of surgery, radiotherapy, chemotherapy, and immunotherapy, and are expected to be a potential tool for individualised HCC diagnosis and treatment (Yuan et al., 2021). Moreover, the TP53 mutation is involved in the expression of specific lnc-RNAs and gains oncogenic function by creating a complex network of interacting pathways (Di Agostino, 2020). For example, lincRNA-p21 serves as a transcriptional repressor in the p53 pathway and was the first lncRNA identified as being transcriptionally induced by wild-type p53 (Huarte et al., 2010). In 26 resected pancreatic cancer specimens, lncRNA1611 was significantly highly expressed in the cancerous pancreatic tissue of 22 patients compared to normal pancreatic tissue and was positively correlated with TP53 mutation .
However, few studies have documented the relationship between TP53 mutation and lncRNA expression (Lin et al., 2019). Furthermore, several recent studies have found that lncRNAs are strongly associated with ferroptosis in cancer (Xie and Guo, 2021;Yao et al., 2021;Zhang et al., 2021;Zuli Wang et al., 2021). However, to date, the number of FR-associated lncRNAs identified in HCC is scarce, and research on FR-associated lncRNAs in HCC is limited. RNA-targeted therapies are developing rapidly, such as the expression of HOTAIR-sbid, a mutant of lncRNA HOTAIR, which has been shown to reduce cell motility, invasiveness, and response to TGF β-induced epithelial-mesenchymal transition Battistelli et al., 2021). In the era of precision medicine, by exploring FR-related lncRNAs, we aimed to obtain a risk model that could predict the efficacy of clinical treatments and patient prognosis. We successfully constructed and validated a new biomarker comprising six FRrelated lncRNAs in patients with HCC based on TCGA dataset. Statistical analysis showed that this risk model has good robustness and predictive power and that it can independently predict OS in patients with HCC. Furthermore, the predictive accuracy of our signature surpassed that of the previously reported prediction signature of HCC dependent on FR-lncRNAs or FRGs (Liang et al., 2020;Chen et al., 2021;Liang Wang et al., 2021;Liang et al., 2021;Nie et al., 2021;Wan et al., 2021;Xu et al., 2021). The nomogram consisting of age, stage, and six FR-related lncRNA risk scores could be used to visually predict the OS of individual patients with HCC at 1, 3, and 5 years. According to the calibration plots, the nomogram has good prediction accuracy and the predicted results match the actual results well. Sorafenib and immunotherapy can inhibit the progression of tumours, including HCC, by inducing ferroptosis (Louandre et al., 2013;Nie et al., 2018;Friedmann Angeli et al., 2019;Wang et al., 2019). Ferroptosis can promote exposure to tumour antigens, thereby enhancing the immunogenicity of the TIME and efficacy of immunotherapy . Meanwhile, ICIs act primarily through the activation of an anti-tumour immune response driven by cytotoxic T cells, and this anti-tumour effect can induce ferroptosis in cancer cells (Chen et al., 2021a). Moreover, CD36 expressed by CD8 + T cells leads to an accumulation of lipid peroxides in CD8 + T cells through the uptake of fatty acids in the tumour environment, which in turn leads to an elevated iron ion content, increased iron death processes, and reduced secretion of cytotoxic cytokines (Ma et al., 2021). Therefore, HCC patients with different ferroptosis characteristics and immunophenotypes may respond differently to sorafenib or immunotherapy. In our study, the proportion of Tregs, M0 macrophages, follicular helper T cells, and memory B cell subsets in the CIBERSORT analysis showed significant infiltration in the high-risk group. The results of the ssGSEA analysis showed that the proportion of aDCs, iDCs, macrophages, and Tregs was significantly increased in the TIME of the high-risk group. In addition, immune checkpoint-associated genes were more highly expressed in patients of the high-risk group compared to the low-risk group, which may provide a basis for identifying patients who may respond to ICI therapy. In contrast, the infiltration of activated NK cells and M1 macrophages, which are immune-promoting cell subsets, was significantly increased in the low-risk group in the CIBERSORT analysis. Similarly, the results of the ssGSEA analysis showed that the low-risk group had a significantly higher proportion of B cells, neutrophils, and NK cells in the TIME, and significantly enhanced type-I IFN response, type-II IFN response, and cytolytic activity. Furthermore, the risk model could predict the sensitivity of HCC patients to sorafenib and ICI treatment. The results indicated that patients in the low-risk group were more sensitive to sorafenib than those in the high-risk group. Conversely, the high-risk group was more sensitive to immunotherapy. Importantly, these results may contribute to personalised immunotherapy and targeted therapy for patients with HCC. Additionally, more patients in the high-risk group had TP53 somatic mutations than in the low-risk group (41% vs. 15%). TP53 mutations affect the cell cycle in approximately 30% of all HCC cases, and patients with this mutation tend to have a poor prognosis (Villanueva, 2019). Interestingly, the effect of TP53 mutations on ferroptosis showed different results depending on the mutation site Jennis et al., 2016). Furthermore, TP53 gene mutations can affect the recruitment and activity of bone marrow cells and T cells, leading to immune evasion and the promotion of cancer progression (Blagih et al., 2020). P53 functions in immune cells, leading to a variety of outcomes that can hinder or support tumour development (Blagih et al., 2020). Thus, a higher frequency of TP53 mutations in the high-risk group of HCC patients may be associated with the status of ferroptosis and an immunosuppressive phenotype. We also found that the FRGs related to our signature were all subject to varying degrees of CNVs, with the highest frequency of loss in the RB1 gene and the highest frequency of gains in YY1AP1. Although "RB1 loss of function" was present in <30% of HCC samples, both genes play an important role in HCC prognosis and treatment.
GSEA analysis results showed that the FR-related signalling pathways, such as P53 and TGF-β, were enriched in the high-risk group significantly, whereas fatty acid metabolism and drug metabolism cytochrome P450 signalling pathways were enriched in the low-risk group (Chen et al., 2021a;Chen et al., 2021b). Immune-related signalling pathways, such as mismatch repair, the Notch, P53, PI3K-AKT-MTOR, TGF-β, and WNT-β-Catenin pathways were significantly enriched in the high-risk group. These results provide further evidence for differences in ferroptosis and TIME characteristics between the two groups, which may serve as new therapeutic targets.
In the constructed FR-related lncRNA risk model, it was demonstrated that the expression of MKLN1-AS and LINC01224 was significantly upregulated in cell lines and HCC tissues and that the higher expression of MKLN1-AS and LINC01224 was associated with poorer prognosis (Dan Gong et al., 2020;Gao et al., 2020). In addition, MKLN1-AS, LINC01063, and PRRT3-AS1 can act as predictors of prognosis in patients with HCC through autophagy-related or immunerelated pathways (Deng et al., 2020;Kong et al., 2020;Yang et al., 2021). However, the roles of MKLN1-AS, LINC01063, PRRT3-AS1, and LINC01224 in the ferroptosis pathway have not yet been reported. LNCSRLR can be used as a biomarker of prognosis in patients with laryngeal squamous cell carcinoma; however, its role in HCC is unknown (Shiqi Gong et al., 2020). The function or pathway of POLH-AS1 has not yet been reported. Furthermore, the expression levels of LNCSRLR, LINC01063, PRRT3-AS1, and POLH-AS1 in HCC tissues have not been determined. We demonstrated that the expression of LNCSRLR, LINC01063, PRRT3-AS1, and POLH-AS1 was higher in HCC tissues than in normal tissues using qPCR.
Our study has certain limitations. First, further basic experiments are needed to validate the relationship between the lncRNAs screened using co-expression analysis and ferroptosis. Our study provides the basis for further in-depth research. Second, we were unable to retrieve a dataset that simultaneously reported six-lncRNA expression levels, clinical characteristics, and survival status of HCC patients. Therefore, we did not perform an external validation of the signature.
In conclusion, we constructed a novel FR-lncRNA signature with favourable specificity and sensitivity for predicting survival time in patients with HCC. A nomogram constructed using age, clinical TNM staging, and risk scores for the six FR-lncRNAs can be a simple tool to predict the survival time of HCC patients. More importantly, our signature can also predict the efficacy of immunotherapy and targeted therapies, which is important for reducing patient suffering, improving the effectiveness of drug treatment, and saving healthcare resources.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the ethic committee of the Chinese PLA General Hospital (Approval No. S2018-111-01). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
ZZ collected and analysed the data. ZZ and WZ were the major contributors in writing the manuscript. YW and XG provided technical support. WZ and CL provided guidance and advice for this article. SL, who provided the financial and ideas support, and critically reviewed the manuscript, was our corresponding author. SL, TW, and BH participated in the operation, took care of patients, and communicated with patients. ZZ and WZ were responsible for writing the manuscript. SL was responsible for the analysis of data, data interpretation, and revision. All authors read and approved the final manuscript.

FUNDING
This study was funded by the National Key R&D Program of China (2017YFA0103003 and 2017YFA0103000), and the National Natural Science Foundation of China (81670950).

ACKNOWLEDGMENTS
We appreciate the TCGA and FerrDb database for providing their platforms and contributors for uploading their meaningful datasets. We appreciate the contributions of Dr. Jin Shang to provide technical guidance and assistance.