Comprehensive analysis of potential prognostic genes for the construction of a competing endogenous RNA regulatory network in hepatocellular carcinoma.

BACKGROUND: Hepatocellular carcinoma (HCC) is an extremely common malignant tumor with worldwide prevalence. The aim of this study was to identify potential prognostic genes and construct a competing endogenous RNA (ceRNA) regulatory network to explore the mechanisms underlying the development of HCC.
METHODS: Integrated analysis was used to identify potential prognostic genes in HCC with R software based on the GSE14520, GSE17548, GSE19665, GSE29721, GSE60502, and the Cancer Genome Atlas databases. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway-enrichment analyses were performed to explore the molecular mechanisms of potential prognostic genes. Differentially expressed miRNAs (DEMs) and lncRNAs (DELs) were screened based on the Cancer Genome Atlas database. An lncRNA-miRNA-mRNA ceRNA regulatory network was constructed based on information about interactions derived from the miRcode, TargetScan, miRTarBase, and miRDB databases.
RESULTS: A total of 152 potential prognostic genes were screened that were differentially expressed in HCC tissue and significantly associated with overall survival of HCC patients. There were 13 key potential prognostic genes in the ceRNA regulatory network: eleven upregulated genes (CCNB1, CEP55, CHEK1, EZH2, KPNA2, LRRC1, PBK, RRM2, SLC7A11, SUCO, and ZWINT) and two downregulated genes (ACSL1 and CDC37L1) whose expression might be regulated by eight DEMs and 61 DELs. Kaplan-Meier curve analysis showed that nine DELs (AL163952.1, AL359878.1, AP002478.1, C2orf48, C10orf91, CLLU1, CLRN1-AS1, ERVMER61-1, and WARS2-IT1) in the ceRNA regulatory network were significantly associated with HCC-patient prognoses.
CONCLUSION: This study identified potential prognostic genes and constructed an lncRNA- miRNA-mRNA ceRNA regulatory network of HCC, which not only has important clinical significance for early diagnoses but also provides effective targets for HCC treatments and could provide new insights for HCC-interventional strategies.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors worldwide.1 However, the early onset of HCC is difficult to identify, and tumors are often found in the later stages of the disease.2,3 Many studies have demonstrated that the occurrence of HCC is closely related to the dysregulation of multiple genes and key signaling pathways.4,5 Therefore, it is important to study the molecular mechanisms of HCC carcinogens to help us explore more effective tumor markers for early detection and key therapeutic targets of HCC.

MicroRNAs are predicted to control the activity of about 50% of all protein-coding genes, and deregulated miRNA expression has been implicated in a variety of cancers.6,7 Accumulating evidence suggests that lncRNAs play important roles in the regulation of cell proliferation, invasion, and angiogenesis through interactions with miRNAs and mRNAs.8–10 The competing endogenous RNA (ceRNA) hypothesis was proposed, which revealed how mRNAs and lncRNAs “talk” to each other using miRNA response elements as letters of a new language, and how ceRNA activity forms a large-scale regulatory network across the transcriptome.11 The interaction between lncRNAs and mRNAs through shared miRNAs represents a novel layer of gene regulation that plays important roles in the physiology and development of cancer.12 However, studies of lncRNA– miRNA–mRNA ceRNA regulatory networks in HCC based on multiple databases have rarely been reported.

In this study, we integrated expression data from multiple databases and screened out potential prognostic genes in HCC by bioinformatic analysis. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway-enrichment analysis was performed to help us understand the molecular mechanisms underlying the development of HCC. Subsequently, we constructed a ceRNA regulatory network to clarify potential relationships among lncRNAs, miRNAs, and mRNAs. Moreover, expression levels of key potential prognostic genes in the ceRNA regulatory network were validated using quantitative real-time (qRT) PCR analysis. In summary, our study not only identified potential molecular markers for the early diagnosis and prognosis of HCC but also provides new insights into HCC-interventional strategies.

