We obtained gene sequence data of 374 patients with HCC and 50 normal samples as well as clinical information from the TCGA database (https://cancergenome.nih.gov/). We obtained metabolism-related genes from genecards (https://www.genecards.org/). The downloaded data from TCGA were processed and analyzed by Perl (https://strawberryperl.com/, version: 5.38.2) and R (https://www.r-project.org/, version: 4.2.0) software. The mRNA expression of these cholesterol metabolism-related genes in normal and HCC tissues was assessed using the “Wilcox.test” method. The R package “ggboxplot” was applied to plot the boxplot.

$|$Log${}_2$(FC)$|$ $>$1 and False discovery rate (FDR) 0.05 were the thresholds for differentially expressed genes (DEGs). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were performed using “ggplot2” package [16].

Screening of lncRNAs associated with cholesterol metabolism-related genes was conducted using the “limma” package with p 0.001 and correlation coefficient filter $>$0.3. Univariate Cox regression analysis was applied to obtain lncRNAs that are correlated with the outcome of HCC (p 0.05). Then a metabolism-related lncRNA signature was created by the least absolute shrinkage and selection operator (LASSO). Subsequently, multivariate Cox regression analysis was performed to obtain the lncRNAs that were used to construct a cholesterol metabolism-related lncRNA prediction model (p 0.05). The formula was calculated as follows:

$$Riskscore=\sum _{i=1}^n\left(Coe{f}_i\times x_i\right)$$

The x denoted the expression value of the cholesterol metabolism-related lncRNAs, and Coef represented the coefficient value. We used the formula to calculate the risk score of each patient. The lncRNA and mRNA interactions graph was produced by cytoscape software. The cytoscape related analysis were followed by a k-means cluster analysis to find associated protein subgroups. Patients with HCC were categorized into high or low-risk groups based on the median value. The two risk groups’ survival status were assessed by KM and tested using the log-rank test. We calculated the patient’s Youden index using the formula Youden = TP – FP. In order to assess the accuracy of the prognostic model, receiver operating characteristic (ROC) curves and calibration curves were performed. Survival curves of the low- and high-risk groups under different clinicopathological conditions were explored by the “survminer” and “survival” packages. Moreover, we internally validated the prognostic model by dividing the downloaded TCGA data into training and validation sets. The Spearman correlation analyses were further performed to validate the between relationship between risk score and clinical characteristics.

Clinicopathologic characteristics such as grade, stage, gender, T-, N-, and M-stage combined with risk scores together constructed a nomogram that could predict survival in HCC patients. Calibration curves were applied to check if predicted survival matched the actual survival [17].

We used GSEA (https://www.gsea-msigdb.org/gsea/index.jsp, Version 4.1.0 [18]) to explore which pathway genes were predominantly enriched in the high-risk group. Nominal p 0.05 and FDR 0.25 were regarded as the thresholds [19].

We employed the “GSVA” package for single-sample gene GSEA (ssGSEA) to calculate the activities of immune-correlated pathways and the infiltration scores of immune cells. Immune checkpoint expression analysis and mapping was performed by the “reshape2”, “limma”, “ggpubr”, and “ggplot2” packages [20].

The half-maximal inhibitory concentration (IC50) of commonly used chemotherapy drugs used in the clinical treatment of HCC was calculated. We compared the IC50 values between low- and high-risk groups using the Wilcoxon signed-rank test. The “limma”, “ggpubr”, “pRRophetic” and “ggplot2” packages [21] were used to predict which high-risk and low-risk groups have different susceptibilities to the drug.

We performed PCA analysis [22] by the “scatterplot3d” and “limma” packages to investigate the distribution of patients with different risk scores.

Statistical analyses were conducted by R software (Version 4.2.0, https://www.r-project.org/). We used the Wilcoxon test to assess the expression of cholesterol metabolism-related DEGs in normal and cancer tissues. We used the univariate Cox regression analysis to explore the association between OS and cholesterol metabolism-related lncRNAs. Multivariate Cox analysis was employed to construct the signature. Furthermore, we used the log-rank test and KM method to explore the patients’ OS different groups. The ROC curves and the AUC values were determined using the “survivalROC” package.

We obtained 59 genes correlated to cholesterol metabolism that were differentially expressed in normal versus HCC tissues, including 41 up-regulated genes and 18 down-regulated genes (Fig. 1A). GO biological process analysis illustrated that DEGs were primarily enriched in the processes of alcohol metabolic process, lipid localization, cholesterol metabolic process, steroid metabolic process, lipid transport and cholesterol transport. GO cellular components category showed the metabolism-related DEGs were mainly located at the lipoprotein particle, plasma lipoprotein particle, high-density lipoprotein particle and protein lipid complex. GO molecular function analysis showed that metabolism-related DEGs were mainly enriched in steroid hydroxylase activity, steroid binding, steroid binding, tetrapyrrole binding and oxidoreductase activity (Fig. 1B). Subsequently, we performed KEGG pathway analysis on these genes, and found that cholesterol metabolism-related DEGs were primarily enriched in the steroid hormone biosynthesis, cholesterol metabolism, metabolism of arachidonic acid, lipid and atherosclerosis, linoleic acid metabolism and steroid biosynthesis pathways (Fig. 1C).

Fig. 1.

Construction of the cholesterol metabolism-related 6 long noncoding RNA (lncRNA) prognostics signature in patients with hepatocellular carcinoma (HCC). (A) Cholesterol metabolism-related differentially expressed genes (DEGs) in normal and HCC tissues. (B) Gene Ontology (GO) analysis of cholesterol metabolism-related DEGs. (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of cholesterol metabolism-related DEGs. (D) 6 lncRNAs (AC124798.1, AL031985.3, AC103760.1, NRAV, WAC-AS1 and AC022613.1) that construct the prognostic signature. (E) The lncRNA-mRNA protein-protein interaction (PPI) network of the 6 lncRNAs. (F) Sankey diagram of the 6 lncRNAs.

We obtained 455 lncRNAs that associated with the cholesterol metabolism-related DEGs (Supplementary Table 1). A univariate Cox regression analysis demonstrated that 52 of these 455 lncRNAs were related to HCC patients’ prognosis (Supplementary Fig. 1). By using multivariate Cox regression analysis, 6 cholesterol metabolism-related lncRNAs were identified (AC124798.1, AL031985.3, AC103760.1, negative regulator of antiviral response (NRAV), WAC antisense RNA 1 (WAC-AS1) and AC022613.1) to construct a predictive signature. These 6 lncRNAs were up-regulated in patients with HCC (Fig. 1D). Subsequently, we constructed the network map of the 6 lncRNAs with associated mRNAs using cytoscape software (version:3.8.2, https://cytoscape.org/) (Fig. 1E). We used an additional k-means cluster analysis to divide the associated protein into 3 subgroups (Supplementary Fig. 2). AC124798.1 was co-expressed with SLC25A19, SMARCA4 and DYNC2LI1. AL031985.3 was co-expressed with NPC1, LBR, SLC25A19, PPARD, VDR, ALDH2, SMARCA4 and DYNC2LI1. AC103760.1 was co-expressed with CYP3A4 and CYP2C9. NRAV was co-expressed with NPC1, ALB, CYP2C9, CAV1, CYP2C8, COG2, LBR, SLC25A19, PPARD, ALDH2, SMARCA4 and DYNC2LI1. WAC-AS1 was co-expressed with NPC1, ALB, COG2, LBR, SLC25A19, PPARD, SMARCA4, DYNC2LI1, STARD3 and HSD11B2. AC022613.1 was co-expressed with HSD11B2. Fig. 1F showed that AC022613.1 and AC103760.1 are protective factors. Besides, AC124798.1, AL031985.3, NRAV and WAC-AS1 were risk factors.

Based on the formula, the risk score of each patient with HCC was measured and the patients were categorized into low- and high-risk groups based on the median. KM analysis of the OS rate of HCC patients shown that high-risk group patients had a markedly worse prognosis than low-risk group (Fig. 2A). Fig. 2B showed the risk scores of the low- and high-risk groups. Moreover, we found the higher the risk score, the more HCC patients die (Fig. 2C). We also performed the KM analysis again by taking the cutoff value of the risk scores at the maximum of the Youden-Index (Supplementary Fig. 3). Furthermore, we performed univariate Cox regression analysis to determine whether the signature could serve as an independent prognostic factor for HCC patients. The result indicated that risk score, T stage, M stage, and stage were correlated with the prognosis of patients with HCC (Fig. 2D). Moreover, risk score was also an independent predictor of HCC patients’ OS (Fig. 2E). The (area under curves) AUCs of the signature at 1, 3, and 5-years survival were 0.795, 0.756 and 0.763, respectively, indicating an accurate predictive value of the metabolism-related lncRNA predictive signature (Fig. 2F). The AUC of the risk score was 0.805, showing a better predictor of prognosis in patients with HCC than other clinicopathologic variables (Fig. 2G). Clinical factors like age and gender that do not have a significant impact on the patient’s prognosis may have smaller AUC values. Fig. 3A showed the heat map of the distribution of clinicopathologic variables and the six prognostic cholesterol metabolism-related lncRNAs in the low- and high-risk groups. AC124798.1, AL031985.3, NRAV, and WAC-AS1 were expressed higher in the high-risk group, whereas AC103760.1 and AC022613.1 were highly expressed in the low-risk group.

Fig. 2.

Validation of the cholesterol metabolism-associated lncRNA prognostic signature as a prognostic factor for HCC patients. (A) Survival analysis of HCC patients in low- and high-risk groups. (B) Distribution of the risk scores in patients with HCC. (C) The number of alive and dead HCC patients with various risk scores. Univariate Cox regression analysis (D) and multivariate Cox regression analysis (E) for the risk score and other clinical factors. (F) Receiver operating characteristic (ROC) curves at 1-, 3- and 5-years survival for the cholesterol metabolism-associated lncRNA prognostic signature. (G) The ROC curve and area under curves (AUCs) of clinicopathological variables and the risk score.

Fig. 3.

The relationship between the prognostic signature and clinicopathologic variables in the low- and high-risk groups. (A) AC124798.1, AL031985.3, NRAV, and WAC-AS1 were highly expressed in the high-risk group, whereas AC103760.1 and AC022613.1 were highly expressed in the low-risk group. (B) Nomograms containing the risk score and clinicopathologic variables predicted 1-, 3-, and 5-year overall survival (OS) in HCC patients. (C–E) Calibration curves examined the consistency between actual OS and predictions of 1-, 3-, and 5-year survival.

A nomogram including the risk score and clinicopathological variables was constructed to predict the 1, 3, and 5-years prognosis of HCC patients (Fig. 3B). The results showed that the higher the risk score, the higher the T, M, grade, and stage rating, the higher the score obtained, and the worse the 1-, 3-, and 5-year prognosis of the patient was. We performed Spearman correlation analyses to validate the between relationship between risk score and clinical characteristics. We found that the risk score was not associated with age (R = 0.078, p = 0.247) and M stage (R = 0.043, p = 0.526). The risk score was negatively correlated with gender (R = –0.143, p = 0.034), survival status (R = –0.406, p 0.001) and survival time (R = –0.365, p 0.001), and positively correlated with grade (R = 0.147, p = 0.029), T stage (R = 0.321, p 0.001) and N (R = 0.134, p = 0.046) stage (Supplementary Fig. 4). The calibration curves revealed good concordance between the predicted survival rates and the actual OS ratio at 1-, 3- and 5- years (Fig. 3C–E).

We grouped HCC patients according to age, stage, grade, T-, N- and M-stage. In HCC patients with age $\le$65 or age $\ge$65, grade 1, 2 and 3, stage 1 and 3, T1 and T3, N0 and M0, high-risk group showed a worse prognosis (Fig. 4). This suggests that this prognosis signature may predict the outcome of HCC patients in these clinicopathologic variables.

Fig. 4.

The OS of patients in low- and high-risk groups categorized according to different clinicopathologic variables. In HCC patients with age $\le$65 or $\ge$65, grade 1, 2 and 3, stage 1 and 3, T1 and T3, N0 and M0, high-risk group had a worse prognosis.

To validate the cholesterol metabolism-related lncRNA predictive signature based on the entire TCGA dataset, we randomized 342 HCC patients into two groups (n = 171). In both cohorts, patients with HCC in the low-risk group had better OS than high-risk group, consisting with the results observed across the entire dataset (Fig. 5A–D,F–I). The 1-, 3-, and 5-year survival’s AUCs in the first cohort were 0.803, 0.771, and 0.805, respectively (Fig. 5E). The AUCs for 1-, 3-, and 5-year survival in the second cohort were 0.783, 0.741, and 0.739, respectively (Fig. 5J). The results demonstrated that the cholesterol metabolism-related lncRNA predictive signature is a promising predictor of prognosis in patients with HCC.

Fig. 5.

Internal validation of the cholesterol metabolism-associated lncRNA prognostic signature. (A) Heatmap of the signature in the first internal cohort. (B) Distribution of the risk scores in patients with HCC in the first internal cohort. (C) The number of alive and dead HCC patients with various risk scores in the first internal cohort. (D) Kaplan-Meier (KM) analysis of OS rates in the first internal cohort. (E) ROC curve and AUCs at 1-year, 3-years and 5-years survival in the first internal cohort. (F) Heatmap of the signature in the second internal cohort. (G) Distribution of the risk scores in patients with HCC in the second internal cohort. (H) The number of alive and dead HCC patients with various risk scores in the second internal cohort. (I) KM analysis of OS rates in the second internal cohort. (J) ROC curve and AUCs at 1-year, 3-years and 5-years survival in the second internal cohort.

In order to investigate the possible mechanisms for the different OS of the low- and high-risk groups, we performed GSEA analysis. The results demonstrated that cell cycle, ERBB signaling pathway, pathways in cancers, lysosome, MAPK pathway, MTOR pathway, VEGF pathway, notch pathway, and RIG I like receptor pathway were significantly enriched in the high-risk group (Fig. 6). These results illustrated that high-risk patients were closely associated with immune- and tumor-related pathways, which may be related to their poor prognosis.

Fig. 6.

Genomic enrichment analysis (GSEA) analysis of low- and high-risk groups of the cholesterol metabolism-associated lncRNA prognostic signature in HCC. This lncRNA prognostic signature is associated with cell cycle, ERBB pathway, pathways in cancer, lysosome, MTOR pathway, MAPK pathway, RIG_I like receptor pathway, VEGF pathway, notch pathway, fatty acid metabolism, PPAR pathway, and biosynthesis of unsaturated fatty acids.

We then explored the correlation between risk scores and immune function and immune cell infiltrates. We found that mast cells were increased in low-risk group, while activated dendritic cells (aDCs), macrophages, Th2 cells and Treg cells were elevated in high-risk group (Fig. 7A). Compared with low-risk group, the immune function scores of MHC class I were higher in the high-risk group, while the immune function scores of type I and type II interferon (IFN) response were lower in high-risk group (Fig. 7B). Additionally, it was revealed that immune checkpoint genes expressed higher in high-risk group patients, indicating that metabolism-related lncRNA predictive signature may be a predictor for the patient’s response to ICIs (Fig. 7C).

Fig. 7.

The association between the cholesterol metabolism-associated lncRNA prognostic signature and the immune status of HCC patients. (A) Relationship between the risk score and immune cell infiltration in HCC. (B) Association of the risk score with immune pathways in HCC. (C) Association between the risk score and the immune checkpoint genes expression in HCC. ns, no significance, *p 0.05, **p 0.01, ***p 0.001.

We investigated the correlation between the prognostic signature of cholesterol metabolism-associated lncRNAs and patients’ response to clinically available drugs for HCC treatment. We found that the riskscore could indicate patient response to relevant drugs and guide more appropriate dosing regimens (Fig. 8). Moreover, PCA maps were used to visualize the distribution of patients according to the whole genome (Fig. 9A), all lncRNAs (Fig. 9B), differential expressed genes (Fig. 9C), and the model-screened cholesterol metabolism-related lncRNAs (Fig. 9D). The results indicated that the two groups were separated to a greater extent by the 6 lncRNAs risk model than by all lncRNAs, 59 DEGs and whole-genome expression profiles; that is, through this lncRNA signature HCC patients could be divided into high- and low-risk groups that differed significantly with respect to cholesterol metabolism status.

Fig. 8.

Analysis of sensitivity to therapeutic drugs in patients in low- and high-risk groups.

Fig. 9.

The construction of principal component analysis (PCA) maps. The PCA maps visualize the distribution of HCC patients according to the whole genome (A), all lncRNAs (B), differential expressed genes (C), and the model-screened cholesterol metabolism-related lncRNAs (D).