This article is mentioned in:

Introduction

Hepatocellular carcinoma is the third leading cause of cancer-associated mortality (1). The most significant characteristics of hepatocellular carcinoma are aggressiveness, invasiveness and frequent recurrence (2). Several risk factors for hepatocellular carcinoma have been identified, including liver cirrhosis, alcohol abuse, steatohepatitis, obesity, diabetes, intake of the fungal metabolite aflatoxin B1 and metabolic syndromes (3–8). In addition, frequent hyper-methylation of tumor suppressor genes, including p16, suppressor of cytokine signaling 1, glutathione S-transferase pi 1, Ras association (RalGDS/AF-6) domain family member 1A and E-cadherin have been involved in the tumorigenesis of hepatocellular carcinoma (9–11). Clinically, surgical resection, interventional therapy, liver transplantation, liver-directed therapy and systemic therapy are common treatment methods (12). However, only liver transplantation and surgical resection are regarded as effective treatments. In addition, only 15% patients are eligible for effective treatments, whereas most patients present with advanced disease at diagnosis (13). Although some new therapeutic methods have been developed, the 5-year survival rate of hepatocellular carcinoma remains poor due to late diagnosis, and the survival rate is currently 7% (14). Therefore, identifying biomarkers for the early diagnosis of hepatocellular carcinoma is required.

MicroRNAs (miRNAs/miRs) play a crucial role in regulating cell proliferation, differentiation, migration and apoptosis (15,16). The expression of miRNAs has been previously investigated in hepatocellular carcinoma, and several miRNAs, such as hsa-miR-21, hsa-miR-223 and hsa-miR-122, were identified as upregulated (17,18); while certain miRNAs, such as hsa-miR-122a, hsa-miR-152 and hsa-miR-22, were downregulated in hepatocellular carcinoma tissues (19–21). A previous study demonstrated that miRNAs may be used as biomarkers and therapeutic targets for the diagnosis and treatment of hepatocellular carcinoma (22). In the present study, in order to identify potential diagnostic biomarkers, differentially expressed miRNAs and mRNAs were investigated in hepatocellular carcinoma based on The Cancer Genome Atlas (TCGA) database. A machine learning approach was used to identify the differentially expressed miRNAs with diagnostic potential for hepatocellular carcinoma. Subsequent analysis was based on these differentially expressed miRNAs with diagnostic potential.

Materials and methods

Data retrieval

