Construction of an HCC recurrence model basedon the investigation of immune-relatedlncRNAs and related mechanisms.

Long noncoding RNAs (lncRNAs) are emerging as critical regulators of gene expression and play fundamental roles in immune regulation. Growing evidence suggests that immune-related genes and lncRNAs can serve as markers to predict the prognosis of patients with cancers, including hepatocellular carcinoma (HCC). This study aimed to contract an immune-related lncRNA (IR-lncRNA) signature for prospective assessment to predict early recurrence of HCC. A total of 319 HCC samples under radical resection were randomly divided into a training cohort (161 samples) and a testing cohort (158 samples). In the training dataset, univariate, lasso, and multivariate Cox regression analyses identified a 9-IR-lncRNA signature closely related to disease-free survival. Kaplan-Meier analysis, principal component analysis, gene set enrichment analysis, and nomogram were used to evaluate the risk model. The results were further confirmed in the testing cohort. Furthermore, we constructed a competitive endogenous RNA regulatory network. The results of the present study indicated that this 9-IR-lncRNA signature has important clinical implications for improving predictive outcomes and guiding individualized treatment in HCC patients. These IR-lncRNAs and regulated genes may be potential biomarkers associated with the prognosis of HCC.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common human malignancies and a leading cause of death worldwide., Although radical resection is the cornerstone of the therapeutic option for early-stage HCC patients, up to 70% of HCC patients present tumor recurrence within 5 years after resection. Early recurrence is one of the most prominent risk factors for HCC prognosis. The current strategies for assessing the risk of relapse include imaging (computed tomography/magnetic resonance imaging), the alpha-fetoprotein (AFP) serum biomarker, and several tumor characteristics, such as microvascular invasion (MVI) and portal vein tumor thrombus. Although these strategies provide some value for predicting the risk of recurrence, they still face great difficulty in providing an accurate evaluation of the risk of recurrence, owing to their limited sensitivity and specificity. Plasma biomarker detection and identification of patients with potential recurrence remain clinical dilemmas. Therefore, it is highly desirable to identify novel molecular biological markers to accurately predict the risk of early recurrence in HCC patients.

Long noncoding RNAs (lncRNAs) are longer than 200 nucleotides in length and have no protein coding functions. lncRNAs play an essential role in cancer development and progression., Evidence suggests that lncRNAs are engaged in epigenetic regulation, gene expression, genetic imprinting, chromatin remodeling, transcription, and posttranscriptional processes., In addition, lncRNAs play key roles in many cellular processes, such as the cell cycle, cell differentiation, and DNA repair., Emerging evidence suggests that lncRNAs, as regulators, play an important role in cancer immunity, such as antigen release, immune cell migration, immune cell infiltration, antigen presentation, and immune activation., Recent studies have suggested that some lncRNAs may serve as potential prognostic biomarkers in HCC. Several studies have reported the potential value of an lncRNA signature for predicting prognosis in HCC.16, 17, 18 However, most studies have focused on overall survival (OS). Currently, studies are lacking regarding the immune-related lncRNA signature in the early recurrence of HCC patients after radical resection.

In the present study, we screened 319 HCC patients under radical resection from The Cancer Genome Atlas (TCGA) to obtain patient clinical symptoms and the expression of lncRNAs. A total of 319 samples were randomly divided into a training cohort (161 samples) and a validation testing cohort (158 samples). In the training cohort, we used Cox regression models and lasso regression models to identify nine immune-related lncRNAs (9-IR-lncRNAs) that are tightly associated with the recurrence of HCC. We developed a model to predict the risk of recurrence (risk score) in HCC patients based on the 9-IR-lncRNA signature. The results revealed that HCC patients with higher risk scores had worse disease-free survival (DFS) than patients with lower risk scores. The finding in the training cohort was validated in the testing cohort.

This 9-IR-lncRNA signature not only improves the ability to predict prognosis in HCC patients with early recurrence and promote better clinical treatment strategies, but also helps elucidate the underlying mechanisms. Finally, we further verified that these hub lncRNAs were significantly differentially expressed between HCC cell lines by quantitative real-time PCR (qRT–PCR).

Results

Clinical characteristics and IR-lncRNA identification of 319 patients in the TCGA-HCC cohort

In the present study, 319 HCC samples with complete DFS information were included after screening 374 TCGA-HCC samples. The detailed analysis overview is shown in Figure 1. The "caret" package in R was used to randomly divide 319 patients into two groups as follows: 161 in the training cohort and 158 in the testing cohort. The clinical characteristics of HCC patients are shown in Table 1. There were no significant differences between the two groups in terms of age, gender, ethnic distribution, serum AFP level, vascular invasion, Child-Pugh grade, American Joint Committee on Cancer (AJCC) stage, recurrence status, and survival status. A total of 344 IR-lncRNAs were identified by expression correlation analysis of immune-related genes (IRGs) and lncRNAs (P < 0.001, cor >0.7).

Figure 1

Analysis overview

Gene expression and clinical data of HCC under radical resection were collected from the TCGA database. To construct a relapse risk prediction model, nine DFS-related IR-lncRNAs were selected using univariate Cox, lasso-Cox, and multivariate Cox analyses. For each HCC sample, tumor purity, TIS, and the fraction of immune cell infiltration were estimated using the ESTIMATE and the gene set variation analysis (GSVA) function packages in R. In the GSVA function, samples were scored according to the level of expression of a comprehensive set of gene signatures related to the immune response. The resulting scores were used to cluster the metastatic samples based on their immune profile.