Methods

Gene-expression data

Five gene-expression-profile groups (GSE14520, GSE17548, GSE19665, GSE29721, and GSE60502) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) from a total of 79 HCC-tissue and 77 adjacent normal liver-tissue samples. Information about the five groups of GEO datasets is shown in Table 1. Whole-transcriptome sequencing (WTS) data of 374 HCC-tissue and 50 adjacent normal liver-tissue samples and clinical data of 374 HCC patients were downloaded from the Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov).

Table 1

Details for GEO data

GEO	Platform	Sample	Normal	Tumor
GSE19665	GPL570	HCC	10	10
GSE17548	GPL570	HCC	20	17
GSE29721	GPL570	HCC	10	10
GSE14520	GPL571	HCC	21	22
GSE60502	GPL96	HCC	18	18

Abbreviation: GEO, Gene Expression Omnibus.

Screening for differentially expressed genes using integrated bioinformatic analysis

Data of the gene-expression profiles (GSE14520, GSE17548, GSE19665, GSE29721, and GSE60502) were calibrated, standardized, and log2-transformed after being processed using the R software package. Differential expression analysis between HCC and adjacent normal liver tissue was performed using R software with the Limma package. Log2 fold change (FC) >1 and P<0.05 were used as cutoff criteria. Subsequently, lists of differentially expressed genes (DEGs) obtained from the differential expression analysis of the five GEO-microarray data-set groups were integrated using R software with the RRA (robust rank aggregation) package, which is openly available on the comprehensive R network (http://cran.r-project.org) and based on the hypothesis that each gene is randomly ordered in each experiment and that there is a greater likelihood of differential gene expression when a gene is ranked high in all experiments. P<0.05 was used as the cutoff criterion.

RNA-sequencing data downloaded from the TCGA database was organized into gene-expression-matrix files. Protein-coding mRNA- and lncRNA-expression data were separated from the RNA-sequencing data, and the extracted data were used without further transformation. Differential mRNA expression was analyzed in HCC-tissue and adjacent normal liver-tissue samples with the R software with the Limma package. Log2FC >1 and P<0.05 were used as cutoff criteria.

Kaplan–Meier curve analysis

Kaplan–Meier curve analysis was performed on the DEGs in HCC-tissue samples compared with adjacent normal liver tissue to obtain potential prognostic genes that were significantly associated with HCC patient prognoses, based on the downloaded clinical data and the standardized gene-expression-matrix file produced by R with the survival package. P<0.05 was used as the cutoff criterion.

GO and KEGG pathway-enrichment analyses of potential prognostic genes

GO and KEGG pathway-enrichment analysis was performed using the DAVID online tool (http://david.abcc.ncifcrf.gov), which provides a comprehensive set of high-throughput functional gene analyses to help researchers understand biological characteristics to analyze potential prognostic genes at the functional level. P<0.05 was used as the cutoff criterion.

Identification of differentially expressed miRNAs and differentially expressed lncRNAs

To clarify the potential interactions among the differentially expressed lncRNAs (DELs), miRNAs (DEMs), and potential prognostic genes and to construct the lncRNA–miRNA– mRNA regulatory networks. DEM and DEL analysis in HCC tissue and adjacent normal liver tissue from the TCGA data sets was performed using R with Limma. Log2FC >2 and P<0.01 were used as the cutoff criteria.

Construction of competing endogenous RNA regulatory network

DEMs that might affect the function of potential prognostic genes were simultaneously searched using the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw), TargetScan (http://www.targetscan.org) and miRDB (http://www.mirdb.org) online databases. Then, the miRcode (http://www.mircode.org) database was used to screen DELs that might affect the function of DEMs. Potential prognostic genes that did not show target DEMs or DELs were deleted to confirm the reliability of the ceRNA regulatory network. The result was visualized using Cytoscape software.

Human specimens

HCC-tissue and adjacent normal liver-tissue sections were collected from the Department of General Surgery at Beijing Tongren Hospital, Capital Medical University. No patients received chemotherapy, radiotherapy, biological treatments, or other treatments before the surgeries were performed. Written informed consent was obtained from all patients, and was in accordance with the Declaration of Helsinki. The study protocol was approved by the institutional review board of Beijing Tongren Hospital, Capital Medical University.

Quantitative real-time polymerase chain reaction

Total RNA was extracted from tissue using Trizol (Thermo Fisher Scientific, Waltham, MA, USA) and cDNA synthesized using a reverse-transcription kit (Thermo Fisher Scientific). RT-PCR was repeated in triplicate using an RT-PCR instrument (Thermo Fisher Scientific) according to the manufacturer’s instructions. The relative expression of each target gene was calculated using the comparative Ct (2−ΔΔCt) method.

Results

Screening for differentially expressed genes from GEO and TCGA databases

A total of 854, 348, 2,120, 820, and 1,052 DEGs were screened from the GSE14520, GSE17548, GSE19665, GSE29721, and GSE60502 data sets, respectively. The differential expression of multiple genes from the two sets of data included in each of the five microarrays is shown in Figure S1, and cluster heat maps of the top-50 DEGs are shown in Figure S2. A total of 429 DEGs were screened using RRA analysis, including 150 upregulated and 279 downregulated genes. The heat map of the top-20 upregulated and downregulated genes is shown in Figure 1. Differential expression analyses of protein-coding mRNAs in HCC-tissue and adjacent normal liver-tissue samples from the TCGA data sets were performed using R with Limma package (Log2FC >1 and P<0.05). A total of 4,844 DEGs were identified, including 3,756 upregulated and 1,088 downregulated mRNAs. DEG volcano plots and cluster heat maps are shown in Figure 2.

Figure 1

Log2FC heat map of the top 20 upregulated and downregulated genes of the five expression-microarray groups downloaded from the GEO data set.

Notes: The differentially expressed gene (DEG) list of each expression microarray from the GEO data was integrated using R software with the RRA package. The abscissa is the GEO ID, and the ordinate is the gene name. Red represents upregulated gene expression, and green represents downregulated gene expression in hepatocellular carcinoma tissues compared with adjacent normal liver tissues. Numbers in the box represent log2FC values that resulted from the integrated analysis.

Abbreviations: FC, fold change; GEO, Gene Expression Omnibus; RRA, robust rank aggregation.

Figure 2

Hierarchical clustering heat map (A) and volcano plot (B) of differentially expressed genes screened based on TCGA HCC WTS data.

Notes: Protein-encoding mRNA-expression data were screened and analyzed using R software and the Limma package. Log2FC >1 and P<0.05 were used as cutoff criteria.

Abbreviations: FC, fold change; TCGA, the Cancer Genome Atlas; HCC, hepatocellular carcinoma; WTS, whole-transcriptome sequencing; FDR, false discovery rate.

Identification of potential prognostic genes

The Venn online tool (http://bioinformatics.psb.ugent.be/webtools/Venn) was used to screen 369 genes, which were differentially expressed in HCC tissue samples compared with adjacent normal liver tissue samples from both the TCGA database and the five GEO database groups. Kaplan–Meier curve analysis was performed on the DEGs from the downloaded clinical data and the standardized gene-expression-matrix files, which were derived from the TCGA data set using R software with the survival package (P<0.05). Finally, a total of 152 potential prognostic genes were screened out, consisting of 86 upregulated genes and 66 downregulated genes, which were differentially expressed in HCC tissue compared with adjacent normal liver tissue and significantly associated with the prognosis of HCC patients, as shown in Table 2.

Table 2

Potential prognostic genes in hepatocellular carcinoma

Expression	Gene names
Upregulated	AKR1B10, ANLN, ASPM, AURKA, BIRC5, BUB1B, CCNA2, CCNB1, CCNE2, CDC6, CDC20, CDCA3, CDCA8, CDK1, CDKN2C, CENPE, CENPF, CENPM, CEP55, CHEK1, COL15A1, DLGAP5, DNAJC6, DTL, E2F8, ECT2, ENAH, EZH2, FAM83D, FANCI, FEN1, FOXM1, G6PD, GINS1, GPSM2, HELLS, HJURP, HMMR, KIF11, KIF14, KIF18A, KIF18B, KIF20A, KIF2C, KIF4A, KNTC1, KPNA2, LRRC1, MAD2L1, MCM10, MCM2, MELK, MKI67, NCAPG, NCAPG2, NCAPH, NDC80, NEK2, NUF2, OIP5, PBK, PLVAP, PRC1, PRIM1, PTTG1, RACGAP1, RAD51AP1, RAD54B, RRM2, SHCBP1, SLC7A11, SPC25, SPDL1, SPP1, STIL, STMN1, SUCO, TK1, TOP2A, TPX2, TRIP13, TTK, UBE2S, UBE2T, ZIC2, ZWINT
Downregulated	ACSL1, ACSM3, ADAMTSL2, ADH1A, AFM, AKR1D1, ALDOB, ANXA10, ASPA, BHMT, C7, C8B, CDC37L1, CLEC1B, COLEC10, CPN2, CPS1, CRHBP, CYP2C8, CYP2C9, CYP3A4, CYP3A43, CYP4A11, CYP4F2, DNASE1L3, EPHX2, ETFDH, ETS2, F9, FAM149A, FTCD, GHR, GNMT, GRAMD1C, GYS2, HGFAC, HMGCS2, HPX, HRG, IGFALS, IL33, KBTBD11, KLKB1, LCAT, LECT2, MT1HL1, NDRG2, NPY1R, PBLD, PDE2A, PON1, PON3, PROZ, PZP, RDH16, SERPINA10, SERPINA4, SLC10A1, SLC16A4, SLC22A1, SLC38A4, SOCS2, SPP2, TAT, VIPR1, XDH

GO and KEGG pathway-enrichment analyses

The DAVID online tool was used to analyze functional and pathway-enrichment to investigate the function of potential prognostic genes. Significant results regarding the GO enrichment analyses of potential prognostic genes are shown in Figure S3. The GO analysis divided potential prognostic genes into three functional groups: molecular function, biological process, and cell composition. In the biological process group, potential prognostic genes were primarily enriched in cell division, mitotic nuclear division, sister-chromatid cohesion, and mitotic spindle organization. In the cell-composition group, potential prognostic genes were mainly enriched in mid-body formations, condensed chromosome kinetochores, and spindle poles. In the molecular function group, potential prognostic genes were mainly enriched in microtubule binding, microtubule-motor activity, and monooxygenase activity. KEGG signaling pathways of potential prognostic genes were mainly enriched in factors associated with the cell cycle, metabolic pathways, retinol metabolism, arachidonic acid metabolism, the p53-signaling pathway, oocyte meiosis, chemical carcinogenesis, drug metabolism, the CYP/CYP450 pathway, complement and coagulation cascades, linoleic acid metabolism, and HTLV1 infection, as shown in Table 3 and Figure S4.

Table 3

KEGG pathway analysis of potential prognostic genes associated with HCC

Pathway	Description	Genes	Gene count	P-value
hsa04110	Cell cycle	CCNB1, CCNE2, CDK1, CDC6, MAD2L1, CDKN2C, TTK, BUB1B, CHEK1, CDC20, MCM2, PTTG1, CCNA2	13	2.88×10−9
hsa01100	Metabolic pathways	XDH, CYP3A4, CYP2C9, CYP2C8, ALDOB, FTCD, EPHX2, ADH1A, CPS1, TAT, TK1, PRIM1, ACSM3, CYP4A11, ASPA, G6PD, ACSL1, HMGCS2, AKR1B10, RRM2, BHMT, PON1, CYP4F2, RDH16, AKR1D1, PON3	26	2.46×10−4
hsa00830	Retinol metabolism	CYP3A4, CYP4A11, CYP2C9, CYP2C8, ADH1A, RDH16	6	4.76×10−4
hsa00590	Arachidonic acid metabolism	CYP4A11, CYP2C9, CYP2C8, EPHX2, CYP4F2	5	0.003462332
hsa04115	p53-signaling pathway	CCNB1, CCNE2, CDK1, RRM2, CHEK1	5	0.0045805
hsa04114	Oocyte meiosis	CCNE2, CDK1, MAD2L1, AURKA, CDC20, PTTG1	6	0.004779011
hsa05204	Chemical carcinogenesis	CYP3A43, CYP3A4, CYP2C9, CYP2C8, ADH1A	5	0.008567382
hsa00982	Drug metabolism – CYP/CYP450	CYP3A4, CYP2C9, CYP2C8, ADH1A	4	0.031237741
hsa04610	Complement and coagulation cascades	C8B, C7, KLKB1, F9	4	0.032426334
hsa00591	Linoleic acid metabolism	CYP3A4, CYP2C9, CYP2C8	3	0.034443742
hsa05166	HTLV1 infection	MAD2L1, CDKN2C, ETS2, BUB1B, CHEK1, CDC20, PTTG1	7	0.044372737

Abbreviations: HCC, Hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of differentially expressed miRNAs and lncRNAs

A total of 126 DEMs were identified, including 123 upregulated and three downregulated miRNAs. Additionally, 1,059 DELs were identified, including 1,002 upregulated and 57 downregulated lncRNAs. Heat maps and volcano plots of the DEMs and DELs are shown in Figure 3.

Figure 3

Hierarchical clustering heat map and volcano plot of screened DEMs (A, B) and DELs (C, D) based on TCGA HCC WTS data.

Notes: miRNA- and lncRNA-expression data were downloaded and analyzed using R software and the Limma package. Log2FC >2 and P<0.01 were used as cutoff criteria.

Abbreviations: DEMs, differentially expressed miRNAs; DELs, differentially expressed lncRNAs; FC, fold change; TCGA, the Cancer Genome Atlas; HCC, hepatocellular carcinoma; WTS, whole-transcriptome sequencing.

We found a total of 13 key potential prognostic genes whose expression might be affected by the eight DEMs. Meanwhile, the functions of these eight DEMs were predicted to be regulated by 61 DELs. Interactions among DELs, DEMs, and the key potential prognostic genes are listed in Table 4. The network of gene nodes and interactions were visualized using the Cytoscape software, as shown in Figure 4. Kaplan–Meier curve analyses were performed on the miRNAs and lncRNAs in the ceRNA regulatory network based on the HCC-patient clinical data from the TCGA database. Nine key candidate lncRNAs (AL163952.1, AL359878.1, AP002478.1, C2orf48, C10orf91, CLLU1, CLRN1-AS1, ERVMER61-1, and WARS2-IT1) were significantly associated with HCC-patient prognoses, as shown in Figure 5. However, no DEMs were identified in the ceRNA regulatory network related to the postoperative survival of HCC patients.

Table 4

Interaction between genes in ceRNA regulatory network

mRNA	miRNA	lncRNA
ZWINT, LRRC1	hsa-miR506	AL033381.1, FAM87A, LINC00221, LINC00501, CLLU1, AL359878.1, AC006305.1, LINC00355, HOTAIR, LINC00200, SACS-AS1, ERVMER61-1, DLX6-AS1, BPESC1, LINC00316, HOTTIP, PVT1, RMST
PBK, KPNA2, SLC7A11, SUCO	hsa-miR373	C2orf48, AC009065.1, AP002478.1, C10orf91, LINC00221, TCL6, AC087392.1, AC061975.6, CLLU1, AL359878.1, SACS-AS1, DLX6-AS1, LINC00462, LINC00494, LINC00322, HOTTIP, PVT1, GRM5-AS1, LINC00485
ACSL1	hsa-miR205	PART1, CCDC13-AS1, AP002478.1, FAM87A, LINC00308, TCL6, DSCR8, AL512652.1, CLLU1, AL163952.1, LINC00488, ZNF385D-AS1, LINC00351, CCDC26, SACS-AS1, ERVMER61-1, BPESC1, MYCNOS, GPC6-AS1, CLRN1-AS1, AC012640.1, HOTTIP, AC011453.1, CRNDE, PVT1, GRM5-AS1, RMST, LINC00485
RRM2, KPNA2, SUCO	hsa-miR519d	C2orf48, AP002478.1, FAM87A, LINC00308, LINC00221, TCL6, AC087392.1, AC061975.6, AL512652.1, TDRG1, CLLU1, AL359878.1, AC006305.1, HOTAIR, LINC00200, DLX6-AS1, LINC00462, SRGAP3-AS4, HOTTIP, PVT1, GRM5-AS1, RMST, AC040173.1
CEP55, CHEK1, CDC37L1	hsa-miR424	PART1, C2orf48, CCDC13-AS1, AL033381.1, AP002478.1, FAM87A, TCL6, AC087392.1, CLLU1, AC006305.1, AC016773.1, LINC00473, WARS2-IT1, SFTA1P, LINC00355, LINC00200, LINC00160, DLX6-AS1, BPESC1, DSCR10, GPC6-AS1, CLRN1-AS1, HOTTIP, GDNF-AS1, PVT1, RMST, LINC00485
KPNA2, SLC7A11, SUCO	hsa-miR372	C2orf48, AC009065.1, AP002478.1, C10orf91, LINC00221, TCL6, AC087392.1, AC061975.6, CLLU1, AL359878.1, SACS-AS1, DLX6-AS1, LINC00462, LINC00494, LINC00322, HOTTIP, PVT1, GRM5-AS1, LINC00485
EZH2	hsa-miR217	AL357153.1, AC024563.1, LINC00221, TCL6, CLLU1, AC016773.1, HOTAIR, LINC00200, MIR137HG, MYCNOS, LINC00494, CLRN1-AS1, CRNDE, PVT1, AC040173.1, NOVA1-AS1
CCNB1, SUCO	hsa-miR183	C2orf48, AL357153.1, AC024563.1, TCL6, LINC00501, LINC00392, LINC00200, MYCNOS, CRNDE, PVT1, AC040173.1

Abbreviation: ceRNA, competing endogenous RNA.

Figure 4

The ceRNA regulatory network in HCC.

Notes: The oval represents mRNA, the diamond represents lncRNA, and the square represents miRNA. Red indicates the gene was upregulated and green indicates the gene downregulated in HCC tissues compared with adjacent normal liver tissues. The line between the genes indicates that there could be regulatory relationships between the two genes.

Abbreviations: ceRNA, the competing endogenous RNA; HCC, hepatocellular carcinoma.

Figure 5

Kaplan–Meier curve analysis was performed on the DELs in the ceRNA regulatory network. P<0.05 was used as the cutoff criterion.

Abbreviations: DELs, differentially expressed lncRNAs; ceRNA, competing endogenous RNA.

Validation of key potential prognostic genes by quantitative real-time polymerase chain reaction

The key potential prognostic genes were subjected to qRT-PCR on 20 pairs of tissue samples. As shown in Figure 6, expression levels of eleven genes (CCNB1, CEP55, CHEK1, EZH2, KPNA2, LRRC1, PBK, RRM2, SLC7A11, SUCO, and ZWINT) were upregulated and those of two mRNAs (ACSL1, CDC37L1) downregulated in HCC tissue compared with adjacent normal liver tissue. (P<0.05). Dysregulation trends were consistent with the results of our previous bioinformatic analysis.

Figure 6

mRNA-expression levels in the ceRNA regulatory network of the 20 pairs of HCC tissue and adjacent normal liver tissue.

Notes: Experiments were repeated three times. *P<0.05.

Abbreviations: ceRNA, competing endogenous RNA; HCC, hepatocellular carcinoma.

Discussion

HCC is one of the leading causes of death worldwide, and is an overwhelming experience for patients, due to its poor prognosis.13 Therefore, exploring the underlying mechanisms and identifying effective prognostic biomarkers for the early detection of and effective therapy for HCC patients has become an urgent mandate.

Several studies have revealed that some dysregulated genes might act as potential biomarkers in HCC. However, previous studies mainly focused on different gene data without clinical data to predict prognosis in HCC.4,14 Several previous studies containing prognostic analyses were mostly based on a single database.15,16 In the present study, we screened and integrated the DEGs of five groups of gene-expression data from the GEO database, and then validated these DEGs based on gene expression and clinical survival data from the TCGA database, which could determine more reliable and accurate potential prognostic genes in HCC. This finding provided numerous potential biomarkers which may be used in diagnosing HCC and predicting patient prognoses with HCC. Considering that these potential prognostic genes are the implementers of molecular function, GO and KEGG pathway-enrichment analyses were performed to help us understand the molecular mechanisms underlying the development of HCC. Remarkably, our study reveals that the cell cycle was the most affected pathway in HCC, which may guide individualized therapy choices.

Currently, cross talk is thought to occur between lncRNAs, miRNAs, and mRNAs during many aspects of tumorigenesis, and the ceRNA hypothesis has been considered a novel measure of gene regulation at the posttranscriptional level in multiple tumor types.17–19 This study attempted to construct a lncRNA–miRNA–mRNA ceRNA regulatory network during HCC tumorigenesis at the transcriptional level using comprehensive analysis. To confirm the reliability of the ceRNA regulatory network, differentially expressed miRNAs and lncRNAs were selected, and potential prognostic genes that did not show target DEMs and DELs were deleted.

In this study, 13 key potential prognostic genes were included to develop the novel ceRNA network, the expression of which might be regulated by eight DEMs and 61 DELs. We conducted qRT-PCR experiments on 20 pairs of HCC-and adjacent liver-tissue samples, and the result showed that expression levels of the 13 key potential prognostic genes were consistent with the result of our previous analysis. Among these key potential prognostic genes, CCNB1,20 CEP55,21 CHEK1,22 EZH2,23 KPNA2,24 LRRC1,25 PBK,26 RRM2,27 SLC7A11,28 and ZWINT29 have been validated as oncogenes in HCC according to previous studies. Interestingly, no studies have reported a relationship between SUCO gene expression and tumorigenesis, and we are planning to conduct experiments to explore the role of SUCO in the development and progression of HCC. It has been reported that ACSL1 was upregulated by HBXIP in breast cancer tissue and might play crucial roles in breast cancer development.30 However, studies regarding the relationship between ACSL1 gene expression and HCC were not reported. In this study, we found that ACSL1 was downregulated in HCC tissues and that patients with high ACSL1 gene expression had a better prognosis, which needs further verification. Another study showed that patients with hepatitis B-associated HCC had higher CDC37L1 gene expression compared with patients without HCC, and CDC37L1 was able to protect HCC patients from persistent drug toxicity.31 In contrast, we found that the CDC37L1 gene was downregulated in HCC tissue compared with adjacent normal liver tissue based on bioinformatic analyses and qRT-PCR experiments.

Recently, dysregulated miRNAs were reported to be associated with many molecular pathways in tumor evolution and progression.32–34 In this ceRNA network, eight DEMs were identified as core elements by which mRNAs and lncRNAs can “talk” to each other in the physiology and development of HCC. The identified DEMs were miR183,35 miR205,36 miR217,37 miR372,38 miR373,39 miR424,40 miR506,41 and miR519d,42 which have been validated previously to play specific roles in HCC tumorigenesis. These DEMs, which were experimentally validated in previous studies, could ensure ceRNA-network stability and were more conducive to the identification of DELs in this study.

Plenty of studies have certified that lncRNAs could be involved in many functional aspects associated with tumor evolution and progression.43–45 In this novel ceRNA network, nine potential prognostic lncRNAs were identified through survival analyses that could participate in the development of HCC. CLLU1 expression was shown to be upregulated in chronic lymphocytic leukemia and has been used for prognostic and diagnostic purposes.46,47 It has been reported that AP002478.1 can function as a prognostic biomarker in patients with Helicobacter pylori-positive gastric cancer.48 Also, C2orf48 and C10orf91 are related to the overall survival time of patients with oral squamous-cell carcinoma.49 However, lncRNAs CLLU1, AP002478.1, C2orf48, and C10orf91 have not previously been reported in other tumor types, including HCC. Moreover, little is known about the other five potential prognostic lncRNAs (AL163952.1, AL359878.1, CLRN1-AS1, ERVMER61-1, and WARS2-IT1). Therefore, further investigations are needed to clarify the function of these potential prognostic lncRNAs in HCC and other cancers.

Several limitations of the current study should be discussed. First, no DEMs related to postoperative HCC-patient survival were identified in the novel ceRNA network. Second, we lacked the experimental validation of lncRNA and miRNA expression in tissue samples to confirm our findings. Therefore, further research is needed to verify our study findings and unravel the mechanism of ceRNA in HCC.

Conclusion

This study provides numerous potential prognostic genes, and the function and signaling pathways they participate in are presented systematically. Moreover, we constructed a lncRNA–miRNA–mRNA ceRNA regulatory network in HCC, which has important clinical significance for the early detection and diagnosis of HCC and for the potential development of therapeutic targets for HCC. However, further biological experiments are needed to assess the function of these potential prognostic genes and the effect of the miRNAs and lncRNAs present within the ceRNA regulatory network.

Differential gene expression between HCC tissue and adjacent normal liver tissue from five expression-microarray groups.

Notes: (A) GSE14520 data; (B) GSE17548 data; (C) GSE19665 data; (D) GSE29721 data; (E) GSE60502 data. The red points represent genes that were upregulated and the green points those that were downregulated in HCC tissue compared with adjacent normal liver tissue, which were analyzed using R software and the Limma package. The black points represent genes without a significant difference between the two tissue types. Log2FC >1 and P<0.05 were used as cutoff criteria. Microarray groups downloaded from the GEO data set.

Abbreviations: FC, fold change; GEO, Gene Expression Omnibus; HCC, hepatocellular carcinoma.

Hierarchical clustering heat maps of differentially expressed genes of five expression-microarray groups, downloaded from the GEO data set.

Notes: (A) GSE14520 data; (B) GSE17548 data; (C) GSE19665 data; (D) GSE29721 data; and (E) GSE60502 data. The color of genes gradually changes with expression-level changes in hepatocellular carcinoma (HCC) tissue compared with adjacent normal liver tissue. Black indicates no significant changes in gene expression. Prominent gene expression in HCC tissue is either upregulated (red color) or downregulated (green color). The depth of color corresponds to expression intensity. Data were analyzed using the R software and the Limma package using. Log2FC >1 and P<0.05 were used as cutoff criteria.

Abbreviations: FC, fold change; GEO, Gene Expression Omnibus.

GO enrichment analysis of the potential prognostic genes in hepatocellular carcinoma.

Notes: (A) The GO analysis divided potential prognostic genes into three functional groups: the biological process, cell composition, and molecular function groups. (B) Significant GO enrichment items relating to the potential prognostic genes in the different functional groups.

Abbreviation: GO, Gene Ontology.

Pathway enrichment of potential prognostic genes.

Note: Blue represents signaling pathways, red represents upregulated genes, and green represents downregulated genes.