TCGA (http://tcga-data.nci.nih.gov/) is a publicly funded project which consists of multidimensional data of for multiple cancer types at DNA, RNA and protein levels. In the database, the clinical data of 377 patients, the mRNA data of 371 patients and the miRNA data of 373 patients were recorded. Details of the dataset are presented in Table SI. The mean age of these patients was 59.45±13.5. In addition, the male:female ratio of the patients was 255:122. The clinical information of these patients in the TCGA is presented in Table SI. The expression data of miRNAs and mRNAs were generated using an RNA sequencing platform. Transcriptome mRNA and miRNA data of hepatocellular carcinoma were all obtained from primary solid tumors and normal solid tissues. Screening requirements for the included samples were as follows: i) Samples without clinical information were excluded; ii) samples with incomplete information about stage and survival time were excluded; and iii) samples with information about both miRNA and mRNA expression levels were reserved. According to the inclusion and exclusion criteria, the mRNA and miRNA data (including 342 cases and 50 controls) were finally used for the following integrated analysis.

Identification of differentially expressed miRNAs and mRNAs

Before identification, data of miRNAs and mRNAs were preprocessed, and miRNAs and mRNAs that were lowly expressed were deleted. miRNAs and mRNAs were considered to have low expression when the number of control samples presenting read counts value of 0 was >20% of the total case sample size or when the number of case samples presenting read counts value of 0 was >20% of the total control sample size. Principal component analysis of these miRNAs and mRNAs was subsequently conducted. Differentially expressed miRNAs and mRNAs were analyzed, as previously described (23). The false discovery rate (FDR) was obtained from multiple comparisons using The Benjamini and Hochberg method (24). Those differentially expressed miRNAs and mRNAs were identified with the criterion of FDR<0.05, abs |(log2FoldChange)|>3 and FDR<0.05, abs |(log2FoldChange)|>1, respectively. In addition, the heat map of differentially expressed miRNAs and mRNAs was performed. Clustering was analyzed using the complete-linkage method together with the Euclidean distance.

Identification of the optimal diagnostic biomarkers based on a machine learning approach

Firstly, the importance value of each differentially expressed miRNA was ranked using a random forests (RF) algorithm. Then, the optimal number of features was identified by subsequently adding one differentially expressed miRNA at a time in a top down forward-wrapper approach. Optimal differentially expressed miRNA with diagnostic value for hepatocellular carcinoma were used to establish classification models, including RF, support vector machine (SVM) and decision tree (DT). The ‘randomForests’ package in R language (https://cran.r-project.org/web/packages/randomForest/), ‘e1071’ package in R language (https://cran.r-project.org/web/packages/e1071/index.html) and ‘rpart’ package in R language (https://cran.r-project.org/web/packages/rpart/index.html) were used to establish the RF model, SVM model and DT models, respectively. Diagnostic ability of classification prediction was evaluated by obtaining specificity, sensitivity and the area under a receiver operating characteristic (ROC) curve (AUC).

Network of differentially expressed miRNAs and mRNAs

The pairwise Pearson correlation coefficients between key differentially expressed miRNAs and mRNAs were calculated. In total, six miRNA-target prediction tools, including miRWalk (version 2.0; http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/index.html), miRanda (http://www.microrna.org/), miRDB (version 2.0; http://mirdb.org/miRDB/), RNA22 (version 2.0; http://cm.jefferson.edu/rna22v2.0/), PICTAR2 (version 2.0; http://pictar.mdc-berlin.de/) and Targetscan (version 6.2; http://www.targetscan.org/) were used to predict the genes targeted by the differentially expressed miRNAs. Subsets of miRNA-target pairs with negative correlations were used to establish the regulatory network using Cytoscape software (version 3.3.0) (25). In addition, concrete ATCG base binding sites in the identified miRNA-mRNA pairs were also detected based on the Starbase database (http://starbase.sysu.edu.cn/index.php).

Functional analysis of the target mRNAs

To understand the biological function of the target mRNAs, Gene Ontology (GO) (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; version 2.2.0; http://www.genome.jp/kegg/pathway.html) pathway enrichment analyses were performed using Metascape (version 3.3.0) (http://metascape.org/gp/index.html). P<0.05 was considered to indicate a statistically significant difference.

In silico validation and in vitro validation

The Gene Expression Omnibus dataset GSE63046 (26) (involving 24 cases and 24 controls) was used to validate the expression of key differentially expressed miRNAs in tumor tissues compared with normal tissues from the same patients. The expression levels of these miRNAs are presented as box plots. Additionally, in vitro validation was performed by reverse transcription-quantitative PCR (RT-qPCR). Tumor and para-carcinoma tissues of seven patients were additionally collected for validation from December 30, 2018 to January 26, 2019 in The Third Hospital of Hebei Medical University. The clinical information (including therapy history, age and sex) of these patients was recorded before therapy.

The present study was approved by The Institutional Ethics Review Board of The Third Hospital of Hebei Medical University (approval no. 2018-025-1). In addition, informed consent was obtained from the individuals. Total RNA was extracted from tissue samples using TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc.), according to the manufacturer's protocols. A total of 2 µg RNA was used to synthesize cDNA using FastQuant Reverse Transcriptase (Sangon Biotech Co., Ltd.) for 60 min at 37°C followed by 5 min at 85°C. qPCR was performed in an ABI 7300 Real-time PCR system with SYBR® Green PCR Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycling conditions were as follows: Initial denaturation for 30 sec at 95°C followed by 40 cycles of 5 sec at 95°C and 30 sec at 60°C. All reactions were performed in triplicate. Hsa-U6 was used as the internal reference. The universal miRNA reverse primer is 5′-AACGAGACGACGACAGAC-3′. The sequences of forward primers for all of the miRNAs analyzed were as follows: 5′-GCAAATTCGTGAAGCGTTCCATA-3′ for Hsa-U6, 5′-UACCCUGUAGAACCGAAUUUGUG-3′ for hsa-miR-10b-5p, 5′-ACAGAUUCGAUUCUAGGGGAAU-3′ for hsa-miR-10b-3p, 5′-UCAAGUCACUAGUGGUUCCGUUUAG-3′ for hsa-miR-224-5p, 5′-UAUGGCACUGGUAGAAUUCACU-3′ for hsa-miR-183-5p, and 5′-UUUGGCAAUGGUAGAACUCACACU-3′ for hsa-miR-182-5p. The experiments were repeated three times. The relative gene expression levels were calculated as fold-changes using the 2−ΔΔCq method (27). The fold change was calculated as the enrichment between tumor tissue and para-carcinoma tissue.

In addition, according to the clinical information, patients with hepatocellular carcinoma were divided into two groups: i) Cirrhosis (75 cases); and ii) without cirrhosis (123 cases) to study whether liver cirrhosis may affect the expression of identified differentially expressed miRNAs.

Diagnosis and prognosis analysis of key differentially expressed miRNAs

ROC analysis was performed to assess the diagnostic value of key differentially expressed miRNAs. In addition, the survival package in R language (https://cran.r-project.org/web/packages/survival/index.html) was used to assess the prognostic value. The 5-year survival curves were plotted according to the clinical information and survival time.

Statistical analysis

All statistical analyses were performed using GraphPad Prism (version 8.0; GraphPad Software, Inc.). For the RT-qPCR experiments, one-way ANOVA, followed by Tukey's test to discriminate among the means, was used to assess statistical significance among two groups. For the box plots, the rank sum test was used to calculate the P-value. P<0.05 was considered to indicate a statistically significant difference. Data are presented as the mean ± SEM. All experiments were repeated independently at least three times.

Results

miRNA and mRNA expression pattern

First, principal component analyses for all miRNAs and mRNAs were performed. The present results demonstrated that these miRNAs and mRNAs were clearly separated according to the type of tissue, normal and tumor (Fig. S1). A total of 14 differentially expressed (all upregulated) miRNAs and 2,982 differentially expressed (1,989 upregulated and 993 downregulated) mRNAs were identified. The 14 differentially expressed miRNAs are presented in Table I. The heat maps corresponding to all miRNAs and top 50 mRNAs are presented in Figs. 1 and 2, respectively.

Figure 1. Heat map of all differentially expressed miRNAs in hepatocellular carcinoma. The diagram shows the result of a two-way hierarchical clustering of all differentially expressed miRNAs and samples. Clustering was analyzed using the complete-linkage method together with the Euclidean distance. Each row represents a differentially expressed miRNA and each column represents a sample. The differentially expressed miRNA color clustering tree is presented on the right. The color scale illustrates the relative expression level of differentially expressed miRNAs. Red indicates below the reference channel. Green indicates above the reference. miRNA, microRNA.

Figure 2. Heat map of the top 100 differentially expressed mRNAs in hepatocellular carcinoma. The diagram shows the result of a two-way hierarchical clustering of the top 100 differentially expressed mRNAs and samples. The clustering was established using the complete-linkage method together with the Euclidean distance. Each row represents a differentially expressed mRNA and each column represents a sample. The differentially expressed mRNA color clustering tree is presented on the right. The color scale illustrates the relative level of differentially expressed mRNA expression. Red indicates below the reference channel. Green indicates above the reference.

Table I. A total of 14 differentially expressed miRNAs in hepatocellular carcinoma, which were all upregulated.

Table I.

A total of 14 differentially expressed miRNAs in hepatocellular carcinoma, which were all upregulated.

miRNA	log2 Fold Change	P-value	FDR
hsa-miR-10b-5p	3.603003957	1.14×10−58	2.99×10−56
hsa-miR-224-5p	3.343475173	1.70×10−47	1.48×10−45
hsa-miR-183-5p	3.870865325	3.51×10−45	2.63×10−43
hsa-miR-1269a	5.556896527	4.26×10−42	2.79×10−40
hsa-miR-182-5p	3.38187178	4.08×10−38	2.14×10−36
hsa-miR-10b-3p	3.549127719	1.34×10−31	4.39×10−30
hsa-miR-96-5p	3.726215737	2.36×10−30	6.51×10−29
hsa-miR-217	4.036865911	2.25×10−27	4.71×10−26
hsa-miR-9-5p	3.161938401	2.98×10−26	5.78×10−25
hsa-miR-196b-5p	3.320156061	1.29×10−25	2.25×10−24
hsa-miR-135a-5p	4.141930478	8.47×10−19	9.87×10−18
hsa-miR-216b-5p	3.822995579	1.83×10−18	2.05×10−17
hsa-miR-216a-5p	3.427536419	7.89×10−17	7.38×10−16
hsa-miR-552-5p	3.736498894	6.81×10−14	4.46×10−13

[i] miRNA/miR, microRNA; FDR, false discovery rate.

Identification of optimal diagnostic biomarkers based on a machine learning approach

The RF feature selection and classification (DT, SVM and RF) procedures were performed for the identification of diagnostic biomarkers. All differentially expressed miRNAs were ranked according to the standardized drop in prediction accuracy (Fig. 3A). Differentially expressed miRNAs, including hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p, were considered as the optimal diagnostic biomarkers for hepatocellular carcinoma after subsequently adding one differentially expressed miRNA at a time in a top-down forward-wrapper approach (Fig. 3B). These five optimal differentially expressed miRNAs with diagnostic value for hepatocellular carcinoma were used to establish various classification models, including DT, SVM and RF. The AUC values in the RF, SVM and DT models were 98.2, 97 and 83.1%, respectively (Fig. 4). The RF model (with the largest AUC value) could effectively predict hepatocellular carcinoma.

Figure 3. Identification of optimal diagnostic miRNAs biomarkers for hepatocellular carcinoma. (A) Ranking of all differentially expressed miRNAs. Differentially expressed miRNAs were ranked according to the standardized drop in prediction accuracy. (B) Tendency chart of the area under the curve along with the increase of differentially expressed miRNAs. miRNA/miR, microRNA.

Figure 4. ROC results of the combination of five optimal diagnostic biomarkers, such as hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p and hsa-miR-183-5p and hsa-miR-224-5p, based on RF, SVM and DT classification models. ROC, receiver operating characteristic; RF, random forests; SVM, support vector machine; DT, decision tree; AUC, area under the curve.

Network of differentially expressed miRNAs and mRNAs

The correlation analysis between the five optimal differentially expressed miRNAs (hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p) and differentially expressed mRNAs was then performed. Following correlation analysis, 3,756 miRNA-mRNA pairs were identified to be negatively correlated (P<0.05; r<0). In the target prediction and negative correlation analyses, 170 miRNA-mRNA pairs, including five miRNAs (upregulated) and 145 mRNA (downregulated) were identified. The established regulatory network of miRNA-targeted mRNAs with negative correlation is presented in Fig. 5. Pairwise Pearson correlation analyses between the five optimal differentially expressed miRNAs (hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p) and their differentially expressed target mRNAs [secreted frizzled related protein 1 (SFRP1), endothelin receptor type B (EDNRB), nuclear receptor subfamily 4 group A member 3 (NR4A3), four and a half LIM domains 2 (FHL2), NK3 homeobox 1 (NKX3-1), interleukin 6 signal transducer (IL6ST) and forkhead box O1 (FOXO1)] are presented in Table II. The present results suggested that there was a significantly correlation between the identified miRNAs and mRNAs.

Figure 5. Network of miRNA-target mRNAs with a negative correlation between five differentially expressed miRNAs and 145 differentially expressed miRNAs in hepatocellular carcinoma. The diamonds and ellipses represent the differentially expressed miRNAs and target mRNAs, respectively. Pink and blue represent upregulated miRNA and downregulated mRNA, respectively. miRNA, microRNA.

Table II. Pairwise Pearson correlation analysis between five optimal differentially expressed miRNAs and their target differentially expressed mRNAs.

Table II.

Pairwise Pearson correlation analysis between five optimal differentially expressed miRNAs and their target differentially expressed mRNAs.

mRNA	miRNA	cor	P-value
SFRP1	hsa-miR-10b-5p	−0.27623	2.70×10−8
EDNRB	hsa-miR-10b-3p	−0.14088	0.005201
NR4A3	hsa-miR-224-5p	−0.21914	1.20×10−5
NKX3-1	hsa-miR-224-5p	−0.21888	1.23×10−5
FHL2	hsa-miR-224-5p	−0.15065	0.002787
IL6ST	hsa-miR-183-5p	−0.12462	0.013548
FOXO1	hsa-miR-182-5p	−0.20816	3.27×10−5

[i] miRNA/miR, microRNA; cor, correlation coefficient; SFRP1, secreted frizzled related protein 1; EDNRB, endothelin receptor type B; NR4A3, nuclear receptor subfamily 4 group A member 3; NKX3-1, NK3 homeobox 1; FHL2, four and a half LIM domains 2; IL6ST, interleukin 6 signal transducer; FOXO1, forkhead box O1.

The Starbase database was used to identify the binding sites between the five miRNA and their target mRNAs, and three ATCG base binding sites were found between hsa-miR-183-5p and IL6ST (Fig. 6A), two ATCG base binding sites between hsa-miR-224-5p and NR4A3 (Fig. 6B), two ATCG base binding sites between hsa-miR-224-5p and FHL2 (Fig. 6C), and two ATCG base binding sites between hsa-miR-182-5p and FOXO1 (Fig. 6D). However, there were no ATCG base binding sites between hsa-miR-224-5p and NKX3-1, hsa-miR-10b-5p and SFRP1, and hsa-miR-10b-3p and EDNRB.

Figure 6. Analysis of the concrete ATCG base binding sites in the miRNA-mRNA pairs identified. (A) ATCG base binding sites between hsa-miR-183-5p and IL6ST. (B) ATCG base binding sites between hsa-miR-224-5p and NR4A3. (C) ATCG base binding sites between hsa-miR-224-5p and FHL2. (D) ATCG base binding sites between hsa-miR-182-5p and FOXO1. miRNA/miR, microRNA; IL6ST, interleukin 6 signal transducer; NR4A3, nuclear receptor subfamily 4 group A member 3; FHL2, four and a half LIM domains 2; FOXO1, forkhead box O1.

Functional analysis of putative miRNA targets

To understand the potential function of the target differentially expressed mRNAs targeted by hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p, GO and KEGG pathway analyses were conducted. The present results suggested that these differentially expressed target mRNAs were most significantly enriched in the GO terms of ‘cellular response to lipid’, ‘fat cell differentiation’ and ‘monocarboxylic acid metabolic process’ (Fig. 7A). Additionally, ‘bile acid biosynthesis, cholesterol=>cholate/chenodeoxycholate’, ‘valine, leucine and isoleucine degradation’ and ‘calcium signaling pathway’ were the most enriched signaling pathways according to the KEGG analysis (Fig. 7B).

Figure 7. Significantly enriched GO and KEGG terms of differentially expressed mRNAs. (A) Top 20 significantly enriched GO terms of differentially expressed mRNAs. (B) Top 18 significantly enriched KEGG terms of differentially expressed mRNAs. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

In silico validation and in vitro validation

The GSE63046 dataset was used to validate the expression levels of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p (Fig. 8). The expression levels of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p were all upregulated in tumor tissues compared with normal tissues from the same patients. In addition, in vitro experiments were performed to further validate the expression level of the five differentially expressed miRNAs in seven patients. The clinical information of the seven patients enrolled in the present study is presented in Table III. The expression levels of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p were all upregulated (Fig. 9). The validation results were consistent with the present bioinformatics analysis.

Figure 8. Expression box plots of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p in the GSE63046 dataset. miR, microRNA.

Figure 9. In vitro reverse transcription-quantitative PCR validation of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p in tumor samples. Fold change>1 and fold change<1 represents upregulation and downregulation, respectively. miR, microRNA.

Table III. Clinical information of seven patients in the reverse transcription-quantitative PCR validation.

Table III.

Clinical information of seven patients in the reverse transcription-quantitative PCR validation.

Patient	Age	Sex	Family history	Other complications	HBV-infected	HCV-infected	Alpha fetoprotein content, 400 µg/l	Imaging examination results	Therapy history
1	55	Male	No	Liver cirrhosis	Yes	No	387.38	Primary hepatocellular carcinoma	Liver transplantation
2	71	Male	No	Liver cirrhosis	No	No	4.83	Primary hepatocellular carcinoma	Hepatectomy
3	74	Male	No	Liver cirrhosis	No	No	1.24	Primary hepatocellular carcinoma	Hepatectomy
4	56	Male	No	Liver cirrhosis	No	No	39.28	Primary hepatocellular carcinoma	Liver transplantation
5	65	Male	No	Liver cirrhosis	No	No	2.48	Primary hepatocellular carcinoma	Hepatectomy
6	50	Male	No	Liver cirrhosis	Yes	No	107.08	Primary hepatocellular carcinoma	Hepatectomy
7	54	Male	No	Liver cirrhosis	Yes	No	1.48	Primary hepatocellular carcinoma	Liver transplantation

[i] HBV, hepatitis B virus; HCV, hepatitis C virus.

In addition, according to the clinical information, patients with hepatocellular carcinoma were divided into two groups: i) Cirrhosis (75 cases); and ii) without cirrhosis (123 cases) to study whether liver cirrhosis may affect the expression of the five differentially expressed miRNAs identified (hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p) and their target mRNAs (SFRP1, EDNRB, NR4A3, FHL2, NKX3-1, IL6ST and FOXO1). The expression levels of these miRNAs and mRNAs are presented as box plots (Fig. 10). The present results suggested that there was no significant difference in these miRNAs and mRNAs, except for FHL2. FHL2 plays a protective mechanistic role during hepatic fibrogenesis (28). Moreover, deficiency in FHL2 aggravates liver fibrosis (28). Collectively, these results suggested that liver cirrhosis may affect special liver cirrhosis-related mRNAs, such as FHL2, without affecting the expression levels of the aforementioned miRNAs and mRNAs identified in patients with hepatocellular carcinoma.