Table 1

Clinical features of HCC patients with radical resection in each cohort

Variables	Training cohort (n = 161)/(%)	Testing cohort (n = 158)/(%)	p
Age			0.633
<50	32 (20)	33 (21)
50–59	51 (32)	41 (26)
60–69	53 (33)	53 (33)
≥70	25 (15)	31 (20)
Gender			0.629
Male	113 (70)	106 (67)
Female	48 (30)	52 (33)
AFP (/ULN)			0.182
<1	58 (36)	47 (30)
1–50	42 (26)	36 (23)
≥50	20 (12)	30 (19)
NA	41 (26)	45 (28)
Vascular Invasion			0.160
None	97 (60)	80 (51)
Yes	42 (26)	50 (32)
NA	22 (14)	28 (18)
Child-Pugh Grade			0.469
A	102 (63)	90 (57)
B	8 (5)	10 (6)
C	1 (1)	0 (0)
NA	50 (31)	58 (37)
AJCC Stage			0.175
Stage I	85 (53)	72 (46)
Stage II	36 (22)	43 (27)
Stage III	33 (20)	40 (25)
Stage IV	4 (23)	0 (0)
NA	3 (2)	3 (2)
Disease-Free Status			0.654
Disease Free	71 (44)	74 (47)
Recurred/Progressed	90 (56)	84 (53)
Survival Status			0.196
Alive	116 (72)	124 (78)
Death	45 (28)	34 (22)

Note: AFP, α-fetoprotein; AJCC, American Joint Committee on Cancer; HCC, hepatocellular carcinoma; NA, not available; ULN, upper limit of normal.

Construction of an HCC DFS prognostic model in the training cohort using nine IR-lncRNAs

In the training cohort, 100 IR-lncRNAs related to DFS were selected through univariate Cox analysis, 13 IR-lncRNAs (Figures S1A–S1C) were further screened using lasso-Cox regression, and nine IR-lncRNAs were included in the multivariate Cox regression. We used the multivariate Cox “stepwise regression method” to incorporate the nine core IR-lncRNAs related to DFS to establish a model (Figure S1D) and calculate the risk score using the following formula: risk score = , where coefi is the regression coefficient, and is the expression of nine IR-lncRNAs (Table 2). We further performed Kaplan-Meier curve survival analysis of the nine IR-lncRNAs. The results showed that high expression of AL031770, AC079328, AC011603, AL031985, AC087294, and LINC02193 significantly correlated with poor prognosis (P < 0.05, Figures S2A–S2G), while high expression of WARS2-IT1 and LINC002481 indicated a better prognosis (P < 0.1, Figures S2H and S2I).

Table 2

Nine immune-related long noncoding RNAs were identified for the construction of a prognostic model by a multivariate Cox regression analysis in training cohort

Variables	Univariate Analysis		Multivariate Analysis
	HR (95% CI)	p	coef	HR (95% CI)	p
LINC002481	0.81 (0.66–0.99)	0.044	−0.42	0.66 (0.52–0.83)	0.001
WARS2-IT1	0.81 (0.67–0.98)	0.028	−0.35	0.71 (0.58–0.87)	0.001
AC079328	1.38 (1.08–1.77)	0.011	0.26	1.29 (0.96–1.74)	0.092
AC011603	1.67 (1.19–2.35)	0.003	0.34	1.41 (0.92–2.16)	0.115
AL031770	1.42 (1.13–1.78)	0.003	0.23	1.25 (0.96–1.64)	0.101
AC087294	1.45 (1.16–1.80)	0.001	0.21	1.23 (0.97–1.56)	0.095
AL031985	1.55 (1.24–1.95)	<0.001	0.30	1.35 (1.02–1.79)	0.038
LINC02193	1.32 (1.12–1.56)	0.001	0.20	1.22 (1.00–1.49)	0.054
AL008718	1.59 (1.25–2.02)	<0.001	0.28	1.33 (1.02–1.72)	0.033

Note: CI, confidence interval; coef, coefficient; HR, hazard ratio.

HCC DFS prediction performance analysis of the 9-IR-lncRNA model

In the training cohort, with an increase in the risk score, the patient's recurrence or mortality gradually increased (Figures 2A and 2B), and the expression levels of WARS2-IT1 and LINC002481 decreased. In addition, the expression levels of AL008718, AL031770, AC079328, AC011603, AL031985, AC087294, and LINC02193 gradually increased with an increase in the risk score (Figure 2C). The 3D-principal components analysis diagram showed that the 9-IR-lncRNA signature distinguished the high-risk and low-risk groups (Figure 2D). Similar phenomena were observed in the verification cohort (Figures 2E–2H).

Figure 2

Relationship between the risk score and 9-IR-lncRNAs in the training cohort and testing cohort

(A) HCC patients were sorted by risk score. Red indicates high risk, and green indicates low risk. (B) DFS status of HCC patients. Dark blue indicates recurrence or dead, and light green indicates disease free. (C) Heatmap of the expression of 9-IR-lncRNAs and (D) PCA in the training cohort. (E–H) Risk score, DFS status, heatmap, and PCA in the testing cohort.

