The RNA and miRNA sequence data of the normal liver tissues and HCC samples were obtained from TCGA database. All file data were downloaded using the GDC Data Transfer Tool. The demographic information (age, gender, and so on), survival endpoint (vital status, days to last follow-up and days to death), and clinical stage and grade of tumor of each sample were also downloaded.

Next, four categories of somatic mutation data of HCC patients were obtained from TCGA portal. We singled out the mutation files, which were obtained through the “SomaticSniper variant aggregation and masking” platform for subsequent analysis.

Then, the DESeq2 method with an adjusted p-value < 0.05 and | log2 fold change (FC)| > 1 setting as the threshold was employed to identify differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs), and lncRNAs (DElncRNAs) between normal and tumor samples. Taking advantage of pheatmap R package (Version: 1.0.12), a hierarchical cluster heatmap representing the expression direction and intensity of DEGs was plotted.

Firstly, interactions between lncRNA and miRNAs were predicted using the miRcode database¹ (Jeggari et al., 2012). After that, the TargetScan, mirTarBase, and miRDB (Hsu et al., 2011; Wong and Wang, 2015) databases were used to retrieve the miRNA–mRNA interaction. Finally, a ceRNA network was visualized using Cytoscape v.3.8.0.

The Entrez ID for each DEmRNA was obtained using R package “org.Hs.eg.db.” To elucidate underlying mechanisms of the hub genes related to DEmRNA in the biological process (BP), GO and KEGG function annotations were analyzed with “ggplot2,” “enrichplot,” and “clusterProfiler” packages.

All components of the ceRNA network associated with prognosis were analyzed through LASSO-penalized Cox regression to assure multifaceted models were not overfitting. Cox proportional hazards model was established using the penalized maximum likelihood algorithm. Ten-fold cross-validation was utilized to derive the best lambda to minimize the mean cross-validated error and predict the regression coefficients (β) of the multivariate Cox regression model. Then, a prognostic model including seven genes was developed, and risk score was calculated with the formula below. Risk score = βgene 1 × expression level of gene 1 + βgene 2 × expression level of gene 2 + ⋯ + βgene n × expression level of gene n. Herein, β was the regression coefficient in the multivariate Cox regression analysis as described previously (Lossos et al., 2004).

According to the previous risk formula, each HCC sample obtained corresponding risk score. All samples were stratified into high- and low-risk clusters when setting the median risk scores as the cutoff point. To visualize the correlation of risk score with clinicopathological variables, R “pheatmap” package was employed and the clinical characteristics between low- and high-risk patients were compared. Next, K-M survival curve was plotted using R package “survival” to identify prognosis difference. Moreover, time-dependent ROC curves were analyzed to validate prognosis predictive performance. Then, univariate and multivariate Cox regression analyses were performed for validity of risk signature as an independent prognostic indicator.

To identify the optimal prognostic indicator, risk score, age, gender, tumor grade, and clinicopathological stage for 1/2/3-year overall survival (OS), ROC analysis was performed (Blanche et al., 2013). To develop a quantitative prognostic pool for HCC patients, a nomogram plot integrating risk score and other clinicopathological features was constructed to predict 1–, 2–, and 3-year OS rate. Then, the calibration curve, which could present predictive validity of the nomogram, was plotted.

To elucidate the potential role of risk score in TIME contexture, seven methods including TIMER, XCELL, MCPcounter, QUANTISEQ, CIBERSORT, CIBERSORT-ABS, and EPIC were employed to estimate immune infiltration. Spearman correlation analysis was performed to investigate the association of risk score with TIME characterization.

Based on published studies, expression levels of ICB-related genes exhibited intimate interaction with immunotherapeutic efficacy (Goodman et al., 2017). Herein, six ICB-related genes were extracted: PD-1 (also known as PDCD1), PD-L1 (also known as CD274), PD-L2 (also known as PDCD1LG2), T-cell immunoglobulin domain and mucin domain-containing molecule-3 (TIM-3, also known as HAVCR2), indoleamine 2,3-dioxygenase 1 (IDO1), and cytotoxic T-lymphocyte antigen 4 (CTLA-4) (Kim et al., 2017; Nishino et al., 2017; Zhai et al., 2018). To reveal the underlying players of risk score in immunotherapy, correlation analysis of risk score with these ICB-related genes expression levels was performed.

The corresponding somatic alteration information of The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) cohort was obtained from the TCGA dataset. TMB was defined as the number of somatic, coding, base replacement, and insert–deletion mutations per megabase of the genome examined using non-synonymous and code-shifting indels under a 5% detection limit. The “maftools” R package (Mayakonda et al., 2018) was employed to detect the number of somatic non-synonymous point mutations within each sample.

