The hepatocellular carcinoma RNA transcriptome raw count data (LIHC) and their clinicopathological characteristic information were downloaded from The Cancer Genome Atlas database (TCGA, https://portal.gdc.cancer.gov/). Those patients were included who (1) only underwent surgical resection; (2) pathologically diagnosed HCC; (3) survival longer than 30 days were included in this study. The patients from TCGA were randomly divided at a ratio of 2:1 into the training set (n = 230) and validation set (n = 118). The clinicopathological characteristic information for the TCGA LIHC including in this study is listed in Table1.

The GENCODE website is an online database that integrates human and mouse gene annotations, including protein-coding RNAs, non-coding RNAs, and pseudogenes (Uszczynska-Ratajczak et al., 2018; Frankish et al., 2019). After downloading the lncRNA annotation file from the GENCODE website and recognizing the Ensemble ID of the gene in the TCGA LIHC dataset, all the lncRNAs were extracted. Differentially expressed lncRNAs were identified with the help of R packages Deseq2 (Love et al., 2014). Cox proportional hazards regression was used to evaluate the prognostic value of these differentially expressed lncRNAs. According to previous articles, 1136 mitochondrial genes and 116 ferroptosis genes were selected (Supplementary Table S1), and their expression values were extracted from the TCGA LIHC dataset (Du and Zhang, 2020; Chen et al., 2021a; Rath et al., 2021; Zhou et al., 20212021). The R package corrplot was used to evaluate the correlations between genes and lncRNAs and calculated the correlation coefficient and p value. In this study, MF-related lncRNAs were defined as absolute correlation coefficient >0.5 and p value < 0.05.

With the help of Glmnet and survival R packages (Engebretsen and Bohlin, 2019), eight MF-related prognostic lncRNAs were identified by least absolute shrinkage and selection operator (LASSO)-Cox regression analysis and RT-qPCR from 48 MF-related lncRNA and a MF-related LncRNA prognostic signature (LPS) was constructed (Tibshirani, 1997). The MF-related LPS risk score calculation formula was listed as follows: R i s k s c o r e = ∑ i = 0 n β i ∗ x i βi was the coefficient and the Xi was the expression of lncRNAi. Based on this formula, the risk score was calculated for each HCC patient.

Multiple Cox regression was utilized to identify the independent risk factors. The discrimination performance of the MF-related LPS was assessed by the area under the receiver operating characteristic (ROC) curve (AUC). Log-rank tests and Kaplan–Meier analyses were performed to assess the predictive ability and stratification ability of MF-related LPS and the independent risk factors. Overfit bias was decreased by bootstrap validation including 1000 resamples.

To explore the mechanisms of different prognosis between low-risk and high-risk groups, differential analysis was performed based on the “Deseq2” package and only when the adjusted p-value <0.05 and | logarithmic fold-change (logFC)| >1.5 were statistically significant (Yang et al., 2019). Gene Ontology (GO), Kyoto Encyclopedia Genes and Genomes (KEGG), and gene set enrichment analysis of the DEGs were performed by R package cluster Profiler (Yu et al., 2012). We investigated the biological meaning of the common DEGs thorough GO analyses. It focused on the biological process and cellular components (CC). Kyoto Encyclopedia Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) helped us identify and evaluate the potential mechanisms and pathways.

To assess the relationship between immune cell infiltration and our MF-related LPS, we performed the ESTIMATE analysis, EPIC analysis, and CIBERSORT algorithm (Chen et al., 2018; Racle and Gfeller, 2020). In addition, the Wilcoxon signed rank test was used to evaluate the proportion of immune cells between the high- and low-risk groups.

Thirty pair samples of hepatocellular carcinoma and adjacent normal tissues were collected from Guangdong Provincial People’s Hospital (GDPH). The samples were immediately frozen after surgical resection and preserved in liquid nitrogen at −80°C. This study is approved by the Research Ethics Committee of Guangdong Provincial People’s Hospital.

Real-time quantitative polymerase chain reaction (RT-qPCR) was utilized to evaluate the lncRNA expression. All the samples were extracted using TRIzol reagent RNA TRIpure reagent (Aidlab Biotechnologies, RN01). cDNA was prepared using RevertAid Reverse Transcriptase (Thermos Fisher, United States). Real-time quantitative PCR was performed by the AceQ Universal SYBR qPCR Master Mix Kit. GAPDH was used for internal control. Each sample was tested in triplicate at least three times. All fold-changes were calculated through relative quantification 2^[(GAPDHCq) − (lncRNACq)]. The detailed primers sequences are provided in Supplementary Table S2.