According to the median risk score of 1.02, the training cohort was divided into a high-risk group (n = 80) and a low-risk group (n = 81). The DFS of the high-risk group was significantly lower than that of the low-risk group (Figure 3A, P < 0.001). The median DFS time of the high-risk group was 0.8 years and that of the low-risk group was 3.8 years. Using the same threshold, the verification cohort was divided into a high-risk group (n = 74) and a low-risk group (n = 84). Similar to the training cohort, the DFS of the high-risk group in the verification cohort was significantly lower than that of the low-risk group (Figure 3B, P < 0.001). The median DFS time of the high-risk group was 1 year and that of the low-risk group was 3 years. The receiver operating characteristic (ROC) curve was used to test the model prediction performance. In the training cohort, the area under the curve (AUC) values of the risk score for 1-, 2-, 3-, and 5-year DFS were 0.784, 0.808, 0.768, and 0.797, respectively (Figure 3B). In the verification cohort, the AUC values of the risk score for 1-, 2-, 3-, and 5-year DFS were 0.773, 0.748, 0.678, and 0.679, respectively (Figure 3D). Furthermore, to test the advantages of our model, we performed ROC curve analysis and compared the 1-year, 2-year, 3-year, and 5-year AUCs of four published lncRNA recurrence models,, which demonstrated that our DFS model had a better performance than previous models. In terms of 1-year recurrence prediction, the AUC value of our model was 0.780, while previous models ranged from 0.53 to 0.610 (Gu JX et al. 3-lncRNA model, AUC = 0.563; Gu JX et al. 6-lncRNA model, AUC = 0.610; Zhang Q et al. 15-lncRNA model, AUC = 0.662; and Zhang ZJ et al. 14-lncRNA model, AUC = 0.563). Similar results were found for the 2-year, 3-year, and 5-year AUC values (Figures 3E–3H). These results indicated that the 9-IR-lncRNA model has moderate prediction accuracy. Next, we analyzed the differential expression of the 9-IR-lncRNAs in cancer and normal tissues in the entire TCGA-HCC cohort. The results indicated that AL008718, AL031770, AC079328, AC011603, AL031985, AC087294, LINC02193, and LINC002481 showed significantly high expression in tumors, while WARS2-IT1 showed low expression in tumors (Figure S3A). Similar results were found in HCC cell lines (Figure S3B). Compared with normal tissues, LINC002481 was significantly upregulated in HCC, leading to a better prognosis. Interestingly, this phenomenon is slightly contradictory, as some genes related to immune regulation have similar phenomena. For instance, CXCL11 expression is significantly upregulated in colon cancer but is associated with better prognosis; high CXCL11 expression results in a higher fraction of antitumor immune cells and a lower fraction of protumor immune cells. These findings suggested that WARS2-IT1 and LINC002481 may play an inhibitory role in tumorigenesis and the immune microenvironment.

Figure 3

Kaplan-Meier curve survival analysis and time-ROC curve analysis of the IR-lncRNA signature

(A) Kaplan-Meier curve survival analysis between the high-risk and low-risk groups of the training cohort. The red line indicates the high-risk group, and the blue line indicates the low-risk group. (B) Time-ROC curve analysis of the 9-IR-lncRNA signature in the training cohort. The red, green, blue, and purple lines represent the 1-, 2-, 3-, and 5-year DFS, respectively. (C and D) Kaplan-Meier curve survival analysis and time-ROC curve analysis of the 9-IR-lncRNA signature in the testing cohort. (E–H) Time-ROC curve analysis of the 9-IR-lncRNA signature compared with other published models. (E) One-year DFS, (F) 2-year DFS, (G) 3-year DFS (H) and 5-year DFS. The red line represents our 9-IR-lncRNA model, and the yellow line represents the 3-lncRNA model by Gu et al. The green line represents the 6-lncRNA model by Gu et al., and the blue line represents the 15-lncRNA model by Zhang et al. The purple line represents the 14-lncRNA model by Zhang et al.

Comparison of the 9-IR-lncRNA risk score to the combined clinical score and each individual clinical feature

Clinical characteristics, such as age, gender, serum AFP level, Child-Pugh grade, vascular invasion, and AJCC staging, were analyzed by univariate and multivariate analyses. In the training cohort, univariate analysis showed that the risk score (hazard ratio [HR] = 1.721, P < 0.001), vascular invasion (HR = 2.348, P < 0.001), and AJCC stage (HR = 1.940, P < 0.001) were significantly different (Figure 4A). Multivariate analysis showed that the risk score was an independent risk factor for the prediction of HCC recurrence (HR = 1.824, P < 0.001) (Figure 4B). In the testing cohort, univariate analysis showed that the risk score (HR = 1.177, P < 0.001), vascular invasion (HR = 1.711, P = 0.037), and AJCC stage (HR = 1.688, P < 0.001) were significantly different (Figure 4C). Multivariate analysis showed that risk score (HR = 1.743, P < 0.001) and AJCC stage (HR = 1.207, P = 0.004) were independent risk factors for the prediction of HCC recurrence (Figure 4D). The clinical risk score was a weighted sum of clinical characteristics, of which the weights were estimates from the multivariable Cox regression analysis. The AUC was used to evaluate the predictive performance of the 9-IR-lncRNA risk score, the clinical risk score, and each clinical feature. In the training cohort, the AUC value of the 9-IR-lncRNA risk score in the ROC curve of the 1-year DFS was the highest (AUC = 0.784) followed by clinical risk score (AUC = 0.719), AJCC stage (AUC = 0.698), AFP (AUC = 0.652), vascular invasion (AUC = 0.645), Child-Pugh grade (AUC = 0.532), gender (AUC = 0.507), and age (AUC = 0.437) (Figure 4E). The AUC values of the 2-year DFS ROC curve had the following order from high to low in the training cohort (Figure 4F): 9-IR-lncRNA risk score (AUC = 0.808) > clinical risk score (AUC = 0.702) > AJCC stage (AUC = 0.702) > vascular invasion (AUC = 0.621) > AFP (AUC = 0.607) > age (AUC = 0.542) > Child-Pugh grade (AUC = 0.525) > gender (AUC = 0.487). Similar results were found in the verification cohort. For the 1-year DFS ROC curve in the validation cohort, the order of AUC values from high to low was as follows: 9-IR-lncRNA risk score (AUC = 0.773) > clinical risk score (AUC = 0.748) > AJCC stage (AUC = 0.711) > AFP (AUC = 0.616) > vascular invasion (AUC = 0.557) > Child-Pugh grade (AUC = 0.532) > sex (AUC = 0.445) > age (AUC = 0.431) (Figure 4G). The order of AUC values from high to low for the 2-year DFS ROC curve in the validation cohort was as follows: clinical risk score (AUC = 0.763) > 9-IR-lncRNA risk score (AUC = 0.748) > AJCC stage (AUC = 0.718) > AFP (AUC = 0.585) > vascular invasion (AUC = 0.537) > age (AUC = 0.505) > Child-Pugh grade (AUC = 0.481) > gender (AUC = 0.431). These results showed that the accuracy of early liver cancer recurrence prediction of the 9-IR-lncRNA model was better than that of the clinical risk score and the existing clinical indicators, such as AJCC stage, AFP, vascular invasion, and Child-Pugh grade.