The Wilcoxon test was employed to compare two groups, whereas the Kruskal–Wallis test was carried out to compare more than two groups. Survival curves were analyzed by the K-M log rank test. The chi-square test was performed to correlate the risk score subgroups with somatic mutation frequency, and the Spearman analysis computed the correlation coefficient. Results of CIBERSORT algorithm with p < 0.05 were adopted in the subsequent analysis. Two-tailed p < 0.05 deemed statistical significance. R software (version 4.0.3) was used for all statistical analyses.

DElncRNAs, DEmiRNAs, and DEmRNAs were analyzed between 374 HCC tissues and 50 adjacent normal liver samples in the TCGA database. After setting an adjusted p-value < 0.05 and | log2 FC| > 1 as cutoff threshold, a total of 136 lncRNAs (104 upregulated and 32 downregulated; Figure 1A), 128 miRNAs (107 upregulated and 21 downregulated; Figure 1C), and 2,028 protein-coding genes (1,222 upregulated and 806 downregulated; Figure 1E) were differently expressed between HCC samples and normal tissues. Clustering analysis of DElncRNAs, DEmiRNAs, and DEmRNAs suggested that HCC samples may be distinguished from normal samples according to expression profiling of DEmRNAs, DEmiRNAs, and DElncRNAs (Figures 1B,D,F). The detailed information of DElncRNAs, DEmiRNAs, and DEmRNAs are listed in Supplementary Table 1.

The interactions among DEmRNA, DElncRNA, and DEmiRNAs were predicted using different databases to construct the ceRNA regulatory network. A ceRNA regulatory network was established including 36 genes (3 lncRNAs, 12 miRNAs, and 21 mRNAs). All these genes made up the interactions of 14 lncRNA–miRNA pairs and 31 miRNA–mRNA pairs (Figure 2A and Supplementary Table 1). To explore the potential role of DEmRNAs in physiological process, GO and KEGG pathway enrichment were analyzed (Supplementary Tables 3, 4). For KEGG analysis, the top enriched terms were Cell cycle, Biosynthesis of cofactors, and Complement and coagulation cascades (Figure 2B). The result of the GO enrichment pathways presented that the DEmRNAs were primarily enriched in small-molecule catabolic process, organelle fission, and nuclear division in BPs; chromosomal region, microtubule, and collagen-containing extracellular matrix in cellular components (CCs); and coenzyme binding, tubulin binding, and organic acid binding in molecular function (MF; Figures 2C–E).

With the help of univariate Cox analysis, 19 ceRNA genes were identified with significant prognostic value (p < 0.05, Supplementary Table 5). In order to avoid overfitting the risk score model, LASSO Cox regression was conducted on the abovementioned hub genes and then recognized 11 ceRNA genes associated with prognosis in HCC (Figure 3A), and the optimal values of the penalty parameter were determined by 10-round cross-validation (Figure 3B). Next, multivariate Cox regression was performed, seven ceRNA genes (RNF24, HMMR, RAP2A, ARL2, S100A10, hsa-miR-421, and hsa-miR-326) were determined as the hub genes, all of which were considered as unfavorable prognostic indicator (all HRs > 1; Figure 3C, Supplementary Table 6). Furthermore, survival analysis showed that abnormal mRNA expression of most hub genes resulted in significant different OS times between low- and high-gene-expression subgroups (all p < 0.05; Figures 3D–J).

Subsequently, these seven hub genes were incorporated into a risk score model, and risk score was computed as follows: risk score = (0.3039 ^∗ expression value of HMMR) + (0.1663 ^∗ expression value of RNF24) + (0.2488 ^∗ expression value of RAP2A) + (0.1732 ^∗ expression value of S100A10) + (0.175 ^∗ expression value of ARL2) + (0.224 ^∗ expression value of hsa-miR-326) + (0.1399 ^∗ expression value of hsa-miR-421). Finally, each HCC sample with corresponding risk score was clustered into high-/low-risk subgroups.

First, the distributions of these seven genes with corresponding groups and samples were delineated in Figure 4A. The distributions of dot pot and risk score of survival status suggested that low-risk patients had longer OS time (Figures 4B,C). Besides, K-M survival curve demonstrated that high-risk samples presented significantly shorter OS time than patients with low-risk (p < 0.001; Figure 4D). Additionally, ROC curves were plotted, and AUC values for the 1–, 2–, and 3-year OS reached 0.784, 0.691, and 0.7, respectively (Figure 4E). Then, univariate Cox analysis pointed out that the hazard ratio (HR) of the risk score was 1.501 (95% confidence interval (CI): 1.368-1.647; Figure 4F). The results of the multivariate Cox proportional hazards model (HR = 1.469, 95% CI: 1.322–1.632; Figure 4G) supported the risk score performed as an independent prognostic indicator in HCC. These results suggested an excellent capacity of our multigene signature for clinical outcome prediction.

Firstly, the distribution of clinical variables with corresponding risk subgroups was visualized (Figure 5A). For early-grade samples and late-grade samples, risk score presented a higher trend in late-grade samples (Figure 5B). We also observed that patients with advanced stage also exhibited a significant increase in risk score (Figure 5C). Similarly, risk score was significantly elevated in T3-4 status (Figure 5D).

