Identification of potential prognostic microRNA biomarkers for predicting survival in patients with hepatocellular carcinoma

Abstract

Background

The aim of the present study was to identify potential prognostic microRNA (miRNA) biomarkers for hepatocellular carcinoma (HCC) prognosis prediction based on a dataset from The Cancer Genome Atlas (TCGA).

Materials and methods

A miRNA sequencing dataset and corresponding clinical parameters of HCC were obtained from TCGA. Genome-wide univariate Cox regression analysis was used to screen prognostic differentially expressed miRNAs (DEMs), and multivariable Cox regression analysis was used for prognostic signature construction. Comprehensive survival analysis was performed to evaluate the prognostic value of the prognostic signature.

Results

Five miRNAs were regarded as prognostic DEMs and used for prognostic signature construction. The five-DEM prognostic signature performed well in prognosis prediction (adjusted P < 0.0001, adjusted hazard ratio = 2.249, 95% confidence interval =1.491–3.394), and time-dependent receiver–operating characteristic (ROC) analysis showed an area under the curve (AUC) of 0.765, 0.745, 0.725, and 0.687 for 1-, 2-, 3-, and 5-year HCC overall survival (OS) prediction, respectively. Comprehensive survival analysis of the prognostic signature suggests that the risk score model could serve as an independent factor of HCC and perform better in prognosis prediction than other traditional clinical indicators. Functional assessment of the target genes of hsa-mir-139 and hsa-mir-5003 indicates that they were significantly enriched in multiple biological processes and pathways, including cell proliferation and cell migration regulation, pathways in cancer, and the cyclic adenosine monophosphate (cAMP) signaling pathway.

Conclusion

Our study indicates that the novel miRNA expression signature may be a potential prognostic biomarker for HCC patients.

Keywords:

• miRNA
• TCGA
• hepatocellular carcinoma
• prognosis
• biomarker

Introduction

Liver cancer is prevalent in People’s Republic of China and has become the fourth highest cancer-related death in Chinese males.¹ The majority of liver cancers are diagnosed as hepatocellular carcinoma (HCC).² Tumorigenesis is derived from environmental and genetic factors; therefore, with the completion of the human genome project, cancer genetic susceptibility has gained more and more attention. The genes in the genome that are dysregulated between tumors and normal tissues are the most promising source of diagnostic and prognostic biomarkers, as well as therapeutic targets.^3–5 Like many other cancers, HCC is characterized by the dysregulation of multiple gene networks and signaling pathways that are ordinarily involved in tissue homeostasis. Both protein-coding genes and noncoding RNA genes are involved in these genetic effects. Although the intense investigation of protein-coding genes has been the focus, microRNAs (miRNAs) are also receiving recognition as contributors to HCC.³ A miRNA is a small noncoding RNA molecule (containing about 22 nucleotides) found in plants, animals, and some viruses that functions in RNA silencing and posttranscriptional regulation of gene expression.^6,7 The miRNA expression profiling of human cancer has identified miRNA expression signatures that are associated with cancer diagnosis, prognosis, staging, and treatment.^8,9 Recently, miRNAs have been reported to play certain roles in HCC development, tumorigenesis, metastasis, clinical implication, and in the risk of hepatitis B virus (HBV) infection.¹⁰ Therefore, a comprehensive analysis of miRNA expression profiling in HCC with a reliable approach to developing clinical applications of miRNAs is highly necessary.

With the development of high-throughput sequencing technology and the completion of The Cancer Genome Atlas (TCGA) program, a large amount of genomic data is shared publicly for cancer research. TCGA has profiled and analyzed the whole-genome sequencing data at the DNA, RNA, protein, and epigenetic levels, and generated comprehensive, multidimensional maps of the key genomic changes in 33 types of cancer.¹¹ These genomic data have been made publicly available, and help the cancer research community improve the prevention, diagnosis, and treatment of cancer.¹² Comprehensive analysis of TCGA HCC whole-genome sequencing data has revealed the genomic characterization of HCC somatic mutations and developed a tumor protein p53 target gene expression signature that correlates with poor survival.¹³ However, comprehensive analysis of the miRNA prognostic signature still needs further investigation. The aim of the current study was to identify potential prognostic miRNA biomarkers for predicting survival in patients with HCC using TCGA datasets and bioinformatics analysis and to develop a miRNA expression-based prognostic signature.

Materials and methods

Data source