Most statistical analyses were performed by R version 4.0.0 software (http://www.r-project.org/). The Deseq2 package was used to identify the differentially expressed lncRNA in the TCGA LIHC dataset. Cutoff values for the MF-related prognostic lncRNA and the nomogram were all defined by the receiver operating characteristic curve (ROC) analyses. Sangerbox, a useful and free online data analysis platform (http://www.sangerbox.com/tool), was used to analyze the correlation between the expression of 8 MF-related lncRNAs and clinicopathological factors. Kaplan–Meier analysis was used to compare the overall survival between the high-risk group and low-risk group. Univariate and multivariate Cox proportional hazards models were used to determine significant prognostic factors. The ROC curve was used to evaluate the predictive efficiency of the 8 MF-lncRNA. The figures of Kaplan–Meier analysis and ROC cures were drawn by Sangerbox platform. A p < 0.05 was considered statistically different.

The flow chart (Figure 1) displays the work process of this project. Based on the criteria, 348 patients from the TCGA LIHC dataset were included in our study. The clinical data showed that the training set included 158 males and 72 females, while the validation set included 78 males and 40 females. The patients’ age of the training set and validation set ranged from 16 to 90. The overall survival of all the patients was longer than 30 days. Approximately, 66.1% patients in training cohorts and 74.6% patients in validation cohorts were at tumor stage I–II. Patient clinicopathologic characteristics are listed in Table 1.

The expression value of 12,079 LncRNAs was extracted from the TCGA LIHC dataset. We identified 2248 lncRNAs (1829 upregulated and 419 downregulated) that were differentially expressed in the TCGA-LIHC dataset (Figure 2A). Univariate Cox regression analysis was performed to further estimate the impact of all the lncRNAs on the overall survival of LIHC patients in TCGA. Approximately, 404 differentially expressed and potential prognostic lncRNAs were screened out. By using R package corrplot, Pearson’s correlation analysis was performed to calculate the correlations between the expression of mitochondrial genes or ferroptosis genes and differentially expressed and prognostic LncRNA in the TCGA dataset. The lncRNA whose | Pearson R| > 0.5 and p < 0.05 was regarded as a mitochondrial-related or ferroptosis-related lncRNA. Subsequently, Venn analysis between the differential and prognostic mitochondrial-related lncRNAs and differential and prognostic ferroptosis-related lncRNAs was performed, generating the common lncRNAs (Figure 2B). Finally, we obtained 48 lncRNAs which were significantly correlated with mitochondrial-related and ferroptosis-related lncRNAs (MF-related lncRNAs, Figure 2C).

Also, these MF-related lncRNAs were identified as potential prognostic indicators and were subjected to least absolute shrinkage and selection operator (LASSO) logistic regression algorithm analysis in the training cohort. Moreover, 10-fold cross-validation was also used for cross-validation. We identified 15 MF-related prognostic lncRNAs based on the min lambda (Figures 3A,B). Some 15 MF-related prognostic lncRNAs were enrolled into multivariate Cox regression analysis to further identify the prognostic lncRNAs. Finally, AC022007.1, AC023090.1, AC099850.4, ACSL6-AS1, CYTOR, LHFPL3-AS2, AL365361.1, LINC02362, and MSC-AS1 were regarded as the MF-related and prognostic lncRNA (Figure 3C).

To further validate these prognostic lncRNA expression experimentally, 30 pairs of hepatocellular carcinoma tissues and adjacent normal liver tissues were selected from Guangdong Provincial People’s Hospital and detected the nine lncRNA expression by RT-qPCR. We found that AC023090.1, AC099850.4, and CYTOR were significantly highly expressed in HCC tissues, while the expression of AL365361.1, AC022007.1, and LINC02362 was largely decreased in HCC tissues. The expression levels of MSC-ASL, LHFPL3-AS2, and ACSL6-AS1 were upregulated in the HCC sample though without statistically significant (Figure 4).

Considering the expression of AC022007.1 in the TCGA cohort was opposite to the validation results of the GDPH cohort, a MF-related LPS was built based on the AC023090.1, AC099850.4, ACSL6-AS1, CYTOR, LHFPL3-AS2, AL365361.1, LINC02362, and MSC-AS1. The risk scores were calculated in the training cohort. The risk score of each patient could be calculated by the following formula: risk score = (−0.16416 × expression value of AC023090.1) + (0.37193 × expression value of AC099850.4) + (−0.29628 × expression value of ACSL6-AS1) + (−0.32870 × expression value of AL365361.1) + (0.34095 × expression value of CYTOR) + (−0.10062 × expression value of LHFPL3-AS2) + (−0.12666 × expression value of LINC02362) + (0.13118 × expression value of MSC-AS1). According to the best cutoff value of the risk scores, the training cohort can be divided into low-risk and high-risk groups (Figure 5A). Kaplan–Meier analysis was utilized to access the stratification power of the MF-related LPS. The OS curve indicated that compared with low-risk group, the high-risk group had poorer OS (Figure 5B). Moreover, by using the time-dependent ROC curves, the MF-related LPS showed excellent accuracy in predicting 1-,3-, and 5-year OS (1-year AUC = 0.80 3-year AUC = 0.76, 5-year OS = 0.84, Figure 5C).

Cox regression analysis was performed to identify the prognostic factors for LIHC patients. After removing the patients with unknown information and performing the univariate Cox regression analysis, MF-related LPS, tumor stage, and T stage were regarded as potential risk factors (Figure 6). All significant univariable predictors were enrolled into multivariate Cox regression analysis. Finally, the results indicated that MF-related LPS and tumor stage were independent prognostic factors for overall survival (OS) in patients with hepatocellular carcinoma who underwent surgical therapy (Table2).

Based on the MF-related LPS and the expression of each lncRNA from the validation cohort, the predicted risk scores were calculated. According to the best cutoff value from the training cohort, the validation population can also be divided into low-risk and high-risk groups (Figure 7A). Kaplan–Meier analysis validated the perfect stratification power of the MF-related LPS. The OS curve indicated that compared with the low-risk group, patients in the high-risk subgroup had worse survival (Figure 7B). As illustrated in Figure 7C, the AUC of MF-related LPS showed excellent accuracy in predicting 1-,3-, and 5-year OS (1-year AUC = 0.72, 2-year AUC = 0.69, 5-year OS = 0.64). After performing the multivariate Cox regression analysis, the MF-related LPS was validated as the independent prognostic factors for overall survival (OS) in the validation cohort (Table 3).

We next performed association analysis to determine whether the MF-related LPS was related to clinicopathological characteristics in whole TCGA HCC patients. As Table 4 indicated, there is no difference in age and gender between the high-risk group and low-risk group. But, there are more patients in the low-risk group with AJCC stage I–II, while there are more patients in the high-risk group with AJCC stage III–IV. As for T stage, 130 patients at T1 stage are in the low-risk group and 77 patients at T2–T3 stage are in the high-risk group. Also, to further investigate the prognostic values of MF-related LPS in different subtypes of the TCGA LIHC population, subgroup analysis was conducted. The results revealed that the high-risk group is related to poorer OS in patients at early AJCC stage (I–II) or advanced AJCC stage (III–IV), patients with tumor stage T1, T2, or T3, and patients with age ≥60 or <60 (Supplementary Figure S1).

To further investigate the differences between high- and low-risk groups, differential expression analysis was conducted by the Deseq2 package. A total of 408 DEGs (upregulation: 308 and downregulation: 100) between two groups were identified (p < 0.05, Log2FoldChange>1.0, Figure 8A). Moreover, the DEGs between the high-risk and low-risk groups were used to perform GO enrichment and KEGG pathway analyses. GO component enrichment analysis showed that upregulated DEGs may be involved in “intrinsic component of plasma membrane, neuron projection, synapse, transporter complex,” etc., while downregulated DEGs may be involved in “synapse, secretory granule cell body,” etc. (Figure 8B). The GO biological process indicated that upregulated DEGs were mainly enriched in “neurogenesis, cell signaling, regulation of membrane potential,” etc., while the downregulated DEGs were mainly enriched in the “G protein-coupled receptor signaling pathway, positive regulation of cell population proliferation, antimicrobial humoral response,” etc. The GO molecular function indicated that upregulated DEGs may be involved in “passive transmembrane transporter activity, gated channel activity, neurotransmitter receptor activity,” etc. The downregulated DEGs were mainly enriched in “signaling receptor binding, receptor regulator activity, and hormone activity.” KEGG pathway analyses indicated that these upregulated genes were mainly involved in “neuroactive ligand-receptor interaction, cAMP signaling pathway, nicotine addiction, and retrograde endocannabinoid signaling.” (Figure 8C). GSEA analysis indicated that the “G2M checkpoint, E2F targets, MYC targets, GLYCOSIS signaling, etc.” were activated (Figure 8D).

To assess the relationship between immune cell infiltration and our MF-related LPS, we performed the ESTIMATE analysis. As shown in Figure 9A, the high-risk group had lower ESTIMATE score and immune score but had higher tumor purity than the low-risk group. Therefore, the EPIC analysis was performed to investigate the immune cell infiltration. As shown in Figure 9B, there were several different immune cell infiltrations between the high-risk and low-risk groups. The low-risk group had much more B cells, macrophages, CD 8+ T cell, and CD 4+ T cell infiltration than the high-risk group. To further explore the correlation between the MF-related LPS and the infiltration of immune cell subtypes, the CIBERSORT algorithm was used. The results showed that the low-risk group had significant correlation with plasma cells, T cells, CD8 T cells, CD4 memory resting, while the high-risk group had significant correlation with regulatory T cells and M0 and M2 macrophages (Figure 9C).