Stratification analysis was employed to validate whether risk score still could identify difference of prognosis when HCC patients were subgrouped into clinical variable groups. When patients were divided based on age, we found that our risk score was still predictive of patient outcomes, with higher scores indicating poorer outcomes (Supplementary Figures 2A,B). Likewise, risk score presented powerful prognostic ability for male or female patients (Supplementary Figures 2C,D), 1–2 or 3–4 clinical grade patients (Supplementary Figures 2E,F), patients in early and late stage (Supplementary Figures 2G,H), patients in T1-2 or T3-4 status (Supplementary Figures 2I,J), patients in N0 category (Supplementary Figure 2K), and patients in M0 category (Supplementary Figure 2L). The above findings, combined with the results of univariable and multivariable regression analyses, emphasized that our risk score was indeed an outstanding prognostic predictor independent from other clinical parameters.

To further investigate whether risk score was the best predictive factor among multiple clinicopathological variables, age, gender, clinical staging, tumor grade, T status, and N status were listed as candidate prognostic indicators. These clinical variables were incorporated to perform the AUC analysis for 1–, 2–, and 3-year OS, and risk score exhibited the most AUC value (Figures 6A–C). Then, a prognostic nomogram consisting of clinical stage and risk score was developed for quantitative prognosis prediction (Figure 6D). Gender, stage, grade, T status, and N status were excluded out of nomogram given, of which AUC values were less than 0.6. In addition, calibrate curves suggested excellent prognosis predictive performance of the nomogram model (Figures 6E–G).

Since HCC progression and infiltration immune cells had intrinsic and intimate connection, we further explored the potential contribution of risk signature in diversity and complexity of TIME. The result presented that risk score was remarkably and negatively correlated with populations of resting NK cells and resting memory CD4 + T cells, whereas positively related with infiltration of cancer associated fibroblast, M2 macrophage, and T regulatory cells (Tregs; Supplementary Figures 3–6). Furthermore, Spearman correlation of risk score with immune infiltration was further analyzed (Figure 7), and the detailed results are provided in Supplementary Table 7.

Given that the information on immunotherapy was not available in the TCGA-LIHC dataset, further analysis was explored for response to immunotherapy. Firstly, correlation of ICB-related gene (CD274, PDCD1, PDCD1LG2, IDO1, HAVCR2, and CTLA-4) (Kim et al., 2017; Nishino et al., 2017; Zhai et al., 2018) mRNA expression level with risk score was performed (Figure 8A). It was discovered that risk score was significantly and negatively correlated with CD274 (r = 0.19; p = 0.00028), CTLA4 (r = 0.26; p = 5.1e–07), HAVCR2 (r = 0.22; p = 2.3e–05), IDO1 (r = 0.12; p = 0.024), PDCD1 (r = 0.14; p = 0.0057), and PDCD1LG2 (r = 0.14; p = 0.0094; Figures 8B–G), suggesting that patients with high risk score may be more blocked by immune checkpoint administration.

Mounting researches have highlighted that tumor burden mutation (TMB) was associated with upregulation of CD8 + T cell infiltration, which could identify cancer cells and then execute antitumor response (Rizvi et al., 2015; McGranahan et al., 2016; Chan et al., 2019). For that, we speculated that TMB might act as a non-negligible prognostic factor of responsiveness to antitumor immunotherapy and aimed to investigate the potential interaction between risk score and TMB to uncover the hereditary variations of risk score subtype. Firstly, the TMB level was detected both in low- and high-risk score subgroups. There was no significant distinction of TMB level between the low-risk score subgroup and the high-risk score subgroup (p = 0.73, Supplementary Figure 3A). Then, the patients were assigned into distinct subtypes on the line of the TMB immune set point, as stated before (Cristescu et al., 2018). Survival curve demonstrated that high TMB value significantly suggested shorter OS time (p < 0.001, Figure 9A). Subsequent correlation analysis further validated that the TMB was positively but not significantly correlated with the risk score (R = 0.027, p = 0.61; Supplementary Figure 3B). To further explore the validity of consistent prognostic significance of risk score and TMB, we validated the cooperative effect of two indicators in prognostic prediction of HCC. As demonstrated in stratified survival curve, there was no interference of TMB status with risk score prognostic predictive performance. Risk score subgroups exhibited evident prognosis distinctions in both low and high TMB status subtypes (p < 0.001; Figure 9B). In summary, these findings indicated that risk score could act as an independent predictor and hold the potential to evaluate the clinical outcome of antitumor immunological treatment.

Besides, we explored and visualized the distribution of gene mutation in both the high-and low-risk score subtypes. The comprehensive landscape of somatic variants visualized the mutation patterns and clinicopathological features of the top 20 driver genes with the most frequent alteration (Figures 9C,D). These findings might contribute novel insight into the intrinsic connection of ceRNA and somatic variants in HCC.