The HCC dataset, including Level 3 miRNA expression data, was downloaded from TCGA data portal (https://portal.gdc.cancer.gov/, accessed November 5, 2017).¹³ The corresponding clinical information was obtained from University of California, Santa Cruz Xena browser (UCSC Xena: http://xena.ucsc.edu/, accessed November 5, 2017). The raw counts of miRNA expression data of 375 tumor samples and 50 adjacent normal liver tissues were downloaded from TCGA dataset. Because the dataset in the current study was downloaded from TCGA, and data acquiring and application complied with the TCGA publication guidelines and data access policies, additional approval by an ethics committee was not needed.

Differentially expressed miRNA (DEM) screening

miRNA expression data were normalized by using the edgeR package on the R platform.¹⁴ The mean value of each miRNA greater than 0 was included in the screening of DEMs. A miRNA with a false discovery rate (FDR) < 0.05 and |log2 fold change (log2FC)| ≥ 1 of expression level between the comparison of HCC tumor and adjacent normal liver tissue was identified as DEM.

DEM expression-based prognostic signature construction

The prognostic value of miRNAs was first assessed by the univariate Cox proportional hazards regression that was performed using the survival package. These DEMs with a P-value <0.001 were regarded as prognostic DEMs whose expression levels were significantly associated with the overall survival (OS) of HCC patients. These prognostic DEMs were assessed by the “step” function to select the optimal combination and then used for prognostic model construction. To evaluate the relative contribution of these prognostic DEMs to HCC survival prediction, they were fitted into a multivariate Cox regression analysis with OS as the dependent variable. A DEM expression-based prognostic risk score model was constructed by the linear combination of the expression levels of DEMs with the multivariate Cox regression coefficient (β) as the weight. The risk score formula was as follows: risk score = expression of DEM₁ × β₁DEM₁ + expression of DEM₂ × β₂ DEM₂ +…expression of DEM_n × β_n DEM_n.^15–17 This DEM prognostic model could divide the HCC patients into high- and low-risk groups using the median risk score that was based on DEM expression. The time-dependent receiver–operating characteristic (ROC) curve was performed using the survivalROC package on the R platform to evaluate the predictive accuracy of this DEM expression-based prognostic signature in HCC OS.¹⁵

Comprehensive analysis of DEM expression-based prognostic signature

To comprehensively evaluate the prognostic model, comprehensive analysis of the DEM expression-based prognostic signature was carried out. In addition to the comparison of clinical outcomes between the low- and high-risk groups, we also investigated the prognosis predictive value of the risk score by nomogram, stratified analysis, and joint effect analysis to evaluate the association between risk score and clinical status in HCC OS. Moreover, the potential application of risk score in clinical status prediction was also explored.

Target prediction and enrichment analysis

Targetscan¹⁸ (http://www.targetscan.org/, accessed November 5, 2017), miRDB¹⁹ (http://www.mirdb.org/, accessed November 5, 2017), and miRTarBase²⁰ (http://mirtarbase.mbc.nctu.edu.tw/, accessed November 5, 2017) were used to predict miRNA target genes. The overlapping target genes among these three databases were identified as miRNA target genes and used for further enrichment analysis. Furthermore, miRNA target gene interaction networks were constructed with Cytoscape v3.4.0 (http://www.cytoscape.org/, accessed January 16, 2017).²¹ The potential function of these miRNA target genes was analyzed with the Database for Annotation, Visualization, and Integrated Discovery 6.8 (DAVID 6.8, https://david.ncifcrf.gov/home.jsp, accessed November 5, 2017)²² and visualized with the ggplot2 package.

Statistical analysis

Adjustment for multiple testing with the Benjamini–Hochberg procedure was used to control the FDR in edgeR.^23–25 Univariate analysis between clinical features and OS were compared using the log-rank test; those with P < 0.05 were entered into the multivariate Cox proportional hazards regression model for adjustment. A volcano plot and heat map were drawn by the ggplot2 package on the R platform. A value of P < 0.05 was considered statistically significant. All statistical analyses were conducted with SPSS version 20.0 (IBM Corporation, Armonk, NY, USA) and R 3.3.0.

Results

DEM screening

A total of 1881 miRNAs were obtained from the Level 3 miRNA expression dataset, and 320 miRNAs (Table S1) were identified as DEMs after edgeR filtering, which meets the criterion of FDR < 0.05 and | log2FC| ≥ 1. The volcano plot and heat map of these 320 DEMs were visualized by the ggplot2 package and shown in Figures 1 and 2, respectively.

Figure 1 Volcano plot of DEMs in HCC.

Abbreviations: DEMs, differentially expressed microRNAs; FDR, false discovery rate; HCC, hepatocellular carcinoma.

Figure 2 Heat map of 320 DEMs in HCC.

Abbreviations: DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma.

Construction of the DEM-based prognostic signature

The univariate Cox regression analysis results of DEMs are shown in Table S2, and suggest that there were 17 DEMs gained significant prognostic value for HCC with a P-value less than 0.001. By analyzing the expression level of these DEMs, we found that five miRNAs (hsa-mir-137, hsa-mir-3680-2, hsa-mir-6728, hsa-mir-3660, and hsa-mir-2682) had low expression in tumor tissue, and more than 20% of samples had an expression of zero. Therefore, only the remaining 12 miRNAs were included in the step function screening to investigate the optimal combination, and all the miRNA expression data were log2 transformed for further analysis. The following five miRNAs were screened through the step function and used to develop a miRNA prognostic model: hsa-mir-139, hsa-mir-101-2, hsa-mir-105-2, hsa-mir-9-3, and hsa-mir-5003. Kaplan–Meier and ROC curves of these prognostic DEMs are shown in Figures 3 and 4, respectively. A multivariate Cox regression analysis was used to assess the relative contribution of these prognostic DEMs in survival prediction. The risk score formula was as follows: risk score = expression of hsa-mir-139 × (−0.1795) + expression of hsa-mir-101-2 × (−0.2396) + expression of hsa-mir-105-2 × (0.0533) + expression of hsa-mir-9-3 × (0.0728) + expression of hsa-mir-5003 × (0.2640). Survival analysis of HCC clinical characteristics and risk scores are summarized in Table 1 and show that tumor stage and radical resection were significantly associated with HCC OS and were included into a multivariate Cox proportional hazards regression model for adjustment. Patients with a high-risk score have a shorter median survival time (MST) than those with a low-risk score (931 vs. 2456 days for high risk vs. low risk, Table 1, Figure 5A and B), and significantly increased risk of death (adjusted P < 0.0001, adjusted hazard ratio = 2.249, 95% confidence interval [CI] = 1.491–3.394 for OS), after being adjusted for tumor stage and radical resection. Time-dependent ROC analysis with the survivalROC package demonstrated that this miRNA expression-based prognostic signature also performed well in HCC OS prediction; the area under the curve (AUC) of the time-dependent ROC curve was 0.765, 0.745, 0.725, and 0.687 for 1-, 2-, 3-, and 5-year survival (Figure 5C). The expression distribution of these miRNAs between tumor tissue and adjacent normal liver tissue is shown in Figure 6A, whereas the expression distribution in low- and high-risk groups is shown in Figure 6B. We also explored the potential use of this risk score in the status of clinical characteristic prediction, and the ROC curves are shown in Figure 7A–G. This 5-miRNA prognostic signature performed well in predicting HCC tumor stage (P = 0.001, AUC = 0.623, 95% CI = 0.552–0.694; Figure 7A), histologic grade (P < 0.001, AUC = 0.650, 95% CI = 0.593–0.706; Figure 7B), serum AFP (P < 0.001, AUC = 0.665, 95% CI = 0.590–0.740; Figure 7C), microvascular invasion (P = 0.011, AUC = 0.587, 95% CI = 0.519–0.654; Figure 7D), and radical resection (P = 0.003, AUC = 0.646, 95% CI = 0.558–0.735; Figure 7G).

Table 1 Distribution of HCC patients’ characteristics and prognosis analysis

Variables	Events/total (n = 371)	MST (days)	Crude HR (95% CI)	Log-rank P
Age (years)				0.188
≤65	72/233	2456	1
>65	56/138	1490	1.265 (0.891–1.796)
Sex				0.283
Female	49/119	1560	1
Male	79/252	2486	0.822 (0.574–1.176)
Alcohol consumptiona				0.733
No	82/235	1852	1
Yes	40/117	1624	1.068 (0.731–1.561)
Ishak fibrosis scoreb				0.794
0 - no fibrosis	30/75	2456	1
1,2 - portal fibrosis	9/30	1372	0.964 (0.451–2.062)
3,4 - fibrous speta	6/30	NA	0.616 (0.254–1.494)
5 - nodular formation and incomplete cirrhosis	2/9	1386	0.754 (0.178–3.184)
6 - established cirrhosis	17/70	NA	0.765 (0.417–1.401)
Tumor stagec				<0.0001
I	42/172	2532	1
II	26/85	1852	1.466 (0.896–2.398)
III+IV	47/90	770	2.808 (1.844–4.278)
Histologic graded				0.858
G1	18/55	2116	1
G2	59/176	1694	1.167 (0.688–1.979)
G3	43/123	1622	1.205 (0.694–2.091)
G4	5/13	NA	1.486 (0.549–4.023)
Serum AFPe				0.788
≤400 ng/mL	62/216	2456	1
>400 ng/mL	22/64	2486	1.070 (0.654–1.750)
Radical resectionf				0.004
R0	108/324	2116	1
R1+R2+RX	17/40	837	2.082 (1.244–3.486)
MVIg				0.162
No	59/206	2131	1
Yes	36/110	2486	1.344 (0.886–2.039)
Child–Pugh scoreh				0.164
A	59/219	3125	1
B+C	9/22	1005	1.644 (0.811–3.331)
Risk score				<0.0001
Low	43/186	2456	1
High	85/185	931	2.667 (1.845–3.856)

Notes:

a Information of alcohol consumption was unavailable in 19 patients;

b information of Ishak fibrosis score was unavailable in 157 patients;

c information of tumor stage was unavailable in 24 patients;

d information of histologic grade was unavailable in 4 patients;

e information of serum AFP was unavailable in 91 patients;

F information of radical resection was unavailable in 7 patients;

g information of MVI was unavailable in 55 patients;

h information of Child–Pugh score was unavailable in 130 patients.

Abbreviations: AFP, α-fetoprotein; CI, confidence interval; HCC, hepatocellular carcinoma; HR, hazard ratio; MST, median survival time; MVI, microvascular invasion; NA, not available.

Figure 3 The Kaplan–Meier curves of five prognostic DEMs in HCC.

Notes: The order of Kaplan–Meier curves of five prognostic DEMs were as follows: hsa-mir-105-2 (A), hsa-mir-5003 (B), hsa-mir-101-2 (C), hsa-mir-139 (D), and hsa-mir-9-3 (E).

Abbreviations: DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma.

Figure 4 ROC curves of five prognostic DEMs to distinguish HCC tissue from adjacent normal liver tissue.

Notes: The order of ROC curves of five prognostic DEMs were as follows: hsa-mir-105-2 (A), hsa-mir-5003 (B), hsa-mir-101-2 (C), hsa-mir-139 (D), and hsa-mir-9-3 (E).

Abbreviations: AUC, area under the curve; CI, confidence interval; DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma; ROC, receiver–operating characteristic.

Figure 5 Prognostic risk score model analysis of five prognostic DEMs in HCC patients.

Notes: (A) From top to bottom are the risk score, patients’ survival status distribution, and five prognostic DEMs’ expression heat map for low- and high-risk groups. (B) Kaplan–Meier curves for low- and high-risk groups. (C) ROC curve for predicting survival in HCC patients by the risk score.

Abbreviations: DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma; ROC, receiver–operating characteristic.

Figure 6 Expression level of five prognostic DEMs in tumor tissue and adjacent normal liver tissue, low- and high-risk groups, respectively.

Notes: (A) Scatter plot of five prognostic DEMs’ expression level between HCC tumor tissue and adjacent normal liver tissue. (B) Scatter plot of five prognostic DEMs’ expression level between low- and high-risk groups; *P<0.05; ****P<0.0001.

Abbreviations: DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma.

Figure 7 The predictive value of the risk score for the HCC clinical parameters.

Notes: ROC curve of risk score for clinical parameters: tumor stage: cutoff by stage I+II and stage III+IV (A); histologic grade: cutoff by G1+G2 and G3+G4 (B); serum AFP: cutoff by 400 ng/mL (C); MVI: cutoff by with and without MVI (D); and Child–Pugh score: cutoff by Child A and B+C (E); Ishak fibrosis score: cutoff by 0 and 1+2+3+4+5+6 (F); radical resection: cutoff by R0 and R1+R2+RX (G).

Abbreviations: AFP, α-fetoprotein; AUC, area under the curve; CI, confidence interval; HCC, hepatocellular carcinoma; ROC, receiver-operating characteristic; MVI, microvascular invasion.

Stratified and joint effects analysis

To further investigate the association between risk score and clinical characteristics in HCC OS, comprehensive survival analysis of nomogram, stratified, and joint effects survival analysis was carried out. The stratified analysis suggests that a high-risk score significantly increased the risk of death in patients with all favorable strata except in patients with G1 grade, without fibrosis and females, and all adverse strata except in patients with Child–Pugh B/C score, without radical resection, G1 grade, tumor stage II, and with alcohol consumption (Figure 8A), indicating that this 5-miRNA prognostic signature was independent of the clinical characteristics. A nomogram was drawn by rms and its auxiliary packages based on the data of 155 HCC patients with complete clinical information in TCGA. It substantiated that the 5-miRNA prognostic signature contributes the most risk points (ranged 0–100), whereas the other clinical characteristics contribute much less (Figure 8B).

Figure 8 The relationship between risk score and clinical information.

Notes: (A) Stratified analysis of association between risk score and OS in HCC. (B) Nomogram for predicting the 1-, 3-, and 5-year event (death) with risk score and clinical information.

Abbreviations: AFP, α-fetoprotein; CI, confidence interval; HCC, hepatocellular carcinoma; MVI, microvascular invasion; OS, overall survival.

Joint effects analysis of the 5-miRNA prognostic signature and clinical parameters indicate that this prognostic signature performed well in OS prediction and in combination with the clinical parameters showed a better predictive value for HCC OS (Figure 9A–G; Table 2).

Table 2 Joint effects survival analysis of clinical factors and the DEMs’ signature risk score with OS in HCC patients

Group	Risk score	Variables	Events/total (n = 370)	MST (days)	Crude HR (95% CI)	Crude P	Adjusted HR (95% CI)	Adjusted P£
		Tumor stagea
A	Low risk	Stage I	14/100	NA	1		1
B	Low risk	Stage II	6/38	3258	1.259 (0.535–2.963)	0.598	1.265 (0.537–2.979)	0.59
C	Low risk	Stage III+IV	13/34	1622	1.888 (0.888–4.012)	0.099	1.894 (0.891–4.027)	0.097
D	High risk	Stage I	28/71	1372	2.182 (1.181–4.031)	0.013	2.160 (1.169–3.992)	0.014
E	High risk	Stage II	20/47	848	2.547 (1.368–4.743)	0.003	2.525 (1.356–4.702)	0.004
F	High risk	Stage III+IV	35/56	419	5.342 (3.112–9.170)	<0.001	4.857 (2.783–8.457)	<0.001
		Histologic gradeb
a	Low risk	G1	8/37	2131	1		1
b	Low risk	G2	18/95	3125	1.063 (0.512–2.207)	0.869	1.107 (0.496–2.472)	0.803
c	Low risk	G3+G4	10/49	NA	0.492 (0.178–1.357)	0.171	0.559 (0.193–1.618)	0.283
d	High risk	G1	10/18	581	2.029 (0.886–4.648)	0.094	1.752 (0.685–4.483)	0.242
e	High risk	G2	42/82	802	2.289 (1.126–4.655)	0.022	1.970 (0.888–4.371)	0.095
f	High risk	G3+G4	38/84	899	2.431 (1.207–4.897)	0.013	2.183 (1.001–4.760)	0.0496
		Serum AFPc
I	Low risk	≤400 ng/mL	20/121	2532	1		1
II	Low risk	>400 ng/mL	6/26	NA	0.622 (0.189–2.049)	0.435	0.439 (0.105–3.403)	0.262
III	High risk	≤400 ng/mL	42/92	1149	2.243 (1.355–3.712)	0.002	1.982 (1.154–3.403)	0.013
IV	High risk	>400 ng/mL	16/38	2486	1.990 (1.105–3.584)	0.022	1.762 (0.945–3.283)	0.075
		Microvascular invasiond
i	Low risk	No	23/119	3125	1		1
ii	Low risk	Yes	11/50	3258	1.022 (0.476–2.191)	0.956	0.727 (0.298–1.774)	0.484
iii	High risk	No	37/87	1372	2.319 (1.381–3.894)	0.001	1.958 (1.128–3.400)	0.017
iv	High risk	Yes	25/58	1149	2.968 (1.714–5.142)	< 0.001	2.236 (1.139–4.391)	0.019
		Child–Pugh scoree
①	Low risk	A	21/116	3258	1		1
②	Low risk	B+C	3/13	NA	2.348 (0.947–5.818)	0.065	2.585 (1.038–6.438)	0.041
③	High risk	A	38/100	1386	2.277 (1.348–3.848)	0.002	2.237 (1.296–3.861)	0.004
④	High risk	B+C	6/9	535	2.842 (0.849–9.515)	0.09	2.201 (0.641–7.554)	0.21
		Ishak fibrosis scoref
1	Low risk	0	13/47	3125	1		1
2	Low risk	1/2/3/4/5/6	13/77	NA	0.912 (0.414–2.010)	0.82	1.006 (0.435–2.324)	0.99
3	High risk	0	17/27	660	2.381 (1.123–5.045)	0.024	2.169 (0.977–4.819)	0.057
4	High risk	1/2/3/4/5/6	21/60	1372	1.836 (0.849–3.974)	0.123	2.079 (0.915–4.722)	0.081
		Radical resectiong
11	Low risk	R0	33/169	3125	1		1
22	Low risk	R1+R2+RX	6/14	837	1.308 (0.315–5.423)	0.712	1.501 (0.359–6.280)	0.578
33	High risk	R0	77/154	802	2.430 (1.645–3.589)	<0.001	2.263 (1.482–3.456)	<0.001
44	High risk	R1+R2+RX	11/26	837	4.154 (2.288–7.542)	<0.001	3.081 (1.533–6.191)	0.002

Notes:

a Information of tumor stage was unavailable in 24 patients;

b information of histologic grade was unavailable in 4 patients;

c information of serum AFP was unavailable in 91 patients;

d information of microvascular invasion was unavailable in 55 patients;

e information of Child–Pugh score was unavailable in 130 patients;

F information of Ishak fibrosis score was unavailable in 157 patients;

g information of radical resection was unavailable in 7 patients;

£ adjusted for tumor stage and radical resection.

Abbreviations: AFP, α-fetoprotein; CI, confidence interval; DEMs, differentially expressed microRNAs; HCC, hepatocellular carcinoma; HR, hazard ratio; MST, median survival time; NA, not available; OS, overall survival.

Figure 9 Joint effects analysis of OS stratified by risk score and HCC clinical parameters.

Notes: Joint effects analysis stratified by risk score and following clinical parameters: tumor stage (A), histologic grade (B), serum AFP (C), microvascular invasion (D), Child–Pugh score (E), Ishak fibrosis score (F), and radical resection (G).

Abbreviations: AFP, α-fetoprotein; HCC, hepatocellular carcinoma; OS, overall survival.

Target prediction and enrichment analysis

To assess the potential biological function of the five miRNAs, the miRNA target genes were predicted by three independent miRNA target gene prediction websites: Targetscan, miRDB, and miRTarBase. The candidate target genes of each miRNA predicted in all three websites were identified as miRNA target genes. Because hsa-mir-9-3, hsa-mir-101-2, and hsa-mir-105-2 were not available in the Targetscan database, only the miRNA target genes of hsa-mir-139 and hsa-mir-5003 were used for further enrichment analysis. A total of 36 genes were regarded as miRNA target genes for these two miRNAs and used for enrichment analysis and miRNA target gene interaction network construction. The miRNA target gene interaction networks are shown in Figure 10A. Functional assessment suggests that these miRNA target genes were significantly enriched in multiple biological processes, such as the regulation of cell proliferation and cell migration, cyclic adenosine monophosphate (cAMP) response element binding, response to drugs, and positive regulation of pri-miRNA transcription from RNA polymerase II promoter (Figure 10B). The Kyoto Encyclopedia of Genes and Genomes enrichment indicates that pathways in cancer, cAMP signaling pathway, proteoglycans in cancer, and focal adhesion were also significantly enriched (Figure 10C).

Figure 10 The prognostic miRNAs-target genes interactions network and their enrichment analysis results.

Notes: (A) The prognostic miRNAs-target genes interactions network. (B) GO term enrichment results. (C) KEGG enrichment results.

Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Discussion

The data mining of TCGA miRNA datasets has been widely used in cancer prognosis prediction. Similar miRNA expression-based prognostic signatures for cancer prognosis prediction have been investigated in numerous studies. For instance, miRNA expression-based prognostic signatures, which were identified by Cox proportional regression models in cancer prognosis prediction using TCGA datasets, have been developed in multiple cancers, such as pancreatic cancer,^26,27 stomach adenocarcinoma,²⁸ head and neck squamous cell carcinoma,^29,30 gliomas,³¹ glioblastoma,^32–34 cervical cancer,³⁵ bladder cancer,³⁶ breast cancer,^37,38 colon cancer,³⁹ lung cancer,^40,41 and HCC.^42–44 In addition, diagnostic^45–47 and metastasis⁴⁸ related miRNA signature biomarkers have been investigated in multiple cancers using TCGA datasets. Updates in the TCGA HCC dataset and the miRNA expression-based prognostic signatures require comprehensive survival analysis to validate; therefore, further investigations of miRNA expression prognostic signatures with a reliable method and comprehensive evaluation is necessary. The advantage of the current study was that the optimum combination of these prognostic miRNAs was investigated using the “step” function of the R platform. Additionally, we also assessed the accuracy of the prognostic signature by the survivalROC package. Furthermore, comprehensive survival analysis of our five-DEM prognostic signature has been evaluated by nomogram, stratified, and joint effects survival analysis to investigate its potential prognosis application.

Previous studies have demonstrated that miRNAs may serve as therapeutic targets and diagnosis and prognostic biomarkers in HCC.^49,50 In the current study, we used five DEMs to construct a miRNA expression-based prognostic signature to predicted HCC clinical outcome. Consistent with our finding, a study by Gu et al also observed that miR-139 was significantly downregulated in HCC tumor tissue and may be a potential diagnostic biomarker of HCC. Furthermore, low miR-139 expression in tumor tissue was significantly associated with increased risk of death.⁵¹ Similar results were also validated by other researchers, and conclude that miR-139 may be a novel prognostic biomarker of HCC.⁵² Other studies have also substantiated that miR-139 was significantly downregulated in HCC tissue^53,54 and serves as a tumor suppressor that inhibits HCC cell metastasis, progression, epithelial–mesenchymal transition, migration proliferation, and invasion by regulated genes and pathways.^53–56 The miR-139 expression level also acts as a prognostic marker in colon cancer and bile duct cancer.^57,58

Previous studies report that miR-9-3 was found to be specifically methylated in tumors of patients with non-small cell lung cancers and serves as a prognostic marker. Moreover, aberrant miR-9-3 methylation was associated with tumor progression.^59,60 Aberrant methylation of miR-9-3 has also been reported in bladder cancer,⁶¹ clear cell renal cell carcinoma,⁶² oral and oropharyngeal squamous cell carcinomas,⁶³ HCC,^64,65 and thymoma and thymic carcinoma.⁶⁶ In addition, hypermethylation of miR-9-3 was a potential diagnostic marker in bladder cancer⁶¹ and HCC,⁶⁵ whereas miR-9-3 methylation status was also significantly associated with OS in HCC⁶⁵ and increased risk of recurrence in clear cell renal cell carcinoma.⁶² In the current study, we focused on investigating the association between the expression of hsa-mir-9-3 and clinical outcome; therefore, the hsa-miR-9-3 methylation status in HCC still needs further analysis. Although the current study did not mine the methylation data of hsa-mir-9-3 using TCGA datasets, our study was the first study to reveal the value of hsa-mir-9-3 expression in HCC prognosis and diagnosis.

Hsa-mir-101-2 is encoded by a gene located in 9p21; studies by Ma et al reveal that miR-101-2-5p was involved in lipid metabolism through binding with the 3′-untranslated region of the ApoB gene in the liver of egg-laying chickens.⁶⁷ Furthermore, work by Bae et al suggests that genetic variations of miR-101-2 had a significant association with the clearance of HBV infection.⁶⁸ In cancer research, genetic variations of miR-101-2 were associated with significantly increased risk of genetic susceptibility in breast cancer.⁶⁹ Functional research in gastric cancer indicates that miR-101-2 was notably downregulated in gastric cancer tumor tissues and acts as a potential tumor suppressor in gastric cancer. Inhibited expression of miR-101-2 significantly reduced cell viability and colony formation, increased cell death, and reduced cell migration and invasion in gastric cancer cells.⁷⁰ A similar data mining study of TCGA in breast cancer demonstrated that miR-101-2 was downregulated in breast cancer tumor tissue, and may be a potential biomarker for the prognosis and diagnosis of breast cancer. Additionally, the expression level of miR-101-2 was associated with the clinical parameters of tumor progression.⁷¹ In the present study, we also observed miR-101-2 downregulation in HCC tumor tissue and low miR-101-2 expression were significantly associated with poor clinical outcome, which is consistent with previous studies. However, presently, the functional verification experiments of miR-101-2 in HCC have not been reported and still need to be further confirmed. Hsa-mir-9-3, hsa-mir-101-2, and hsa-mir-105-2 have not been reported in previous studies, and their function still needs further investigation.

Functional assessment of the target genes of hsa-mir-139 and hsa-mir-5003 suggests that these target genes were significantly enriched in the regulation of cell-based biological processes, such as cell proliferation and cell migration. Recently, miR-139 has been reported as a tumor suppressor in various types of cancer and can regulate many cancer-related genes and pathways that affect the cell function.^72–74 Previous experimental studies have substantiated that overexpression of miR-139 expression in HCC cell lines can inhibit cell proliferation, migration, and invasion through regulating multiple genes and pathways.^52,53,55,56 MiR-139 downregulation increased the invasive abilities of HCC cells in vitro and HCC metastasis in vivo.⁵⁴ These reports conclude that the function of miR-139 may be a potential tumor suppressor in HCC and may provide a therapeutic strategy for the treatment of HCC patients, especially for HCC tumor metastasis.^54–56 By reviewing the literature, we failed to retrieve the relevant reports on the function of hsa-mir-5003 in cell proliferation and migration in biological processes.

In the current study, through a comprehensive survival analysis of the prognostic signature, we have developed an accurate risk score model for HCC prognosis prediction based on the expression of five miRNAs. The risk score of patients with HCC can be divided into high- and low-risk groups, and the clinical outcome of HCC patients was significantly different between high- and low-risk groups. Time-dependent ROC analysis also suggested that this risk score model has long-term prognosis prediction and has relatively good accuracy. Stratified analysis revealed that this risk score model was an independent factor of HCC, while joint effects analysis and nomogram indicated that the risk score model performed better than traditional prognostic markers. Although this model performs well in HCC prognosis prediction, there are some limitations in the current study that still need clarification. First, because the clinical parameters from the TCGA database were not complete, such as the detailed information of treatment after the surgery was not available, we therefore could not perform a comprehensive survival analysis of the miRNA expression-based prognostic signature that considered all the potential prognostic factors of HCC. Second, patients in the current study were exclusively from a single cohort with a relatively small sample size, which may serve to bias our results. In addition, our sample size was not large enough to accurately validate the impact of the prognostic signature on OS in all strata of the stratified analysis. Therefore, our findings still need further verification in independent larger cohorts.

Despite these limitations, our current study has identified numerous DEMs, which may have potential value in HCC diagnosis and prognosis. We also constructed an effective prognostic signature to predict HCC prognosis that may have value in potential clinical applications.

Conclusion

Our study has identified 320 DEMs between HCC tumor tissue and adjacent normal liver tissue and assessed their predictive value in HCC prognosis. Five of the DEMs (hsa-mir-139, hsa-mir-101-2, hsa-mir-105-2, hsa-mir-9-3, and hsa-mir-5003) were used for prognostic signature construction based on their expression level. Survival analysis indicates that this five-DEM prognostic signature was significantly associated with HCC OS, and time-dependent ROC analysis also revealed that the prognostic signature showed a good performance in HCC survival prediction with the AUC of 0.765, 0.745, 0.725, and 0.687 for 1-, 2-, 3-, and 5-year survival, respectively. Stratified analysis substantiated that the prognostic signature serves as an independent factor of HCC and performs better in prognosis prediction than other traditional clinical indicators. Joint effects analysis demonstrated that a combination of this prognostic signature and other clinical indicators could obviously improve the performance of HCC prognosis prediction, and the nomogram also revealed that the prognostic signature contributed more than other traditional clinical indicators in HCC prognosis prediction. In addition, these five DEMs may serve as a potential independent risk factor and prognosis indicator for HCC patients and may also have potential diagnostic and prognostic value in HCC. However, due to limitations in our study, validation of our findings is necessary with independent cohorts. Moreover, additional findings obtained in our study are also worthy of further investigation.

Acknowledgments

This work was supported in part by the National Nature Science Foundation of China (81560535, 81072321, 30760243, 30460143, and 30560133), 2009 Program for New Century Excellent Talents in University (NCET), Guangxi Nature Sciences Foundation (GuiKeGong 1104003A-7), and Guangxi Health Ministry Medicine Grant (Key-Scientific Research-Grant Z201018). The present study is also partly supported by Scientific Research Fund of the Health and Family Planning Commission of Guangxi Zhuang Autonomous Region (Z2016318). We would also acknowledge the supported provided by the National Key Clinical Specialty Programs (General Surgery and Oncology) and the Key Laboratory of Early Prevention and Treatment for Regional High-Incidence-Tumor (Guangxi Medical University), Ministry of Education, China. The authors thank the contributors of TCGA (https://cancergenome.nih.gov/) and UCSC Xena (http://xena.ucsc.edu/) for sharing the HCC data on open access. In addition, we would like to acknowledge the helpful comments on this article received from our reviewers.

Supplementary materials

Table S1 DEMs identified by the edgeR package.

Abbreviation: DEMs, differentially expressed microRNAs.

Table S2 The univariate survival analysis results of the DEMs.

Abbreviation: DEMs, differentially expressed microRNAs.

Disclosure

The authors report no conflicts of interest in this work.