Development and validation of ferroptosis-related lncRNAs signature for hepatocellular carcinoma

Data source and clinical information

The RNA sequencing (RNA-seq) data together with relevant clinical data were accessed from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). An overview of the clinical information and source file of the patients with HCC can be found in Table 1 and Table S1. Notably, those patients with follow-up time greater than one month were used for the study. Totally, 374 Hepatocellular carcinoma (HCC) patients were available for further analysis. The GSE14520 dataset was acquired from GEO database (http://www.ncbi.nlm.nih.gov/geo/), containing 488 patients with HCC. A total of 259 ferroptosis-related genes (Marker: 111; Driver: 108; Suppressor: 69) were identified from FerrDb Database (Zhou & Bao, 2020) (FerrDb, http://www.zhounan.org/ferrdb/; Table S2). The TCGA dataset was utilized for training cohort while the GSE14520 dataset was validation cohort.

LncRNAs and ferroptosis-related genes data processing

The “limma” package was employed to select differentially expressed ferroptosis-related genes, which were visualized through the volcano and heatmaps. Then, we carried out functional enrichment analysis (Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)) to determine the major biological attributes. The “GOplot” package was utilized to visualize enrichment terms.

To calculate the correlation between candidate ferroptosis-related lncRNAs (FRlncRNAs) and differentially expressed ferroptosis-related genes utilizing Pearson correlation. The coefficient P < 0.001 and —R²—> 0.3 were regarded to be FRlncRNAs. Finally, Cytoscape software was utilized to draw co-expression network of prognostic FRlncRNAs and ferroptosis-related genes.

Construction of prognostic FRlncRNAs signature

First, the FRlncRNAs associated with prognosis were assessed using univariate Cox regression in training cohort. Then, FRlncRNAs with P ≤ 0.05 were included into multivariate Cox regression for the construction of FRlncRNAs signature. The formula utilized was as follows: risk score of FRlncRNAs signature = ${\sum }_{i}Coefficient\left(FRlncRNA{s}_i\right)\ast Expression\left(FRlncRNA{s}_i\right)$. To stratify patients into low- or high-risk groups, the best cut-off of the FRlncRNAs signature was identified applying receiver operating characteristic (ROC) curve at 1 year for overall survival (OS). Survival analysis between the two risk groups were assessed by Kaplan–Meier (KM) and compared using log-rank statistical methods.

Nomogram was utilized to predict 1, 3 years, and 5 years survival for the patients with HCC. ROC curves and Calibration curves were utilized to explore the accuracy of model based on training cohort. Then, we adjusted other clinical features in independent prognostic analysis in order to confirm whether the FRlncRNAs signature was an independent indicator to predict the prognosis of patient with HCC.

Validation of the FRlncRNAs signature

The GEO cohort was enrolled to verify the robustness of model established from training cohort. The FRlncRNAs signature was calculated based on validation cohort. Then, survival analysis and Cox regression were utilized to evaluate whether the FRlncRNAs signature was significantly connected with OS in validation cohort. The ROC curves were established to assess whether the novel model could accurately predict patient survival.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA, http://www.broadinstitute.org/gsea/index.jsp) was utilized to investigate functional phenotypes differences between the two risk groups (high- and low-risk groups). In this research, we carried out functional enrichment of FRlncRNAs signature, and visualized the pathway closely related to immune and tumorigenesis development. The reference gene sets contained “c7.all.v7.2.symbols.gmt [Immunologic signatures] and h.all.v7.2.symbols.gmt [cancer hallmarks]”.

Immune correlation analysis

The CIBERSORT (Charoentong et al., 2017; Newman et al., 2015), EPIC (Racle et al., 2017), ESTIMATE (Yoshihara et al., 2013), MCP counter (Shi et al., 2020), QUANTISEQ (Finotello et al., 2019), TIMER (Li et al., 2017), and single-sample gene set enrichment analysis (ssGSEA) (Yi et al., 2020) algorithms were used to infer the relative content of tumor-infiltrating immune cells (TIICs) between two risk group based on FRlncRNAs signature. The heatmap was utilized to visualize the differences of TIICs abundance under different algorithms. Besides, the correlation analysis between the abundance of TIICs and FRlncRNAs signature was utilized to reveal the potential role of FRlncRNAs signature on the immunologic features based on TIMER results.

The expression of immune checkpoint gene might be related to treatment responses of immune checkpoint inhibitors (ICIs) (Goodman, Patel & Kurzrock, 2017). Thus, we investigated six ICIs: programmed death 1 (PD-1) and its ligand 1 (PD-L1), ligand 2 (PD-L2), indoleamine 2,3-dioxygenase 1 (IDO1), cytotoxic T-lymphocyte antigen 4 (CTLA-4), and T cell immunoglobulin and mucin domain-containing protein-3 (TIM-3) in HCC (Kim et al., 2017; Nishino et al., 2017; Zhai et al., 2018). We analyzed the Spearman correlation between the ICIs and the signature, which aimed to investigate the potential role of FRlncRNAs signature in immune checkpoint blockade therapy.