Figure 10. Box plots of selected miRNAs and mRNAs. (A) Box plots of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p in liver cirrhosis. Box plots of selected miRNAs and mRNAs. (B) Box plots of SFRP1, EDNRB, NR4A3, FHL2, NKX3-1, IL6ST and FOXO1 in liver cirrhosis. *P<0.05 vs. cirrhosis. miR, microRNA; SFRP1, secreted frizzled related protein 1; EDNRB, endothelin receptor type B; NR4A3, nuclear receptor subfamily 4 group A member 3; FHL2, four and a half LIM domains 2; NKX3-1, NK3 homeobox 1; IL6ST, interleukin 6 signal transducer; FOXO1, forkhead box O1; ns, not significant.

Diagnosis and survival prediction of key differentially expressed miRNAs

ROC curve analysis was performed to assess the diagnosis ability of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p (Fig. 11A). The AUC values of hsa-miR-10b-5p (0.889), hsa-miR-10b-3p (0.871) and hsa-miR-224-5p (0.859) were all >0.8. For hepatocellular carcinoma diagnosis, the specificity and sensitivity of hsa-miR-10b-5p was 96.0 and 75.7%, respectively; the specificity and sensitivity of hsa-miR-10b-3p was 98.0 and 69.9%, respectively; and the specificity and sensitivity of hsa-miR-224-5p was 96.0 and 70.5%, respectively. In addition, the potential prognostic values of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p were analyzed using the online software survival package in R (https://cran.r-project.org/web/packages/survival/index.html). The present results demonstrated that hsa-miR-10b-5p and hsa-miR-10b-3p were considered to be significantly negatively associated with survival (P<0.05) in patients with hepatocellular carcinoma. The survival curves of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p are presented in Fig. 11B.

Figure 11. Diagnostic and survival analysis of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p. (A) ROC curves of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p between patients with hepatocellular carcinoma and healthy controls. ROC curves were used to show the diagnostic ability with 1-specificity and sensitivity. (B) Survival curves of hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-182-5p, hsa-miR-183-5p and hsa-miR-224-5p in patients with hepatocellular carcinoma. Survival curves were used to show the survival ability with time and survival rate. ROC, receiver operating characteristic; miR, microRNA; AUC, area under the curve.

Discussion

In the present study, hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p were considered as optimal diagnostic biomarkers for hepatocellular carcinoma based on machine learning approaches. hsa-miR-10b-5p was upregulated in advanced liver fibrosis (29). In addition, higher levels of expression of hsa-miR-10b-5p were identified in hepatocellular carcinoma cell lines (30). A previous study suggested that hsa-miR-10b-5p may be associated with the invasion and migration in hepatocellular carcinoma (31). In the present study, it was additionally identified that hsa-miR-10b-5p upregulated in hepatocellular carcinoma, which was consistent with a previous study (30). Notably, it was identified that hsa-miR-10b-5p exhibited a diagnostic and prognostic value in patients with hepatocellular carcinoma. In addition, SFRP1 was one of the target mRNAs of hsa-miR-10b-5p. SFRP1, a putative tumor suppressor protein, is decreased in hepatocellular carcinoma (32). A previous study suggested that SFRP1 expression may be downregulated by methylation levels, which may activate the Wnt signaling pathway, and increase cell growth and proliferation in hepatocellular carcinoma (33). Huang et al (34) identified that SFRP1 was a potential diagnostic biomarker for hepatocellular carcinoma. The present results suggested that hsa-miR-10b-5p may regulate cell growth and proliferation of hepatocellular carcinoma by targeting SFRP1.

A previous study demonstrated that hsa-miR-10b-3p plays an important role in tumor growth and metastasis (35). hsa-miR-10b-3p was upregulated in hepatocellular carcinoma (36). Moreover, the upregulation of hsa-miR-10b-3p was associated with the diagnosis and prognosis of hepatocellular carcinoma (36,37). Similarly, the present study identified that hsa-miR-10b-3p was upregulated in hepatocellular carcinoma. Furthermore, hsa-miR-10b-3p had a significant diagnostic and prognostic value for patients with hepatocellular carcinoma. The present results suggested the important role of hsa-miR-10b-3p in the development of hepatocellular carcinoma. In the target analysis, EDNRB was identified as a target mRNA of hsa-miR-10b-3p. EDNRB, a tumor suppressor, is very frequently methylated in hepatocellular carcinoma tissues (14). EDNRB was previously identified as a biomarker for hepatocellular carcinoma (38). This suggested that hsa-miR-10b-3p may serve a crucial role in the process of hepatocellular carcinoma by regulating EDNRB.

hsa-miR-224-5p wass upregulated in both tumor tissues and blood in patients with hepatocellular carcinoma (39–41). In addition, hsa-miR-224-5p was significantly associated with survival rate in patients with hepatocellular carcinoma (42). Similarly, the present study identified that hsa-miR-224-5p was upregulated in hepatocellular carcinoma tumor. In addition, NR4A3, FHL2 and NKX3-1 were three of the target mRNAs of hsa-miR-224-5p. NR4A3, a transcription factor, is a regulator of hepatoma cell and is associated with survival time in patients with hepatocellular carcinoma (43,44). FHL2 is an anti-proliferative- and metastasis-associated gene (45,46). It was identified that FHL2 plays a protective mechanistic role during hepatic fibrogenesis (28). Moreover, deficiency in FHL2 aggravates liver fibrosis (28). The expression level of FHL2 is downregulated in most patients with hepatocellular carcinoma (47). NKX3-1, a tumor suppressor, is associated with liver fibrosis (48). Aberrant methylation of NKX3-1 was observed in hepatocellular carcinoma (49). It was previously suggested that NKX3-1 is a potential predictor of patients with hepatocellular carcinoma recurrence (50). Investigating the regulation between hsa-miR-224-5p and NR4A3, FHL2 and NKX3-1 may provide insight for the understanding of the molecular mechanism underlying hepatocellular carcinoma.

hsa-miR-183-5p was upregulated in advanced liver fibrosis and hepatocellular carcinoma tissue (26,51). Leung et al (52) identified that high expression of hsa-miR-183-5p was significantly associated with invasion and metastasis, and may be a potential biomarker for the survival time of patients with hepatocellular carcinoma. In the present study, it was demonstrated that hsa-miR-183-5p was upregulated in hepatocellular carcinoma, in agreement with the aforementioned previous studies. Moreover, IL6ST was one of the target mRNAs of hsa-miR-183-5p. IL6ST is a gene involved in liver development (53). Changes in IL6ST were significantly associated with hepatotoxicity, including liver damage, inflammation and fibrosis (54). Alterations in IL6ST were frequently observed in hepatocellular adenoma (55). In addition, frequent upregulation and mutations in IL6ST were also detected in hepatocellular carcinoma (56). The present study suggested that hsa-miR-183-5p served a crucial role in the development of hepatocellular carcinoma by targeting IL6ST.

hsa-miR-182-5p was upregulated in hepatocellular carcinoma tissues and cell lines (51,57). hsa-miR-182-5p was associated with hepatocellular carcinoma metastasis, and could be a potential diagnostic and prognostic biomarker in patients with hepatocellular carcinoma (57,58). In addition, hsa-miR-182-5p may be a predictor of early recurrence in patients with hepatocellular carcinoma undergoing surgery (57). In the present study, the expression of hsa-miR-182-5p was increased in hepatocellular carcinoma. Additionally, FOXO1 was predicted to be regulated by hsa-miR-182-5p. FOXO1 is one of the most abundantly expressed genes in the liver, and regulates the expression of genes involved in cell cycle, metabolism and differentiation (59). Calvisi et al (60) demonstrated that the expression level of FOXO1 was downregulated in hepatocellular carcinoma. FOXO1 could reverse epithelial-interstitial transformation by inhibiting invasion and metastasis in hepatocellular carcinoma cells (61). FOXO1 was previously considered as a prognostic biomarker and potential target for hepatocellular carcinoma (62). The present findings suggested that the assocation between hsa-miR-182-5p and FOXO1 was associated with hepatocellular carcinoma.

According to the KEGG pathway analysis performed on the mRNAs targeted by the five key miRNAs, ‘bile acid biosynthesis, cholesterol’ was the most enriched signaling pathway. Bile acids are essential for protecting the liver from cholesterol. Bile acid metabolism is significantly regulated by enzymes involved in the liver. Bile acids could be highly toxic if accumulated in high concentrations in the liver. Additionally, a previous study demonstrated that bile acids are promoters of hepatocarcinogenesis (63–65). In addition, cholesterol intake is an independent risk factor for hepatocellular carcinoma (66–68).

Collectively, a number of differentially expressed miRNAs and mRNAs were identified in the present study. According to a machine learning approach, hsa-miR-10b-5p, hsa-miR-10b-3p, hsa-miR-224-5p, hsa-miR-183-5p and hsa-miR-182-5p were considered as optimal diagnostic biomarkers for hepatocellular carcinoma. The genes targeted and downregulated by these five miRNAs, including SFRP1, EDNRB, NR4A3, FHL2, NKX3-1, IL6ST and FOXO1, may be involved in hepatocellular carcinoma tumorigenesis. However, there are certain limitations in the present study. The sample size in the RT-qPCR experiments was small and larger numbers of tumor tissues are required for validating the data of the present study. The molecular mechanisms of the differentially expressed miRNAs and mRNAs identified in hepatocellular carcinoma tumorigenesis were not investigated. Additional in vitro experiments, such as cell culture and establishment of an animal model, are required to further investigate the potential mechanisms underlying the disease. Furthermore, the minimally invasive diagnostic methods for hepatocellular carcinoma were lacking and the potential use of miRNAs as blood/serum markers of hepatocellular carcinoma requires further examination. The present study may provide research basis for the diagnosis and treatment of hepatocellular carcinoma.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was funded by Provincial Outstanding Clinical Medicine Talents (Medical Leading Talents).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

XZ, JD, JC, YW, QG, QZ, WL, BL, ZC, LT and JZ analyzed and interpreted the data. XZ and CZ wrote, edited and revised the manuscript. CZ proposed the conception fo the study and designed the project. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The present study was approved by The Institutional Ethics Review Board of The Third Hospital of Hebei Medical University (approval no. 2018-025-1). In addition, informed consent was obtained from the individuals.

Patient consent for publication

Informed written consent was obtained from all subjects.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References