Figure 4

Univariate and multivariate Cox regression analyses of clinical characteristics and the 9-IR-lncRNA risk score in the training cohort and testing cohort

(A) Univariate and (B) multivariate Cox regression forest plots of clinical characteristics and the 9-IR-lncRNA risk score. (C) Univariate and (D) multivariate Cox regression forest plots of clinical characteristics and the 9-IR-lncRNA risk score. (E) One-year and (F) 2-year DFS ROC curve analysis of the 9-IR-lncRNA risk score and six clinical characteristics in the training cohort. The red, yellow, light green, green, light blue, blue, purple, and purplish red lines represent the 9-IR-lncRNA risk score, clinical risk score, AJCC stage. vascular invasion, Child-Pugh grade, AFP, age, and gender, respectively. (G) One-year and (H) 2-year DFS ROC curve analysis of the 9-IR-lncRNA risk score and six clinical characteristics in the testing cohort.

To further explore the potential relationships between the risk score and multiple clinicopathologic factors, correlation analysis was conducted via independent t tests. In the training and testing cohorts, the risk scores of AJCC stage II and III/IV were higher than that of stage I (P < 0.05, Figure S4A), and the risk scores of T2 and T3/4 were higher than that of T1 (P < 0.05, Figure S4B). Moreover, the risk scores of N0, N1, and NX were not significantly different in the training and testing cohorts (P > 0.05, Figure S4C). The patients with macro- and microvascular invasion had higher risk scores than patients with no vascular invasion (P < 0.05, Figure S4D). A higher risk score was more commonly detected in patients with increased AFP (P < 0.05, Figure S4E), but the risk scores of Child-Pugh grades A and B were not significantly different (P > 0.05, Figure S4F).

Construction of HCC DFS prognostic nomogram

In the entire TCGA-HCC cohort, we selected three parameters, including risk score, AJCC stage, and vascular invasion, which exhibited significant differences (P < 0.05) in univariate Cox regression analysis (both training and testing cohorts), to construct the nomogram model for predicting the combined risk of HCC recurrence (Figure 5A). The calibration chart showed that the nomogram model predicted DFS at 1, 2, and 3 years, exhibiting high consistency with the actual observed DFS (Figures 5B–5D), and the C-index of the DFS nomogram model was 0.732 (95% confidence interval 0.686–0.778). The AUC value of the time ROC curve is used to evaluate the combined risk score, risk score, AJCC stage, and vascular invasion of the nomogram model to predict the accuracy of DFS for 1, 2, and 3 years. The results showed that the highest AUC value of the 1-year DFS ROC curve was the combined risk score (AUC = 0.805) followed by the risk score (AUC = 0.782), AJCC stage (AUC = 0.658), and vascular invasion (AUC = 0.615) (Figure 5E). In the 2-year DFS ROC curve, the highest AUC value was also the combined risk score (AUC = 0.822) followed by the risk score (AUC = 0.777), AJCC stage (AUC = 0.696), and vascular invasion (AUC = 0.614) (Figure 5F). The same trend was also found in the 3-year DFS ROC curve with the combined risk score AUC value being the highest (AUC = 0.750) followed by the risk score (AUC = 0.718), AJCC stage (AUC = 0.645), and vascular invasion (AUC = 0.592) (Figure 5G). These results showed that the combined risk score of the nomogram model is generally more accurate than the risk score, AJCC stage, and vascular invasion in the prediction of HCC DFS.

Figure 5

Construction of the nomogram based on clinical characteristics and IR-lncRNA risk score in the entire TCGA cohort

(A) Nomogram for predicting 1-year, 2-year, and 3-year DFS in the entire TCGA cohort. (B–D) Calibration curves of the nomogram for the consistency between predicted and observed 1-year, 2-year, and 3-year DFS in the entire TCGA cohort. The dashed line at 45° indicates a perfect prediction, and the actual performances of our nomogram are shown as blue lines. (E–G) One-year, 2-year, and 3-year DFS ROC curve analysis of the combined risk score, 9-IR-lncRNA risk score, AJCC stage, and vascular invasion in the entire TCGA cohort.