Statistical analysis

All the statistical analyses were conducted using R language (version 4.0). KM analysis with log-rank test from “survival” package was applied to compare the survival difference among two risk groups. In order to evaluate the prognostic value of the novel FRlncRNAs signature, univariate and multivariate analyses were conducted through Cox proportional hazards regression model. Stratification analysis was implemented based on age (≥65 and <65 years), gender (male and female), stage (stage 1–2 and stage 3–4), grade (grade 1–2 and grade 3–4), etiology (Alcoholic liver disease and nonalcoholic fatty liver disease; HBV/HCV and non-HBV/HCV), Radiotherapy (receive radiotherapy and no receive radiotherapy), and family history (have family history of cancer and no family history). GSEA was applied to differentiate between two risk groups of functional annotations. Statistical tests were bilateral, with P value ≤ 0.05 indicated statistically significant differences.

Differentially expressed ferroptosis-related genes

After extracting the expression values of 259 ferroptosis-related genes in the patient with HCC, 69 up-regulated genes and 13 down-regulated genes were authenticated (FDR <0.05, log2FC >1; Table S3). The differentially expressed ferroptosis-related genes were visualized by volcano and heatmaps (Figs. 1A–1C). The GO enrichment revealed that these differentially expressed ferroptosis-related genes mainly participated in cellular response to chemical stress, response to oxidative stress, and carboxylic acid transport among others (Fig. 1D). The major KEGG pathways included Ferroptosis, VEGF signaling pathway, Arachidonic acid metabolism and some cancer-related signaling pathways (Fig. 1E; Table S4).

Identification of prognostic FRlncRNAs and an FRlncRNAs signature

A total of 1,184 FRlncRNAs were identified in the training cohort through the correlation analysis between differentially expressed ferroptosis-related genes and lncRNAs (—R²—>0.3 and P < 0.001; Table S5). Among them, 32 FRlncRNAs co-existed in training cohort and validation cohort were used for subsequent model construction and validation. The batch effects from different cohorts were corrected by combat function in “sva” package. Cox regression were applied to screen prognostic FRlncRNAs. In accordance with univariate cox results, eight FRlncRNAs had potential prognostic value for the patients with HCC (P < 0.05, Fig. 2A, Table S6). As shown in Table S6, eight FRlncRNAs (RHPN1 −AS1, SNHG17, DGCR11, MAPKAPK5 −AS1, PVT1, PART1, LINC00339 and HCG18) were discovered to be harmful prognostic indicators.

Subsequently, multivariate Cox regression found that three FRlncRNAs were associated with prognostic for the patients with HCC (Table S6). The three FRlncRNAs were utilized to establish a FRlncRNAs signature. The risk score was estimated using the following formula: risk score of FRlncRNAs signature = (0.89773*RHPN1-AS1) + (0.48558* MAPKAPK5-AS1) + (0.55674*PART1). The optimal cut-off point of risk score was considered to be 1.038 through ROC curve. Based on this cut-off point, the patients with HCC were classified into high- or low-risk group (Fig. 2B). The relationship between the three prognostic FRlncRNAs and co-expressed mRNA was shown in Fig. 2C and Fig. 2D. The risk score was significantly relevant to OS of patients, where OS in high-risk group possessed shorter than those in low-risk group (P < 0.001, Fig. 3A). Concurrently, the AUC of the FRlncRNAs signature was 0.719, showing great performance in contrast to other traditional clinical pathological features in predicting the prognosis of patients with HCC (Fig. 3B). The survival status plot showed that risk score of patients was inversely proportional to their survival rate. Besides, the risk heatmap demonstrated that the expression of three FRlncRNAs had positive correlation with the risk levels (Fig. 3C). The areas under the ROC (AUC) values corresponding to 1, 3, and 5 years were 0.711, 0.649 and 0.632, respectively (Fig. 3D).

The FRlncRNAs signature was an independent prognostic indicator

Univariate independent prognostic analysis revealed that risk score was a prognostic factor and significantly associated with worse survival (HR = 1.732, 95% CI [1.324–2.267]; P < 0.001) (Table 2, Fig. 4A). Moreover, after adjusting other available clinical parameters such as age, gender, stage, grade, etiology, radiotherapy, and family history, our signature still maintained an independent prognostic factor in multivariate independent analysis (HR = 1.721, 95% CI [1.280–2.314]; P < 0.001) (Table 2, Fig. 4B). As indicated in the nomogram, three-FRlncRNA-based signature was the largest contribution to OS of each period in HCC (Fig. 4C). The calibration displayed that the FRlncRNAs signature possessed high accuracy (Figs. 4D–4F).

Stratification analyses

In the training cohort, stratification analysis was conducted based on the clinicopathological features of HCC (e.g., gender, age, grade, stage, etiology, radiotherapy, and family history). As a result, the FRlncRNAs signature was still closely associated with worse survival in male, older (≥65 years) or younger (<65 years), advanced grade (Grade 3–4) or early grade (Grade 1–2) and advanced stage (Stage 3–4) or early stage (Stage 1–2) patients (all P < 0.05; Figs. 5A–5F), which indicated that FRlncRNAs signature based on risk grouping could serve as a useful tool for predicting HCC survival among each stratum of gender, age, grade, and stage. Meanwhile, this signature can also be utilized as a potential prognostic tool for HCC patients with alcoholic or non-alcoholic liver disease, no family history of cancer, HBV-positive, HBV- and HCV- negative, and no radiotherapy (Figs. 5G–5M).

Validation of the FRlncRNAs signature

In validation cohort, the risk score of FRlncRNAs signature was estimated refer to the previous formula. The cut-off point of the risk score of the validation cohort was consistent with that of the training cohort (Cutoff = 1.038). The signature was also statistically associated with OS of patients with HCC (P = 4.504e−02, Fig. 6A). Univariate independent prognostic analysis revealed that FRlncRNAs signature acted as an independent prognostic factor (HR of risk score = 1.410, 95% CI [1.054–1.887], P < 0.05, Table 2, Fig. 6B). Notably, after controlling gender, age, TNM stage, and CLIP stage, FRlncRNAs signature was no longer a prognostic factor in multivariate analysis (HR = 1.049, 95% CI [0.764–1.439], P > 0.05), which indicated that more independent cohorts needed to be included for validation (Table 2, Fig. 6C). The AUC values were more than 0.65 at 1 year, 2 years and 3 years, which showed that the FRlncRNAs signature established from training cohort had powerful accuracy and robustness (Fig. 6D).

The FRlncRNAs signature mediated DNA repair, glycolysis, MYC targets, P53 pathway and PI3K AKT mTOR signaling

GSEA was utilized to explore the potential biological mechanisms of the FRlncRNAs signature involved in HCC progression. The results of cancer hallmarks indicated that DNA repair, glycolysis, MYC targets, P53 pathway and PI3K AKT mTOR signaling were activated by the high-risk group of the FRlncRNAs signature. While adipogenesis, bile acid metabolism, coagulation, fatty acid metabolism and xenobiotic metabolism were activated by the low-risk group (Fig. 7A). Moreover, the signature also regulated many immunologic features within immune system, such as BCL6 high versus low TFH CD4 T cell down, cent memory CD4 T cell versus TH2 down etc., which indicated that the FRlncRNAs signature was implicated in immunity-related regulation (Fig. 7B).

Correlation of FRlncRNAs signature with TICIs and immune checkpoint blockade (ICB) molecule

Many immunologic features were found to be regulated by FRlncRNAs signature according to the above results. Therefore, the signature was further explored whether it was associated with TIICs based on TIMER results. The results indicated that this signature was most significantly positive correlation with immune infiltration of Neutrophil cells (COR = 0.312, P < 0.001), Myeloid cells (COR = 0.268, P < 0.001), CD4+ T cells (COR = 0.203, P < 0.001), B cells (COR = 0.130, P = 0.016), and Macrophage cells (COR = 0.124, P = 0.021; Figs. 8A–8F). The heatmap of immune responses based on CIBERSORT, EPIC, ESTIMATE, MCP counter, QUANTISEQ, TIMER and ssGSEA algorithms was displayed in Fig. 9. These findings powerfully indicated that this FRlncRNAs signature was related to immune cell infiltration in HCC.

Tumor immunotherapy utilizing ICB had gradually become a promising strategy for therapy of advanced HCC (Sangro et al., 2021; Wing-Sum Cheu & Chak-Lui Wong, 2021; Zongyi & Xiaowu, 2020). In training cohort, we carried out the association between the FRlncRNAs signature and six common ICB therapy-related targets (PD-1, PD-L1, PD-L2, TIM-3, IDO1, and CTLA-4) to explore the potential role of FRlncRNAs signature in the immunotherapy of ICB in the patients with HCC. The results showed that the FRlncRNAs signature was positively related to PD −1 (R = 0.17, P = 0.0019), CTLA-4 (R = 0.19, P = <0.001), and TIM-3 (R = 0.16, P < 0.001; except for IDO1, PD-L1, and PD-L2), revealing that the FRlncRNAs signature might play vital roles in assessment of response to ICB immunotherapy in the patients with HCC (Figs. 8G–8I and Fig. S1). Meanwhile, the expression levels of PD-1, CTLA-4 and TIM-3 were significantly higher in high-risk group contrasted with those in low- group (Figs. 8J–8L).

CONCLUSIONS

In conclusion, we identified a novel FRncRNAs signature related to prognosis of the patients with HCC, which could be utilized as a powerful tool in predicting the prognosis of patients with HCC. The FRlncRNAs signature could divide clinical characteristic subgroups according to survival. In addition, the signature was connected with ICB targets and TICIs. Hence, our study afforded a possible strategy for individualized risk stratification of the patients with HCC and evaluation response to ICB immunotherapy. The three FRncRNAs might be potential therapeutic targets of HCC.