Gene set enrichment analysis and immune-related characteristic analysis of the high- and low-risk groups of HCC

To further explore the potential molecular mechanisms of immune-related lncRNA signaling related to the risk of HCC recurrence, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set enrichment analysis (GSEA) between the high-risk and low-risk groups. The results showed that the high-risk group was mainly enriched in nine KEGG pathways, including DNA replication, spliceosome, mismatch repair, homologous recombination, nucleotide excision repair, base excision repair, proteasome, RNA degradation, and cell cycle (Figures S5A–S5I).

To explore immune-related characteristics, we compared the expression of immune checkpoint molecules, ESTIMATE ((Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) score, T cell inflammatory signature (TIS), and mRNA-based stem-like index (mRNAsi) between the high-risk and low-risk groups. Compared with the low-risk group, the high-risk group showed a trend of decreased expression of seven genes (CTLA-4, PD-L1, PD-L2, TIGIT, VISIR, LAG3, and TIM3), with two of these genes (PD-L2 and TIGIT) being significantly different (p < 0.05) (Figure 6A). Similarly, the stromal score, immune score, and ESTIMATE score in the high-risk group were significantly lower than those in the low-risk group (Figure 6B). The TIS, which is used to predict the response to anti-PD-L1 drugs, was significantly lower in the high-risk group than that in the low-risk group (Figure 6C). However, the mRNAsi in the high-risk group was significantly higher than that in the low-risk group (Figure 6D).

Figure 6

Comparison of immune checkpoint gene expression, ESTIMATE score, TIS, and mRNAsi between the high- and low-risk groups

Violin plot showing the expression of immune checkpoint genes (A), ESTIMATE score (B), TIS (C), and mRNAsi (D) between the high- and low-risk groups. The green bars represent the low-risk group, and the red bars represent the high-risk group. ∗p < 0.05, ∗∗p < 0.01 and ∗∗∗p < 0.001.

IR-lncRNAs regulate the infiltration of a variety of immune cells. We used the CIBERSORT method to calculate the infiltration fraction of 22 immune cells between the high-risk and low-risk groups, and we analyzed the differences (Figures 7A and 7B). The fraction of M0 macrophages in the high-risk group was higher than that in the low-risk group (P < 0.05), and the fraction of M1 macrophages, CD8+ T cells, resting memory CD4+ T cells, activated memory CD4+ T cells, and naive CD4+ T cells was lower in the high-risk group than in the low-risk group (P < 0.05). These results suggested that the high-risk group is dominated by an increase in suppressor immune cells and M0 macrophages and that the low-risk group is dominated by an increase in killer and antigen-presenting immune cells.

Figure 7

Immune cells infiltration landscape of HCC patients in high- and low-risk groups

(A) Heatmap of 22 immune cell ratios (log10) between high- and low-risk patients. The red and light green bars represent the IR-lncRNA low- and high-risk groups, respectively. (B) Correlation analysis of 22 immune cell infiltrates between the high-risk and low-risk groups. The red and dark green bars represent the IR-lncRNA low- and high-risk groups, respectively.

Construction of ceRNA regulatory network in HCC

To better understand the role of the nine-hub IR-lncRNAs screened out in HCC and to further clarify the interaction between IR-lncRNAs and miRNAs, we constructed a lncRNA-miRNA-mRNA related ceRNAs network. First, we retrieved 70 interacting lncRNA and miRNA pathways from two lncRNAs (WARS2-IT1 and AC087294) and 67 miRNAs from the miRcode database. First, we used TargetScan target gene data to identify a total of 163,558 pathways of 27 miRNAs targeting 17,781 mRNAs. Cytoscape software (version 3.7.3; https://www.cytoscape.org/) was used to construct the ceRNA regulatory network (Figure 8A), which consisted of two lncRNAs, 27 miRNAs, and 69 mRNAs. To ensure the reliability of the results, we then crossed the interaction pathways of the three databases (miRTarBase, miRDB, and TargetScan) and identified a total of 566 pathways of 17 miRNAs targeting 509 mRNAs (Figure 8B). Finally, Cytoscape software was used to construct a stricter ceRNA regulatory network (Figure 8C), which consisted of lncRNAs (WARS2-IT1), miRNAs (hsa-miR-24-3p, hsa-miR-27a-3p, and hsa-miR-140-5p), and target genes (SEMA7A, HDAC7, and CD34). These results suggested that SEMA7A, HDAC7, and CD34 may be core genes that affect the recurrence of liver cancer. We further analyzed the expression of the three target genes and three miRNA genes between cancer and normal samples. Compared with normal samples, the expression of hsa-miR-24-3p and hsa-miR-27a-3p was significantly lower in tumors, and the expression of hsa-miR-140-5p was significantly higher in tumors (Figure S6A). In addition, the expression of CD34, HDAC7, and SEMA7A was significantly higher in tumors compared with normal samples (Figure S6B). We verified the expression of seven ceRNA regulatory network genes in human normal hepatocellular (MIHA) cells and two HCC cell lines (SNU-739 and SNU- 368). Compared with normal hepatocellular cell lines, the expression of SEMA7A, HDAC7, and CD34 was higher in HCC cell lines, and the expression of hsa-miR-24-3p, hsa-miR-27a-3p, and hsa-miR-140-5p was lower in HCC cell lines (Figure 8D). Due to different origins and molecular backgrounds in cell lines, the expression of hsa-miR-24-3p was higher in SNU-368 cell lines than in MIHA normal hepatocellular cell lines and the SNU-739 HCC cell lines.

Figure 8

The ceRNAs regulatory network in hepatocellular carcinoma

(A) Plot of the ceRNA network based on the TargetScan database. The lncRNAs, miRNAs, and mRNAs are represented by nodes in pink, green, and blue. (B) Venn diagram of miRNA-mRNA pathways among the miRTarBase, miRDB, and TargetScan databases. (C) Plot of the hub ceRNA network based on the miRTarBase, miRDB, and TargetScan databases. (D) Bar plot of the expression of seven ceRNA network genes in a human normal hepatocellular cell line and HCC cell lines. Gray bars indicate the normal liver cell line (MIHA), the red and green bars indicate SNU-739 and SNU-368 HCC cell lines, respectively. Data are expressed as the mean ± SEM.

DFS and OS Kaplan-Meier curve survival analysis showed that the prognosis of patients was differentially correlated with the expression levels of the three hub miRNAs and three mRNAs. Patients with high expression of hsa-miR-24-3p had poor prognosis, as indicated by low DFS and OS (Figure S7A and S7D; DFS, p = 0.005; OS, p = 0.046). However, patients with high expression of hsa-miR-140-5p had low DFS and better OS (Figures S7B and S7E; DFS, p = 0.184; OS, p < 0.001). Patients with high expression of hsa-miR-27a-3p had poor prognosis as indicated by low DFS and OS (Figures S7C and S7F; DFS, p = 0.017; OS, p = 0.016). Patients with high CD34 expression had a better prognosis as indicated by a better DFS and OS (Figures S8A and S8D; DFS, p = 0.122; OS, p = 0.040). Patients with high expression of HDAC7 had poor prognosis with low DFS and OS (Figures S8B and S8E; DFS, p < 0.001; OS, p = 0.001), and patients with high expression of SEMA7A had better prognosis with better DFS and OS (Figures S6C and S6F; DFS, p = 0.001; OS, p = 0.026).

Finally, we used the GSE76427 dataset to perform DFS and OS survival analysis verification of CD34, HDAC7, and SEMA7A (Figures S7A–S7F), and the ICGC-LIRI cohort was used to verify the OS analysis of CD34, HDAC7, and SEMA7A (Figures S7G–S7I). Patients with high CD34 expression had a better prognosis as indicated by a better DFS and OS (Figures S8A, S8D, and S8G; GSE76427S cohort DFS, p = 0.122; GSE76427S OS, p = 0.040; ICGC-LIRI cohort OS, p = 0.004). Patients with high expression of HDAC7 had poor prognosis with low DFS and OS (Figures S8B, S8E, and S8F; GSE76427S cohort DFS, p = 0.014; GSE76427S OS, p = 0.932; ICGC-LIRI cohort OS, p = 0.288). Patients with high expression of SEMA7A had better prognosis as indicated by better DFS and OS (Figures S8A, S8D, and S8G; GSE76427S cohort DFS, p = 0.343; GSE76427S OS, p = 0.056; ICGC-LIRI cohort OS, p = 0.005).

Discussion

Recurrence of HCC is a critical cause of death worldwide in HCC patients due to the lack of a reliable method of prediction. Although many studies in the past few decades have been conducted to find reliable markers to predict the risk of recurrence, more research is required to construct a reliable system to assess the risk of recurrence in HCC. The majority of studies on predicting risk factors are focused on proteins encoding mRNAs or miRNAs., Currently, lncRNAs, especially IR-lncRNAs, which have not been previously studied, have generated broad interest. An increasing number of studies have revealed that the aberrant expression of lncRNAs is closely associated with the tumorigenesis, development, and metastasis of HCC and could become a potential biological marker for early diagnosis, immune therapy, and assessment of the risk of recurrence in HCC. Therefore, to further determine the value of IR-lncRNAs in the prediction of the recurrence of HCC, we employed univariate, multivariate Cox regression and lasso-Cox regression analyses to thoroughly screen the recurrence-related IR-lncRNAs in HCC patients with radical resection from TCGA database. Based on the screen, we constructed a 9-IR-lncRNA model to predict recurrence in HCC and further validated the model. The results indicated that the recurrence risk prediction model has good reliability in predicting recurrence in HCC and is a better predictor than currently used clinical markers, such as AJCC, MVI, and AFP. To further improve the prediction, we constructed a nomogram by combining the 9-IR-lncRNA risk score and clinical characteristics.

Our model differed from a previous lncRNA HCC prediction model.29, 30, 31 First, we used a gene expression correlation analysis to confirm IR-lncRNA. And followed by multiple step screening using univariate, lasso, and multivariate Cox regression analyses, which ensured the reliability of the lncRNAs that were selected for the model. In particular, the use of lasso regression allowed the variables to be simultaneously screened, excluded variables that contributed less to the model and selected the variables that had the most impact on the model. Therefore, our approach, which efficiently eliminated the collinearity of different variables and prevented overfitting of the model, is suitable for analyzing gene databases with large variables.

Second, we used multiple approaches, including GSEA, ESTIMATE, TIS, mRNAsi, and immune cell infiltration, to evaluate the 9-IR-lncRNA risk score. We found that the high-risk group had a significant increase in nine molecular pathways, including DNA replication, spliceosome, mismatch repair, homologous recombination, nucleotide excision repair, base excision repair, proteasome, RNA degradation, and cell cycle. Subsequently, we compared the differences between the high- and low-risk groups regarding the relative expression of immune checkpoint genes, immune score, stromal score, TIS, and stem-like index. We found that the high-risk group had a lower level of expression of immune checkpoint genes, lower immune score, lower stromal score, lower TIS, and a higher stem-like index. These findings indicated that high-risk patients have a higher purity of tumor tissue with more stem-like features and less potential benefit from immune therapy. Following immune infiltration analysis, we found that the high-risk group had significantly lower levels of CD8+ or CD4+ T cells. Third, while previous studies primarily focused on OS, we used DFS to evaluate the risk of recurrence. The use of DFS allowed us to evaluate the risk of recurrence in HCC patients and intervene promptly for high-risk patients after surgery, leading to improved survival post recurrence.

Finally, after predicting target genes using three databases (miRTarBase, miRDB, and TargetScan) and investigating the potential mechanisms of HCC recurrence with a ceRNA network, we constructed a ceRNA regulatory network consisting of lncRNAs (WARS2-IT1), miRNAs (hsa-miR-24-3p, hsa-miR-27a-3p, and hsa-miR-140-5p), and target genes (CD34, SEMA7A and HDAC7). Ye et al. and Yue et al. reported that lncRNA WARS2-IT1 affects the prediction of HCC recurrence., However, the mechanism behind lncRNA in HCC recurrence has not been investigated. Previous studies have indicated that the hsa-miR-24-3p miRNA is associated with many diseases, including atherosclerosis, ischemic chest pain, diabetes, and Parkinson disease. miR-24-3p plays a role in the development of multiple cancers. The expression of miR-24-3p in breast cancer has been found to be positively correlated with the expression of genes involved in hypoxia, and miR-24-3p has been reported to inhibit the growth of nasopharyngeal carcinoma by directly regulating Jab1/CSN5. Chen et al. found that exosomes from cells with miR-24-3p deletion effectively block EMT. All the data suggest that miR-24-3p is associated with poor prognosis in multiple cancers, including HCC, which is consistent with our finding that patients with high expression of miR-24-3p had lower DFS and OS than those with low expression. CD34 is a downstream target of miR-24-3p. As a characteristic marker of blood vessel endothelium, elevated levels of CD34 in cancer generally indicate increased angiogenesis in tumor tissues. Elevated angiogenesis may promote tumor growth and inhibit tumor growth by recruiting inflammatory cells to mediate inflammation and tumor immunity. Therefore, CD34 is a double-edged sword in cancer.

OIP5-AS1 regulates the HDAC7 epigenetic regulatory factor through miR-140-5p to promote metastasis of non-small cell lung cancer. miR-140-5p plays a role in the recurrence of glioblastoma, and the expression of miR-140-5p gradually decreases during the transformation of normal colon mucosa to tumor tissue with a further reduction in liver metastasis. As such, miR-140-5p has been considered a critical regulator in the development and metastasis of colon cancer. Although miR-140-5p is associated with the poor prognosis of colon cancer, the relationship between miR-140-5p and the prognosis of HCC has not been reported. HDAC7 is a histone deacetylase that regulates tumorigenesis and tumor development in multiple cancers. In thoracic cancer, HDAC7 promotes tumorigenesis through downregulation of miR-4465, leading to upregulation of EphA2. HDAC7 has become a potential therapeutic target in thoracic cancer. After analyzing TCGA data, Kim Freese et al. found that high expression levels of HDAC7 in HCC correlate with poor OS. However, there is no mechanistic study on the function of HDAC7 in HCC.

The expression of miR-27a-3p is elevated in multiple cancers, including triple-negative breast cancer, glioblastoma, and osteosarcoma, and the expression level of miR-27a-3p is correlated with OS in these cancers. From the TCGA database, Xu et al. found that the expression of miR-27a-3p is elevated in HCC, but they did not further perform prognostic analysis or investigate the related mechanism. Semaphorin 7A (SEMA7A) is a highly effective immune regulator that is expressed in active lymphocytes and bone marrow cells. SEMA7A is significantly elevated in natural killer (NK) cells when stimulated by cytokines, specifically tagging active NK cells, leading to a strong capability of interferon-γ release. In lung adenocarcinoma clinical samples, high levels of SEME7A predict poor response to EGFR-TKI therapy. In breast cancer, SEMA7A is expressed in breast cells during gland degeneration, leading to changes in macrophage biology and the formation of lymphatic vessels, thereby promoting breast cancer metastasis. Overexpression of SEMA7A in estrogen receptor-positive (ER+) cells promotes cell growth. Crump et al. reported that SEMA7A confers drug resistance to Fulvestrant in primary tumors and induces lung metastasis in ER+ breast cancer. However, the study of SEMA7A in HCC has not been reported.

In conclusion, our research showed that the 9-IR-lncRNA model predicts the risk of HCC recurrence and that the lncRNA(WARS2-IT1)-miRNA(hsa-miR-24-3p, hsa-miR-27a-3p, and hsa-miR-140 -5p)-mRNA (CD34, SEMA7A, and HDAC7) regulatory network is closely related to the prognosis of HCC, which is expected to become a potential target for HCC recurrence prediction and immunity therapy.

Materials and methods

Data acquisition and processing

RNA-sequencing (RNA-seq) data of HCC patients and clinical data of HCC samples were downloaded from the TCGA program (https://portal.gdc.cancer.gov). The RNA-seq data were normalized using the trimmed mean of M-values (TMM) method in the “edgeR” package. The mean gene expression levels with a logarithmic fold change of >1 in all samples were retained. The normalized expression levels were further transformed using log2 (x+1) transformation. The initial dataset contains 374 HCC tissues and 50 adjacent normal tissues. Samples without radical resection treatment and complete survival data were removed, resulting in 319 radical resection HCC samples with complete clinical information for subsequent analysis in our study.

Identification of IR-lncRNAs

The IRGs were obtained from the Molecular Signatures Database v7.0 (Immune system process M13664, Immune response M19817)., The IR-lncRNAs were identified by a Pearson correlation analysis between the IRGs and lncRNA expression levels in the HCC samples (Pearson correlation coefficient >0.7, p < 0.001). The “caret” package in R was used to randomly divide the samples into the training cohort (161 samples) and the testing cohort (158 samples). The expression levels of the IR-lncRNAs were generated from both the training cohort and testing cohort.

Construction of an IR-lncRNA signature and combined clinical risk score associated with DFS

The “survival,” “survminer,” and “forestplot” packages in R were used to perform univariate and multivariate Cox regression analyses. The “glmnet” package in R was used to perform a lasso regression analysis. IR-lncRNAs with a significant prognostic value (p < 0.05) were screened. In the multivariate Cox regression analysis, the IR-lncRNA signature was expressed as follows: risk score =. The median risk score served as a cutoff value to classify the patients in the training cohort into the high- and low-risk groups. The testing cohort was also divided into high- and low-risk groups using the same cutoff value. We constructed a combined clinical characteristic risk score weighted sum of those clinical characteristics, of which the weights were estimates from the multivariable Cox regression analysis. The clinical risk score formula was as follows:

Clinical_risk_score = 0.190∗AJCC_stage+0.767∗Vascular_invasion+0.500∗Child_grade+0.332∗AFP+0.128∗Age-0.188∗ Gender. The “survminer” package in R was used to plot Kaplan-Meier survival curves. The “survivalROC” package in R was used to investigate the prognostic value of the IR-lncRNA risk score over time. A two-sided log rank p < 0.05 was considered significant in the survival analysis.

Construction of recurrence risk prediction nomogram based on clinical characteristics and IR-lncRNA risk score

A prognostic model was constructed by obtaining the independent prognostic factors and IR-lncRNA risk score. A time-dependent ROC curve was employed to compare the models. The most appropriate independent prognostic factors were selected for the prognostic model to construct a nomogram in the entire TCGA cohort. The calibration plot and concordance index (C-index) were used to evaluate the capability of the nomogram to calibrate and discriminate against different factors (by a bootstrap method with 1,000 resamples).

Gene set enrichment analysis

To understand the KEGG pathways of the IR-lncRNA signature, a GSEA was used to analyze the enrichment terms in the entire TCGA cohort. GSEA software version 4.1.0 (Cambridge, MA) was used to perform the analysis, and p < 0.05 was considered statistically significant.

ESTIMATE score, TIS, mRNAsi, and immune cell infiltration

ESTIMATE is a tool to predict the tumor purity from gene signatures and calculate the following three scores: (1) stromal score, which predicts the presence of stromal cell types in tumor bulk; (2) immune score, which infers the infiltration of immune cells in tumor tissue; and (3) estimate score, which estimates the tumor purity. All patients were divided into high-risk and low-risk groups. Gene set variation analysis from the GSVA package in R was performed to obtain the immune profile of the HCC samples. This function performs a nonparametric, unsupervised analysis for estimating variation of the given gene sets through the samples in the expression matrix, returning an enrichment score for each sample. The TIS scored with the GSVA method was used to predict the putative response to an anti-PD-L1 drug. This is a genetic profile reported as an acceptable predictor of clinical response to pembrolizumab across a wide variety of tumor types. In addition, the mRNAsi, which reflects gene expression, was adopted by means of one class logistic regression machine learning to quantify cancer stem-like indices., The fraction of 22 immune cell types infiltrated in tumor tissue was assessed using the “CIBERSORT” algorithm.

Construction of the ceRNA network

We used an miRcode database to predict the interaction pairs between lncRNAs and miRNAs. mRNAs targeted by the miRNAs were retrieved from the miRDB, miRTarBase, and TargetScan databases.66, 67, 68 Candidate mRNAs were those identified by all three databases. Based on miRNA-lncRNA and miRNA-mRNA interaction pairs, the ceRNAs network was constructed and visualized using Cytoscape software (version 3.7.3; https://www.cytoscape.org/).

Expression level of seven ceRNA network genes in cell lines as detected by a qRT–PCR assay

Total RNA was extracted from the SNU-739 and SNU-368 HCC cell lines as well as the normal human hepatocellular cell line (MIHA cells). Total RNA was prepared using Total RNA Kit I (Omega, R6834-01), and cDNA was synthesized using a reverse transcription kit (TaKaRa Biotechnology, Shiga, Japan) according to the manufacturer's protocol. Quantitative real-time PCR (qRT–PCR) was performed using SYBR Green Mix (TaKaRa RR820A, Japan) and a C1000 system (Bio-Red, Hercules, CA). Target gene expression values were normalized to human GAPDH, and target miRNA expression values were normalized to human U6.

Statistical analysis

Statistical analyses were performed using SPSS software version 26.0 (Armonk, NY) and R software version 4.0.2 (Auckland City, New Zealand). The χ2 test or Fisher's exact test was used to evaluate the qualitative variables as appropriate. All comparisons between the high-risk and low-risk groups were analyzed using Wilcoxon tests, as variables were either not normally distributed or variances were not equally distributed between groups. For all tests applied, differences were considered significant when p < 0.05.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article.